5 Session 5: Many variables at once: PCA, heatmaps and PLS-DA

5.1 Overview

Quick poll (Poll Everywhere)

Biological question

Drought and the heatwave reduced the plants’ biomass (Session 4). Do they also change the plants’ chemistry, the same way in the root, the leaf, and in what the roots release into the soil?

In Session 4, you analysed one variable at a time: the biomass, then the Shannon diversity. Today, each sample comes with dozens, then thousands of variables at once: the concentrations of 15 to 37 small molecules (metabolites) measured in the root, the stem, the leaf and the root exudates of the same soybean plants, then 2,500 untargeted signals from the exudates. One boxplot per molecule would give you dozens of plots and no overview. You will learn the classic tools to see the whole table at once:

PCA and heatmaps with clustering, which look for structure without being told which sample is which (unsupervised);
PLS-DA, which is told the groups and looks for what separates them (supervised), and the trap that comes with it;
the volcano plot: one test per molecule, as in Session 4, and why you must correct the p-values when you run many tests.

How you work and submit today

As in Session 4, everything happens in WORKBOOK.qmd and you click Render. New today: Render makes a PDF, WORKBOOK.pdf, and that is what you submit. Blocks A to D are guided; blocks E (the leaf) and F (the exudates) you do on your own. Every block counts in the grade. Submit on Moodle by Sunday evening (25 October).

5.1.1 Learning outcomes

By the end of this session, you should be able to:

explain why concentrations are log-transformed and scaled before a multivariate analysis;
turn a long table into a sample × variable matrix with pivot_wider();
run a PCA, read the share of variance of each axis, the scores (samples) and the loadings (variables);
draw a heatmap with clustering, and read its trees;
run a PLS-DA, and explain why a supervised method must always be cross-validated;
run one test per variable, correct the p-values for multiple testing, and draw a volcano plot;
say which method answers which question.

5.1.2 Session plan

Time	Block	What happens
14:15-14:20	Retrieval	Session 4 quiz.
14:20-14:35	Theory 1	Metabolomics: targeted and untargeted, and the data you get.
14:35-14:50	Theory 2	From one variable to many: PCA, clustering, PLS-DA, volcano.
14:50-15:00	Practical briefing	What today’s workbook will ask of you.
15:00-15:15	Practical A	Install mixOmics (once); same project habits as Session 4, new data.
15:15-15:35	Practical B	One molecule at a time: ABA, proline, and why we take the log.
15:35-15:45	Break
15:45-16:00	Warm-up	Guided: what a PCA does, on decathlon athletes.
16:00-16:35	Practical C	Guided: the whole root metabolome, without labels (PCA, heatmap).
16:35-17:10	Practical D	Guided: PLS-DA and its trap, then one test per molecule.
17:10-17:15	Wrap-up	Render, and how blocks E (the leaf) and F (the exudates) work.

5.2 Theory (≈45 min)

5.2.1 1. From one variable to many

Each sample now has one value per metabolite. For the analyses of today, the table must have one row per sample and one column per metabolite: a matrix. Your data file is long (one row per sample and metabolite, as in Session 4), so the first step is pivot_wider(), the reverse of the pivot_longer() you used in Session 4.

Two more steps before any multivariate method:

log-transform: concentrations go from 0.02 to about 190,000. On the raw scale, a few very abundant molecules would dominate everything. log10() turns ratios into differences: a molecule that doubles moves by the same amount, whether it is rare or abundant.
scale: each metabolite is centred and divided by its standard deviation, so that each one has the same weight in the analysis (scale. = TRUE in prcomp()).

5.2.2 2. Unsupervised: PCA and clustering

A PCA (principal component analysis) summarises the 34 columns into a few new axes, the principal components. PC1 is the direction along which the samples differ the most, PC2 the next one, at right angles, and so on. Three things to read:

the share of variance of each axis: how much of the differences between samples it captures;
the scores: where each sample falls on the axes (the usual PCA plot);
the loadings: how much each metabolite contributes to each axis, and in which direction.

A heatmap with clustering shows the whole matrix as colours, and reorders rows and columns so that similar ones sit next to each other. The trees (dendrograms) show the groups found.

In the practical, you first run a PCA on data you know without any biology: 41 decathlon performances, 10 events each. You will see the axes rediscover the best athletes, and why scaling matters when the variables are in different units (seconds and metres).

Both are unsupervised: they never see your labels (water, heat). If the samples still group by condition, the condition is a major source of variation in the data.

5.2.3 3. Supervised: PLS-DA, and why it must be validated

A PLS-DA (partial least squares discriminant analysis) is told the groups, and looks for the combination of metabolites that best separates them. It is useful to find which molecules distinguish two conditions, but it is dangerous: with many variables and few samples, it can separate any groups, even groups drawn at random. A separation on the data the model learned from proves nothing.

The answer is cross-validation: hide part of the samples, train the model on the rest, predict the hidden samples, repeat. Only the accuracy on samples the model did not see tells you whether the separation is real. You will see this trap with your own eyes in block D.

5.2.4 4. Many tests at once: the volcano plot

You can also test each metabolite on its own: one t-test per molecule, as in Session 4. With 34 tests at the 5% level, you expect about 1 or 2 “significant” results by chance alone. Correcting the p-values (here with the Benjamini-Hochberg method, which controls the false discovery rate) keeps that under control.

A volcano plot shows all the tests at once: the effect on the x axis (log2 fold change: positive means higher under drought), the evidence on the y axis (−log10 of the adjusted p-value: the higher, the stronger).

Which method, when?

Your question	Method	In R	Supervised?
What are the main patterns in the samples? Do they group by condition?	PCA	`prcomp()`	no
Which samples and which variables look alike?	heatmap with clustering	`pheatmap()`	no
Which combination of variables separates known groups?	PLS-DA, always with cross-validation	`plsda()`, `perf()`	yes
Which variables change between two groups, one by one?	one test per variable + FDR, volcano plot	`t.test()`, `p.adjust()`	(it is a test)

5.2.5 5. Worked example: the metabolome of the soybean plants

Today’s data come from the plants you met in Session 4, measured with two kinds of metabolomics:

	Targeted	Untargeted
What is measured	a chosen list of known molecules	every signal the machine detects
Each variable is	a named molecule, quantified	a “feature”: a peak with a retention time and a mass (m/z)
How many variables today	15 to 37 per organ	2,500 in the exudates
Identity	known for sure	unknown for many; a putative name from a database search, with a confidence score
Good for	testing precise hypotheses, comparing known pathways	discovering what you did not think of looking for

Root, stem, leaf (targeted): the omics_and_soil plants, frozen at harvest and ground; 4 samples per cultivar and condition, each a pool of two plants.
Exudates (targeted and untargeted): what the roots of the exudates plants released into a solution, collected at harvest; one sample per plant.

The exudates are the molecules the root gives to the soil: the food and the signals of the bacteria you studied in Session 4.

5.2.6 Practical briefing

Blocks A-D in pairs, blocks E and F on your own. Block B looks at single molecules you can understand without any statistics. The warm-up runs a first PCA on data you can check with common sense: decathlon athletes and their 10 events. Blocks C and D analyse the root, guided: first without labels (PCA, heatmap), then with labels (PLS-DA, volcano). In block E, you analyse the leaf: same reasoning, a different answer. In block F, the exudates, first with the targeted list, then with the untargeted scan: same plants, two very different tables. Your submission is the whole rendered workbook: every block counts.

5.3 Your R toolbox

New tools today, on top of the Session 4 toolbox (still valid: see the Session 4 R toolbox).

pivot_wider() is the reverse of pivot_longer(): it spreads one column into many.

root_wide <- root_long %>%
  pivot_wider(names_from = metabolite, values_from = log_conc)

names_from: the column whose values become the new column names (here: the metabolite names);
values_from: the column that fills the new columns (here: the log concentrations).

Then select(-sample, ...) and as.matrix() keep the numbers only, as a matrix, and rownames() puts the sample names back on the rows.

Common mistake. Keeping a column that differs inside a sample (such as the raw concentration next to log_conc): pivot_wider() then cannot put one sample on one row. Select only the columns you need before pivoting.

prcomp() is part of R itself: no package needed.

Code	What it gives
`prcomp(matrix, scale. = TRUE)`	the PCA, each variable scaled
`summary(pca)$importance[2, ]`	the share of variance of each axis
`pca$x`	the scores: one row per sample, one column per axis
`pca$rotation`	the loadings: one row per variable, one column per axis

Reading the loadings. A metabolite with a large negative loading on PC1 is high in the samples on the left of the PCA plot; a large positive loading, high on the right.

Common mistake. Forgetting scale. = TRUE: the variables with the largest numbers then decide everything.

pheatmap(t(root_matrix), scale = "row", annotation_col = annotation,
         show_colnames = FALSE, fontsize_row = 7)

t() turns the matrix around: metabolites in rows, samples in columns;
scale = "row": each metabolite is shown as “above or below its own average” (z-score), otherwise the most abundant ones would fill the colour scale;
annotation_col: a small table, one row per sample (with the sample names as row names), whose columns become coloured bars on top;
the trees on the left and on top are the clustering: they are computed from the data, not from the annotation.

The reference package for PLS-DA (and much more) in biology.

Code	What it does
`plsda(matrix, groups, ncomp = 2)`	fits a PLS-DA with 2 components
`plotIndiv(model, ellipse = TRUE, legend = TRUE)`	the sample plot, with one ellipse per group
`predict(model, matrix)$class$max.dist`	the predicted group of each sample
`perf(model, validation = "Mfold", folds = 4, nrepeat = 10)`	cross-validation: 4 folds, repeated 10 times
`perf(...)$error.rate$BER`	the balanced error rate, per component (accuracy = 1 − error)

Install once, load every time. BiocManager::install("mixOmics", update = FALSE) runs once, in the Console (block A0). library(mixOmics) goes at the top of every report. Never put an install line in your workbook: every Render would install it again.

Common mistake: load it first. mixOmics loads another package (MASS) that also has a function called select(). Whatever is loaded last wins. So load mixOmics before the tidyverse; otherwise select() stops working, with a confusing error (“unused argument”).

Documentation. mixomics.org

geom_text_repel() writes text labels on a plot and moves them so that they do not overlap, unlike geom_text().

geom_text_repel(data = slice_min(root_tests, p_adjusted, n = 10),
                aes(label = metabolite), size = 3)

data = ... labels only some of the points: here, the 10 with the smallest adjusted p-values (slice_min(), from dplyr).

5.4 Practical (≈1 h 50)

Tested on LEGcompute

pheatmap and ggrepel are already installed in RStudio Server’s R. mixOmics is not: each of you installs it once, at the start of block A, in about 1 to 2 minutes (tested: 83 seconds, mixOmics 6.36.0).

5.4.1 A. Install mixOmics, same habits, new data (15:00-15:15)

A0. Install mixOmics, once. In the Console (not in the workbook), run these three lines:

dir.create(Sys.getenv("R_LIBS_USER"), recursive = TRUE, showWarnings = FALSE)
.libPaths(c(Sys.getenv("R_LIBS_USER"), .libPaths()))
BiocManager::install("mixOmics", update = FALSE, ask = FALSE)

line 1 creates your personal R library, a folder in your home where you can install packages (you cannot write in the shared one);
line 2 tells R to use it, from now on in this session;
line 3 installs mixOmics from Bioconductor (a repository of biology packages), with the few packages it still needs. update = FALSE is important: otherwise R offers to update dozens of packages you already have.

It takes 1 to 2 minutes, and prints a lot of compilation messages: that is normal. While it runs, go on with the Terminal commands below.

Checkpoint

The Console ends with * DONE (mixOmics). Then library(mixOmics) in the Console loads it without error (a few welcome lines are printed). You will never need to install it again: from now on, R finds it in your personal library.

A1. The project. In the Terminal tab:

mkdir -p ~/bioinfo_course/session_05/data ~/bioinfo_course/session_05/img ~/bioinfo_course/session_05/results
cp /home/cmaslard/bioinfo_tools_course/session_05/soybean_metabolites.csv ~/bioinfo_course/session_05/data/
cp /home/cmaslard/bioinfo_tools_course/session_05/soybean_exudates_untargeted.csv ~/bioinfo_course/session_05/data/
cp /home/cmaslard/bioinfo_tools_course/session_05/soybean_exudates_annotation.csv ~/bioinfo_course/session_05/data/
cp /home/cmaslard/bioinfo_tools_course/session_05/decathlon.csv ~/bioinfo_course/session_05/data/
cp /home/cmaslard/bioinfo_tools_course/session_05/experimental_design.png ~/bioinfo_course/session_05/img/
cp /home/cmaslard/bioinfo_tools_course/session_05/WORKBOOK.qmd ~/bioinfo_course/session_05/
ls -R ~/bioinfo_course/session_05

Then File → New Project → Existing Directory, pick ~/bioinfo_course/session_05, open WORKBOOK.qmd, write your name at the top and click Render once.

Checkpoint

WORKBOOK.pdf appears (RStudio opens it), with the design figure as Figure 1, and the first chunk lists the columns sample, organ, genotype, water, heat, panel, metabolite, kegg_id, concentration. A few chunks further down show Error: unexpected input: those are the ___ you will fill.

panel is the analytical panel that measured each compound (amino acids, organic acids, flavonoids, phytohormones, cytokinins), as in the lab’s files: a practical grouping, not always the strict chemical family. concentration was cleaned before class: values below the detection limit were replaced by half of the smallest value measured for that molecule, a common rule in metabolomics.

5.4.2 B. One molecule at a time (15:10-15:35)

B1. What is in the table?

# how many samples per organ? (count by the column that holds the organ)
metabolites %>% distinct(organ, sample) %>% count(___)
metabolites %>% distinct(organ, metabolite) %>% count(organ)

Checkpoint

Two small tables, each with one row per organ: exudate, leaf, root, stem. If you see one row per sample instead, check what you put in count().

B2. The drought hormone. Abscisic acid (ABA) makes the leaves close their stomata, the small pores through which they lose water. Look at it in the three organs (it is not on the exudate list):

aba <- metabolites %>%
  filter(metabolite == "Abscisic Acid")

# one panel per organ
ggplot(aba, aes(x = water, y = concentration, colour = heat)) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(position = position_jitterdodge(jitter.width = 0.1)) +
  facet_wrap(~ ___, scales = "free_y")

# the mean per organ and per water regime
aba %>%
  group_by(organ, ___) %>%
  summarise(mean_aba = mean(concentration))

scales = "free_y" gives each organ its own y axis: the levels are very different between organs.

B3. Proline, a small amino acid that protects cells when water is short. Same plot, with one ___ to fill.

B4. Why the log?

range(metabolites$___)

ggplot(metabolites, aes(x = concentration)) +
  geom_histogram(bins = 50)

# the same histogram, on the log10 scale
ggplot(metabolites, aes(x = log10(___))) +
  geom_histogram(bins = 50)

Checkpoint

range() prints two numbers. The first histogram is one tall bar near zero; the second one is spread out and readable.

5.4.3 15:35-15:45: Short break

5.4.4 Warm-up: what a PCA does, on athletes (15:45-16:00, pairs)

Before 34 metabolites, a PCA on data you can check with common sense: 41 decathlon performances (Decastar meeting and Athens Olympic Games, 2004; data of the FactoMineR package). Each athlete did 10 events: four runs (times in seconds, small is good) and six jumps and throws (metres, large is good).

W1. The table.

decathlon <- read_csv("data/decathlon.csv")
glimpse(decathlon)

# how many athletes per competition?
decathlon %>% count(___)

W2. Without and with scaling. The same PCA twice: first with each event in its own unit, then with every event scaled.

# the 10 events only, as a matrix, with the athletes' names as row names
events <- decathlon %>%
  select(run_100m:run_1500m) %>%
  as.matrix()
rownames(events) <- decathlon$___

# how much does each event vary between athletes? (standard deviation, in its own unit)
decathlon %>% summarise(across(run_100m:run_1500m, sd))

# a PCA WITHOUT scaling: each event keeps its own unit (seconds, metres)
pca_raw <- prcomp(events, scale. = FALSE)
round(pca_raw$rotation[, 1], 2)   # the loadings on PC1

# the same PCA WITH scaling: every event gets the same weight
decathlon_pca <- prcomp(events, scale. = ___)
decathlon_variance <- round(100 * summary(decathlon_pca)$importance[2, ], 1)
decathlon_variance

across(run_100m:run_1500m, sd) applies sd() to every column from run_100m to run_1500m.

W3. Arrows and points. The loadings become arrows (one per event), the scores points (one per athlete):

# the events: one arrow per event, from 0 to its loadings on PC1 and PC2
decathlon_loadings <- as_tibble(decathlon_pca$rotation[, 1:2], rownames = "event")

ggplot(decathlon_loadings) +
  geom_segment(aes(x = 0, y = 0, xend = ___, yend = ___),
               arrow = arrow(length = unit(0.2, "cm"))) +
  geom_text_repel(aes(x = PC1, y = PC2, label = ___)) +
  coord_equal() +
  labs(x = paste0("PC1 (", decathlon_variance[1], "%)"),
       y = paste0("PC2 (", decathlon_variance[2], "%)"))

# the athletes: one point per athlete, coloured by the final score
decathlon_scores <- as_tibble(decathlon_pca$x[, 1:2]) %>%
  bind_cols(select(decathlon, athlete, competition, points))

ggplot(decathlon_scores, aes(x = PC1, y = PC2, colour = ___)) +
  geom_point(size = 3) +
  geom_text_repel(aes(label = athlete), size = 3) +
  labs(x = paste0("PC1 (", decathlon_variance[1], "%)"),
       y = paste0("PC2 (", decathlon_variance[___], "%)"))

# does PC1 follow the final score?
cor(decathlon_scores$PC1, decathlon_scores$___)

Checkpoint

Two plots: one with 10 arrows (one per event), one with 41 points. Some names appear twice, because 9 athletes competed in both meetings. On another computer, the whole plot may be flipped left-right or up-down: the direction of an axis is arbitrary, only the angles between arrows matter.

5.4.5 C. Guided: the whole root metabolome, without labels (16:00-16:35, pairs)

C1. From long to wide. One row per sample, one column per metabolite, log10 of the concentration:

# one row per sample, one column per metabolite, log10 of the concentration
root_wide <- metabolites %>%
  filter(organ == ___) %>%
  mutate(log_conc = ___(concentration)) %>%
  select(sample, genotype, water, heat, metabolite, log_conc) %>%
  pivot_wider(names_from = ___, values_from = ___)   # which column gives the new column names? which one fills them?

dim(root_wide)

# the numbers only, as a matrix, with the sample names as row names
root_matrix <- root_wide %>%
  select(-sample, -genotype, -water, -heat) %>%
  as.matrix()
rownames(root_matrix) <- root_wide$sample

dim(root_matrix)

Checkpoint

dim() prints two numbers for each table: rows, then columns. root_matrix[1:3, 1:3] shows numbers only, with the sample names as row names.

C2. The PCA.

# scale. = TRUE or FALSE? (see C2 b)
root_pca <- prcomp(root_matrix, scale. = ___)

variance <- round(100 * summary(root_pca)$importance[2, ], 1)
variance

C3. Where are the samples? The PCA never saw water and heat: you only use them now, to colour the points.

root_scores <- as_tibble(root_pca$x[, 1:2]) %>%
  bind_cols(select(root_wide, sample, genotype, water, heat))

# colour by water regime, shape by heat regime; the variance of each axis in its title
ggplot(root_scores, aes(x = PC1, y = PC2, colour = ___, shape = ___)) +
  geom_point(size = 3) +
  labs(x = paste0("PC1 (", variance[1], "%)"),
       y = paste0("PC2 (", variance[___], "%)"))

C4. Which molecules drive the axes?

# sort the metabolites by their loading on the first axis
root_loadings <- as_tibble(root_pca$rotation[, 1:2], rownames = "metabolite") %>%
  arrange(___)

head(root_loadings, 5)   # most negative on PC1
tail(root_loadings, 5)   # most positive on PC1

C5. The heatmap. A small annotation table puts the conditions on top of the heatmap:

# a small table: one row per sample (sample names as row names), the two conditions
annotation <- root_wide %>%
  select(sample, ___, ___) %>%
  column_to_rownames("sample")

# scale each metabolite (each row) on its own; put the conditions on top
pheatmap(t(root_matrix), scale = ___, annotation_col = ___,
         show_colnames = FALSE, fontsize_row = 7)

5.4.6 D. Guided: with labels, then one test per molecule (16:35-17:10, pairs)

D1. PLS-DA with the real labels.

# the labels to learn: the water regime; 2 components
root_plsda <- plsda(root_matrix, root_wide$___, ncomp = ___)

plotIndiv(root_plsda, ellipse = TRUE, legend = TRUE,
          title = "Root, PLS-DA with the real labels")

D2. The trap. Now give the PLS-DA labels that mean nothing: the same “drought” and “watered” labels, mixed at random with sample(). Let it use 10 components, then ask how well it classifies the samples it learned from:

# the trap: the same labels, mixed at random, and 10 components
set.seed(1)
shuffled <- sample(root_wide$___)
root_plsda_shuffled <- plsda(root_matrix, ___, ncomp = ___)

# how often does it predict the shuffled label of the samples it learned from?
predicted <- predict(root_plsda_shuffled, root_matrix)$class$max.dist[, ___]
mean(predicted == ___)

set.seed(1) makes the random draw the same for everyone, so that you can compare your numbers.

D3. Cross-validation. The same two models, now judged on samples they did not learn from:

# cross-validation: read the error on the line of the last component of each model
set.seed(1)
cv_shuffled <- perf(root_plsda_shuffled, validation = "Mfold", folds = 4, nrepeat = 10)
1 - cv_shuffled$error.rate$BER[___, "max.dist"]

set.seed(1)
cv_real <- perf(root_plsda, validation = "Mfold", folds = 4, nrepeat = 10)
1 - cv_real$error.rate$BER[___, "max.dist"]

Checkpoint

D2 prints one number between 0 and 1 (the share of samples classified correctly), D3 prints two. subscript out of bounds means you asked for a component the model does not have: look at ncomp. The cross-validation numbers change a little from one run to the next: that is normal.

D4. One test per molecule. The t-test of Session 4, run on each metabolite, then the correction:

# one t-test per metabolite: drought vs watered
root_tests <- metabolites %>%
  filter(organ == "root") %>%
  group_by(___, panel) %>%
  summarise(log2_fold_change = log2(mean(concentration[water == ___]) /
                                    mean(concentration[water == ___])),
            p_value = t.test(log10(concentration) ~ ___, var.equal = TRUE)$p.value,
            .groups = "drop") %>%
  mutate(p_adjusted = p.adjust(p_value, method = ___))   # Benjamini-Hochberg

# how many are significant after the correction?
sum(root_tests$p_adjusted < ___)
root_tests %>% arrange(p_adjusted)

log2_fold_change: 1 means twice as high under drought, −1 half as high, 0 no change.

D5. The volcano plot.

# x: the effect, y: the evidence; a dashed line at the 5% level
ggplot(root_tests, aes(x = ___, y = -log10(___), colour = panel)) +
  geom_point(size = 2) +
  geom_hline(yintercept = -log10(___), linetype = "dashed") +
  geom_vline(xintercept = 0) +
  geom_text_repel(data = slice_min(root_tests, p_adjusted, n = 10),
                  aes(label = metabolite), size = 3, show.legend = FALSE)

5.4.7 E. On your own: the leaf

Same reasoning, on the leaf: the wide table and the matrix (E1, adapt C1), the PCA and its plot, saved and inserted into your report (E2), the PLS-DA with its cross-validated accuracy (E3, adapt D1 and D3), and one test per molecule (E4, adapt D4). Then compare the leaf with the root (E5). The workbook gives you the skeleton of E1 to E4, with more ___ than in C and D: fill them from your own code.

ggsave("results/leaf_pca.png", width = 6, height = 4)

Then insert results/leaf_pca.png into your report, as in Session 4, with your own caption and the label #fig-leaf, and cite it as @fig-leaf in E5.

5.4.8 F. On your own: the exudates, targeted then untargeted

The same plants seen by two methods. The workbook gives you the skeleton of the code: by now, most of the ___ are the choices you made in C and D, so reuse your own code:

F1, the targeted list (15 molecules): the PCA of C2-C3, on the exudates. Which stress dominates now?
F2, the untargeted table: one row per plant, one column per feature (F0001 to F2500), read directly in the wide format that metabolomics software produce. Some intensities are exactly 0, hence log10(x + 1). Does the PCA tell the same story as the targeted list?
F3, the PLS-DA trap again, with 2,500 variables and only 2 components.
F4, one test per feature: 2,500 tests, so the correction matters even more; then the names proposed for the significant features, in soybean_exudates_annotation.csv, and how much you can trust them; and the volcano plot of the 2,500 features, as in D5.
F5, putting the three compartments and the two methods together, and F6, one limitation.

5.4.9 Wrap-up (17:10-17:15)

Click Render. Check that your PDF shows all the plots of blocks B to D, warm-up included. Blocks E and F follow the same logic as C and D. Submit WORKBOOK.pdf on Moodle by Sunday evening (25 October). If you see errors in your PDF, or if it does not render, submit in addition your WORKBOOK.qmd and your results/ folder (your saved plots).

How this is graded (10 points)

Criterion	Points	What it looks for
Reproducible workflow	2	The workbook renders (or the code is complete and in order), inside a clean project.
Working analysis	3	All the code of blocks A to F (warm-up included) ran, and all your answers use the right numbers.
Required result	2	The leaf PCA (E2): coloured by water, shapes by heat, the variance of each axis in the axis titles, saved and inserted with a caption.
Biological interpretation	2	Your answers on the biology, block by block, and especially E5 and F5, with your own numbers.
Reflection	1	Your answers on the methods (the trap, the corrections, the putative names) and your limitation (F6).

Every question counts, except the exit questions.

5.5 Checkpoint

You are done when all of the following are true:

you can explain why concentrations are log-transformed and scaled before a PCA;
you turned a long table into a matrix with pivot_wider();
you can read a PCA: the share of variance, the scores and the loadings;
you can read a heatmap with clustering, and say what its trees are computed from;
you can explain why a PLS-DA that separates groups proves nothing until it is cross-validated, with the numbers of D2 and D3;
you can read a volcano plot, and say why the p-values were adjusted;
you can compare the leaf with the root, and say what is different;
you can say what targeted and untargeted metabolomics each see, and why a putative name from an untargeted search is not an identification;
you can say what the heatwave does to the exudates, and why it matters for the soil bacteria of Session 4.

5.6 Optional challenge

For your own interest, not graded.

A. The stem. Redo the PCA and the volcano plot for the stem. Is it closer to the root or to the leaf?

B. A heatmap of the leaf. Adapt C5 to the leaf. Does the first split of the top tree separate the water regimes, or the climates?

C. The heatwave, alone. Run the PLS-DA of the root on heat instead of water, with its cross-validation. Is the heatwave as easy to predict as drought? Look back at C3 to explain why.

D. Drought in the exudates. Redo F4 on water instead of heat. How many features change under drought, compared with the heatwave? Is it the same ranking of the two stresses as in the targeted list (F1)?

5.7 Going further

Principal component analysis and Clustering, Points of Significance (Nature Methods): two pages each.
mixOmics: tutorials on PCA, PLS-DA and their variants.
Modern Statistics for Modern Biology, chapters on multiple testing and multivariate analysis.

5.8 Before the next session

Session 6 scales up again: from 34 molecules to thousands of genetic markers, in a GWAS (genome-wide association study) on Arabidopsis. You will run thousands of the linear models of Session 4, correct for multiple testing as in today’s volcano plot, and use a PCA once more, this time to correct for the genetic structure of the population.

--- title: "Session 5: Many variables at once: PCA, heatmaps and PLS-DA" engine: knitr --- ## Overview [Download the lecture slides (PDF)](slides/05_lecture_slides.pdf) ::: callout-tip ## Quick poll (Poll Everywhere) Respond at [pe.app/cmaslard4](https://pe.app/cmaslard4) ::: ::: callout-note ## Biological question Drought and the heatwave reduced the plants' biomass (Session 4). Do they also change the plants' chemistry, the same way in the root, the leaf, and in what the roots release into the soil? ::: In Session 4, you analysed **one** variable at a time: the biomass, then the Shannon diversity. Today, each sample comes with **dozens, then thousands of variables at once**: the concentrations of 15 to 37 small molecules (metabolites) measured in the root, the stem, the leaf and the root exudates of the same soybean plants, then 2,500 untargeted signals from the exudates. One boxplot per molecule would give you dozens of plots and no overview. You will learn the classic tools to see the whole table at once: - **PCA** and **heatmaps with clustering**, which look for structure **without** being told which sample is which (unsupervised); - **PLS-DA**, which **is** told the groups and looks for what separates them (supervised), and the trap that comes with it; - the **volcano plot**: one test per molecule, as in Session 4, and why you must correct the p-values when you run many tests. ::: callout-important ## How you work and submit today As in Session 4, everything happens in `WORKBOOK.qmd` and you click **Render**. New today: Render makes a **PDF**, `WORKBOOK.pdf`, and that is what you submit. Blocks A to D are guided; blocks E (the leaf) and F (the exudates) you do on your own. Every block counts in the grade. Submit on Moodle **by Sunday evening (25 October)**. ::: ### Learning outcomes By the end of this session, you should be able to: - explain why concentrations are log-transformed and scaled before a multivariate analysis; - turn a long table into a sample × variable matrix with `pivot_wider()`; - run a PCA, read the share of variance of each axis, the scores (samples) and the loadings (variables); - draw a heatmap with clustering, and read its trees; - run a PLS-DA, and explain why a supervised method must always be cross-validated; - run one test per variable, correct the p-values for multiple testing, and draw a volcano plot; - say which method answers which question. ### Session plan | Time | Block | What happens | |---:|---|---| | 14:15-14:20 | Retrieval | Session 4 quiz. | | 14:20-14:35 | Theory 1 | Metabolomics: targeted and untargeted, and the data you get. | | 14:35-14:50 | Theory 2 | From one variable to many: PCA, clustering, PLS-DA, volcano. | | 14:50-15:00 | Practical briefing | What today's workbook will ask of you. | | 15:00-15:15 | Practical A | Install mixOmics (once); same project habits as Session 4, new data. | | 15:15-15:35 | Practical B | One molecule at a time: ABA, proline, and why we take the log. | | 15:35-15:45 | Break | | | 15:45-16:00 | Warm-up | Guided: what a PCA does, on decathlon athletes. | | 16:00-16:35 | Practical C | Guided: the whole root metabolome, without labels (PCA, heatmap). | | 16:35-17:10 | Practical D | Guided: PLS-DA and its trap, then one test per molecule. | | 17:10-17:15 | Wrap-up | Render, and how blocks E (the leaf) and F (the exudates) work. | ## Theory (≈45 min) ### 1. From one variable to many Each sample now has one value per metabolite. For the analyses of today, the table must have **one row per sample and one column per metabolite**: a **matrix**. Your data file is long (one row per sample *and* metabolite, as in Session 4), so the first step is `pivot_wider()`, the reverse of the `pivot_longer()` you used in Session 4. Two more steps before any multivariate method: - **log-transform**: concentrations go from 0.02 to about 190,000. On the raw scale, a few very abundant molecules would dominate everything. `log10()` turns ratios into differences: a molecule that doubles moves by the same amount, whether it is rare or abundant. - **scale**: each metabolite is centred and divided by its standard deviation, so that each one has the same weight in the analysis (`scale. = TRUE` in `prcomp()`). ### 2. Unsupervised: PCA and clustering A **PCA** (principal component analysis) summarises the 34 columns into a few new axes, the **principal components**. PC1 is the direction along which the samples differ the most, PC2 the next one, at right angles, and so on. Three things to read: - the **share of variance** of each axis: how much of the differences between samples it captures; - the **scores**: where each sample falls on the axes (the usual PCA plot); - the **loadings**: how much each metabolite contributes to each axis, and in which direction. A **heatmap with clustering** shows the whole matrix as colours, and reorders rows and columns so that similar ones sit next to each other. The trees (dendrograms) show the groups found. In the practical, you first run a PCA on data you know without any biology: 41 decathlon performances, 10 events each. You will see the axes rediscover the best athletes, and why scaling matters when the variables are in different units (seconds and metres). Both are **unsupervised**: they never see your labels (`water`, `heat`). If the samples still group by condition, the condition is a major source of variation in the data. ### 3. Supervised: PLS-DA, and why it must be validated A **PLS-DA** (partial least squares discriminant analysis) is told the groups, and looks for the combination of metabolites that **best separates** them. It is useful to find which molecules distinguish two conditions, but it is dangerous: with many variables and few samples, it can separate **any** groups, even groups drawn at random. A separation on the data the model learned from proves nothing. The answer is **cross-validation**: hide part of the samples, train the model on the rest, predict the hidden samples, repeat. Only the accuracy on samples the model did not see tells you whether the separation is real. You will see this trap with your own eyes in block D. ### 4. Many tests at once: the volcano plot You can also test each metabolite on its own: one t-test per molecule, as in Session 4. With 34 tests at the 5% level, you expect about 1 or 2 "significant" results **by chance alone**. Correcting the p-values (here with the Benjamini-Hochberg method, which controls the **false discovery rate**) keeps that under control. A **volcano plot** shows all the tests at once: the effect on the x axis (log2 fold change: positive means higher under drought), the evidence on the y axis (−log10 of the adjusted p-value: the higher, the stronger). **Which method, when?** | Your question | Method | In R | Supervised? | |---|---|---|---| | What are the main patterns in the samples? Do they group by condition? | PCA | `prcomp()` | no | | Which samples and which variables look alike? | heatmap with clustering | `pheatmap()` | no | | Which combination of variables separates known groups? | PLS-DA, **always with cross-validation** | `plsda()`, `perf()` | yes | | Which variables change between two groups, one by one? | one test per variable + FDR, volcano plot | `t.test()`, `p.adjust()` | (it is a test) | ### 5. Worked example: the metabolome of the soybean plants Today's data come from the plants you met in Session 4, measured with two kinds of metabolomics: | | Targeted | Untargeted | |---|---|---| | What is measured | a chosen list of known molecules | every signal the machine detects | | Each variable is | a named molecule, quantified | a "feature": a peak with a retention time and a mass (m/z) | | How many variables today | 15 to 37 per organ | 2,500 in the exudates | | Identity | known for sure | unknown for many; a **putative** name from a database search, with a confidence score | | Good for | testing precise hypotheses, comparing known pathways | discovering what you did not think of looking for | - **Root, stem, leaf** (targeted): the `omics_and_soil` plants, frozen at harvest and ground; 4 samples per cultivar and condition, each a pool of two plants. - **Exudates** (targeted **and** untargeted): what the roots of the `exudates` plants released into a solution, collected at harvest; one sample per plant. The exudates are the molecules the root gives to the soil: the food and the signals of the bacteria you studied in Session 4. ### Practical briefing Blocks A-D in pairs, blocks E and F on your own. Block B looks at single molecules you can understand without any statistics. The warm-up runs a first PCA on data you can check with common sense: decathlon athletes and their 10 events. Blocks C and D analyse the **root**, guided: first without labels (PCA, heatmap), then with labels (PLS-DA, volcano). In block E, you analyse the **leaf**: same reasoning, a different answer. In block F, the **exudates**, first with the targeted list, then with the untargeted scan: same plants, two very different tables. Your submission is the whole rendered workbook: every block counts. ## Your R toolbox New tools today, on top of the Session 4 toolbox (still valid: see the [Session 4 R toolbox](04_biological_tables_and_visualisation.qmd#your-r-toolbox)). ::: panel-tabset ## tidyr: long → wide `pivot_wider()` is the reverse of `pivot_longer()`: it spreads one column into many. ```r root_wide <- root_long %>% pivot_wider(names_from = metabolite, values_from = log_conc) ``` - `names_from`: the column whose values become the **new column names** (here: the metabolite names); - `values_from`: the column that fills the new columns (here: the log concentrations). Then `select(-sample, ...)` and `as.matrix()` keep the numbers only, as a matrix, and `rownames()` puts the sample names back on the rows. **Common mistake.** Keeping a column that differs inside a sample (such as the raw `concentration` next to `log_conc`): `pivot_wider()` then cannot put one sample on one row. Select only the columns you need before pivoting. ## PCA (base R) `prcomp()` is part of R itself: no package needed. | Code | What it gives | |---|---| | `prcomp(matrix, scale. = TRUE)` | the PCA, each variable scaled | | `summary(pca)$importance[2, ]` | the share of variance of each axis | | `pca$x` | the **scores**: one row per sample, one column per axis | | `pca$rotation` | the **loadings**: one row per variable, one column per axis | **Reading the loadings.** A metabolite with a large negative loading on PC1 is high in the samples on the left of the PCA plot; a large positive loading, high on the right. **Common mistake.** Forgetting `scale. = TRUE`: the variables with the largest numbers then decide everything. ## pheatmap ```r pheatmap(t(root_matrix), scale = "row", annotation_col = annotation, show_colnames = FALSE, fontsize_row = 7) ``` - `t()` turns the matrix around: metabolites in rows, samples in columns; - `scale = "row"`: each metabolite is shown as "above or below its own average" (z-score), otherwise the most abundant ones would fill the colour scale; - `annotation_col`: a small table, one row per sample (with the sample names as row names), whose columns become coloured bars on top; - the trees on the left and on top are the **clustering**: they are computed from the data, not from the annotation. ## mixOmics The reference package for PLS-DA (and much more) in biology. | Code | What it does | |---|---| | `plsda(matrix, groups, ncomp = 2)` | fits a PLS-DA with 2 components | | `plotIndiv(model, ellipse = TRUE, legend = TRUE)` | the sample plot, with one ellipse per group | | `predict(model, matrix)$class$max.dist` | the predicted group of each sample | | `perf(model, validation = "Mfold", folds = 4, nrepeat = 10)` | cross-validation: 4 folds, repeated 10 times | | `perf(...)$error.rate$BER` | the balanced error rate, per component (accuracy = 1 − error) | **Install once, load every time.** `BiocManager::install("mixOmics", update = FALSE)` runs once, in the Console (block A0). `library(mixOmics)` goes at the top of every report. Never put an install line in your workbook: every Render would install it again. **Common mistake: load it first.** mixOmics loads another package (MASS) that also has a function called `select()`. Whatever is loaded **last** wins. So load `mixOmics` **before** the tidyverse; otherwise `select()` stops working, with a confusing error ("unused argument"). **Documentation.** [mixomics.org](https://mixomics.org/) ## ggrepel `geom_text_repel()` writes text labels on a plot and moves them so that they do not overlap, unlike `geom_text()`. ```r geom_text_repel(data = slice_min(root_tests, p_adjusted, n = 10), aes(label = metabolite), size = 3) ``` `data = ...` labels only some of the points: here, the 10 with the smallest adjusted p-values (`slice_min()`, from dplyr). ::: ## Practical (≈1 h 50) ::: callout-tip ## Tested on LEGcompute `pheatmap` and `ggrepel` are already installed in RStudio Server's R. `mixOmics` is not: each of you installs it once, at the start of block A, in about 1 to 2 minutes (tested: 83 seconds, mixOmics 6.36.0). ::: ### A. Install mixOmics, same habits, new data (15:00-15:15) **A0. Install mixOmics, once.** In the **Console** (not in the workbook), run these three lines: ```r dir.create(Sys.getenv("R_LIBS_USER"), recursive = TRUE, showWarnings = FALSE) .libPaths(c(Sys.getenv("R_LIBS_USER"), .libPaths())) BiocManager::install("mixOmics", update = FALSE, ask = FALSE) ``` - line 1 creates your **personal R library**, a folder in your home where you can install packages (you cannot write in the shared one); - line 2 tells R to use it, from now on in this session; - line 3 installs mixOmics from Bioconductor (a repository of biology packages), with the few packages it still needs. `update = FALSE` is important: otherwise R offers to update dozens of packages you already have. It takes 1 to 2 minutes, and prints a lot of compilation messages: that is normal. While it runs, go on with the Terminal commands below. ::: callout-note ## Checkpoint The Console ends with `* DONE (mixOmics)`. Then `library(mixOmics)` in the Console loads it without error (a few welcome lines are printed). You will never need to install it again: from now on, R finds it in your personal library. ::: **A1. The project.** In the **Terminal** tab: ```bash mkdir -p ~/bioinfo_course/session_05/data ~/bioinfo_course/session_05/img ~/bioinfo_course/session_05/results cp /home/cmaslard/bioinfo_tools_course/session_05/soybean_metabolites.csv ~/bioinfo_course/session_05/data/ cp /home/cmaslard/bioinfo_tools_course/session_05/soybean_exudates_untargeted.csv ~/bioinfo_course/session_05/data/ cp /home/cmaslard/bioinfo_tools_course/session_05/soybean_exudates_annotation.csv ~/bioinfo_course/session_05/data/ cp /home/cmaslard/bioinfo_tools_course/session_05/decathlon.csv ~/bioinfo_course/session_05/data/ cp /home/cmaslard/bioinfo_tools_course/session_05/experimental_design.png ~/bioinfo_course/session_05/img/ cp /home/cmaslard/bioinfo_tools_course/session_05/WORKBOOK.qmd ~/bioinfo_course/session_05/ ls -R ~/bioinfo_course/session_05 ``` Then **File → New Project → Existing Directory**, pick `~/bioinfo_course/session_05`, open `WORKBOOK.qmd`, write your name at the top and click **Render** once. ::: callout-note ## Checkpoint `WORKBOOK.pdf` appears (RStudio opens it), with the design figure as Figure 1, and the first chunk lists the columns `sample`, `organ`, `genotype`, `water`, `heat`, `panel`, `metabolite`, `kegg_id`, `concentration`. A few chunks further down show `Error: unexpected input`: those are the `___` you will fill. ::: `panel` is the analytical panel that measured each compound (amino acids, organic acids, flavonoids, phytohormones, cytokinins), as in the lab's files: a practical grouping, not always the strict chemical family. `concentration` was cleaned before class: values below the detection limit were replaced by half of the smallest value measured for that molecule, a common rule in metabolomics. ### B. One molecule at a time (15:10-15:35) **B1. What is in the table?** ```r # how many samples per organ? (count by the column that holds the organ) metabolites %>% distinct(organ, sample) %>% count(___) metabolites %>% distinct(organ, metabolite) %>% count(organ) ``` ::: callout-note ## Checkpoint Two small tables, each with one row per organ: `exudate`, `leaf`, `root`, `stem`. If you see one row per sample instead, check what you put in `count()`. ::: **B2. The drought hormone.** Abscisic acid (ABA) makes the leaves close their stomata, the small pores through which they lose water. Look at it in the three organs (it is not on the exudate list): ```r aba <- metabolites %>% filter(metabolite == "Abscisic Acid") # one panel per organ ggplot(aba, aes(x = water, y = concentration, colour = heat)) + geom_boxplot(outlier.shape = NA) + geom_point(position = position_jitterdodge(jitter.width = 0.1)) + facet_wrap(~ ___, scales = "free_y") # the mean per organ and per water regime aba %>% group_by(organ, ___) %>% summarise(mean_aba = mean(concentration)) ``` `scales = "free_y"` gives each organ its own y axis: the levels are very different between organs. **B3. Proline**, a small amino acid that protects cells when water is short. Same plot, with one `___` to fill. **B4. Why the log?** ```r range(metabolites$___) ggplot(metabolites, aes(x = concentration)) + geom_histogram(bins = 50) # the same histogram, on the log10 scale ggplot(metabolites, aes(x = log10(___))) + geom_histogram(bins = 50) ``` ::: callout-note ## Checkpoint `range()` prints two numbers. The first histogram is one tall bar near zero; the second one is spread out and readable. ::: ### 15:35-15:45: Short break ### Warm-up: what a PCA does, on athletes (15:45-16:00, pairs) Before 34 metabolites, a PCA on data you can check with common sense: 41 decathlon performances (Decastar meeting and Athens Olympic Games, 2004; data of the FactoMineR package). Each athlete did 10 events: four runs (times in seconds, small is good) and six jumps and throws (metres, large is good). **W1. The table.** ```r decathlon <- read_csv("data/decathlon.csv") glimpse(decathlon) # how many athletes per competition? decathlon %>% count(___) ``` **W2. Without and with scaling.** The same PCA twice: first with each event in its own unit, then with every event scaled. ```r # the 10 events only, as a matrix, with the athletes' names as row names events <- decathlon %>% select(run_100m:run_1500m) %>% as.matrix() rownames(events) <- decathlon$___ # how much does each event vary between athletes? (standard deviation, in its own unit) decathlon %>% summarise(across(run_100m:run_1500m, sd)) # a PCA WITHOUT scaling: each event keeps its own unit (seconds, metres) pca_raw <- prcomp(events, scale. = FALSE) round(pca_raw$rotation[, 1], 2) # the loadings on PC1 # the same PCA WITH scaling: every event gets the same weight decathlon_pca <- prcomp(events, scale. = ___) decathlon_variance <- round(100 * summary(decathlon_pca)$importance[2, ], 1) decathlon_variance ``` `across(run_100m:run_1500m, sd)` applies `sd()` to every column from `run_100m` to `run_1500m`. **W3. Arrows and points.** The **loadings** become arrows (one per event), the **scores** points (one per athlete): ```r # the events: one arrow per event, from 0 to its loadings on PC1 and PC2 decathlon_loadings <- as_tibble(decathlon_pca$rotation[, 1:2], rownames = "event") ggplot(decathlon_loadings) + geom_segment(aes(x = 0, y = 0, xend = ___, yend = ___), arrow = arrow(length = unit(0.2, "cm"))) + geom_text_repel(aes(x = PC1, y = PC2, label = ___)) + coord_equal() + labs(x = paste0("PC1 (", decathlon_variance[1], "%)"), y = paste0("PC2 (", decathlon_variance[2], "%)")) # the athletes: one point per athlete, coloured by the final score decathlon_scores <- as_tibble(decathlon_pca$x[, 1:2]) %>% bind_cols(select(decathlon, athlete, competition, points)) ggplot(decathlon_scores, aes(x = PC1, y = PC2, colour = ___)) + geom_point(size = 3) + geom_text_repel(aes(label = athlete), size = 3) + labs(x = paste0("PC1 (", decathlon_variance[1], "%)"), y = paste0("PC2 (", decathlon_variance[___], "%)")) # does PC1 follow the final score? cor(decathlon_scores$PC1, decathlon_scores$___) ``` ::: callout-note ## Checkpoint Two plots: one with 10 arrows (one per event), one with 41 points. Some names appear twice, because 9 athletes competed in both meetings. On another computer, the whole plot may be flipped left-right or up-down: the direction of an axis is arbitrary, only the angles between arrows matter. ::: ### C. Guided: the whole root metabolome, without labels (16:00-16:35, pairs) **C1. From long to wide.** One row per sample, one column per metabolite, log10 of the concentration: ```r # one row per sample, one column per metabolite, log10 of the concentration root_wide <- metabolites %>% filter(organ == ___) %>% mutate(log_conc = ___(concentration)) %>% select(sample, genotype, water, heat, metabolite, log_conc) %>% pivot_wider(names_from = ___, values_from = ___) # which column gives the new column names? which one fills them? dim(root_wide) # the numbers only, as a matrix, with the sample names as row names root_matrix <- root_wide %>% select(-sample, -genotype, -water, -heat) %>% as.matrix() rownames(root_matrix) <- root_wide$sample dim(root_matrix) ``` ::: callout-note ## Checkpoint `dim()` prints two numbers for each table: rows, then columns. `root_matrix[1:3, 1:3]` shows numbers only, with the sample names as row names. ::: **C2. The PCA.** ```r # scale. = TRUE or FALSE? (see C2 b) root_pca <- prcomp(root_matrix, scale. = ___) variance <- round(100 * summary(root_pca)$importance[2, ], 1) variance ``` **C3. Where are the samples?** The PCA never saw `water` and `heat`: you only use them now, to colour the points. ```r root_scores <- as_tibble(root_pca$x[, 1:2]) %>% bind_cols(select(root_wide, sample, genotype, water, heat)) # colour by water regime, shape by heat regime; the variance of each axis in its title ggplot(root_scores, aes(x = PC1, y = PC2, colour = ___, shape = ___)) + geom_point(size = 3) + labs(x = paste0("PC1 (", variance[1], "%)"), y = paste0("PC2 (", variance[___], "%)")) ``` **C4. Which molecules drive the axes?** ```r # sort the metabolites by their loading on the first axis root_loadings <- as_tibble(root_pca$rotation[, 1:2], rownames = "metabolite") %>% arrange(___) head(root_loadings, 5) # most negative on PC1 tail(root_loadings, 5) # most positive on PC1 ``` **C5. The heatmap.** A small annotation table puts the conditions on top of the heatmap: ```r # a small table: one row per sample (sample names as row names), the two conditions annotation <- root_wide %>% select(sample, ___, ___) %>% column_to_rownames("sample") # scale each metabolite (each row) on its own; put the conditions on top pheatmap(t(root_matrix), scale = ___, annotation_col = ___, show_colnames = FALSE, fontsize_row = 7) ``` ### D. Guided: with labels, then one test per molecule (16:35-17:10, pairs) **D1. PLS-DA with the real labels.** ```r # the labels to learn: the water regime; 2 components root_plsda <- plsda(root_matrix, root_wide$___, ncomp = ___) plotIndiv(root_plsda, ellipse = TRUE, legend = TRUE, title = "Root, PLS-DA with the real labels") ``` **D2. The trap.** Now give the PLS-DA labels that mean **nothing**: the same "drought" and "watered" labels, mixed at random with `sample()`. Let it use 10 components, then ask how well it classifies the samples it learned from: ```r # the trap: the same labels, mixed at random, and 10 components set.seed(1) shuffled <- sample(root_wide$___) root_plsda_shuffled <- plsda(root_matrix, ___, ncomp = ___) # how often does it predict the shuffled label of the samples it learned from? predicted <- predict(root_plsda_shuffled, root_matrix)$class$max.dist[, ___] mean(predicted == ___) ``` `set.seed(1)` makes the random draw the same for everyone, so that you can compare your numbers. **D3. Cross-validation.** The same two models, now judged on samples they did not learn from: ```r # cross-validation: read the error on the line of the last component of each model set.seed(1) cv_shuffled <- perf(root_plsda_shuffled, validation = "Mfold", folds = 4, nrepeat = 10) 1 - cv_shuffled$error.rate$BER[___, "max.dist"] set.seed(1) cv_real <- perf(root_plsda, validation = "Mfold", folds = 4, nrepeat = 10) 1 - cv_real$error.rate$BER[___, "max.dist"] ``` ::: callout-note ## Checkpoint D2 prints one number between 0 and 1 (the share of samples classified correctly), D3 prints two. `subscript out of bounds` means you asked for a component the model does not have: look at `ncomp`. The cross-validation numbers change a little from one run to the next: that is normal. ::: **D4. One test per molecule.** The t-test of Session 4, run on each metabolite, then the correction: ```r # one t-test per metabolite: drought vs watered root_tests <- metabolites %>% filter(organ == "root") %>% group_by(___, panel) %>% summarise(log2_fold_change = log2(mean(concentration[water == ___]) / mean(concentration[water == ___])), p_value = t.test(log10(concentration) ~ ___, var.equal = TRUE)$p.value, .groups = "drop") %>% mutate(p_adjusted = p.adjust(p_value, method = ___)) # Benjamini-Hochberg # how many are significant after the correction? sum(root_tests$p_adjusted < ___) root_tests %>% arrange(p_adjusted) ``` `log2_fold_change`: 1 means twice as high under drought, −1 half as high, 0 no change. **D5. The volcano plot.** ```r # x: the effect, y: the evidence; a dashed line at the 5% level ggplot(root_tests, aes(x = ___, y = -log10(___), colour = panel)) + geom_point(size = 2) + geom_hline(yintercept = -log10(___), linetype = "dashed") + geom_vline(xintercept = 0) + geom_text_repel(data = slice_min(root_tests, p_adjusted, n = 10), aes(label = metabolite), size = 3, show.legend = FALSE) ``` ### E. On your own: the leaf Same reasoning, on the **leaf**: the wide table and the matrix (E1, adapt C1), the PCA and its plot, saved and inserted into your report (E2), the PLS-DA with its cross-validated accuracy (E3, adapt D1 and D3), and one test per molecule (E4, adapt D4). Then compare the leaf with the root (E5). The workbook gives you the skeleton of E1 to E4, with more `___` than in C and D: fill them from your own code. ```r ggsave("results/leaf_pca.png", width = 6, height = 4) ``` Then insert `results/leaf_pca.png` into your report, as in Session 4, with your own caption and the label `#fig-leaf`, and cite it as `@fig-leaf` in E5. ### F. On your own: the exudates, targeted then untargeted The same plants seen by two methods. The workbook gives you the skeleton of the code: by now, most of the `___` are the choices you made in C and D, so reuse your own code: - **F1**, the targeted list (15 molecules): the PCA of C2-C3, on the exudates. Which stress dominates now? - **F2**, the untargeted table: one row per plant, one column per feature (`F0001` to `F2500`), read directly in the wide format that metabolomics software produce. Some intensities are exactly 0, hence `log10(x + 1)`. Does the PCA tell the same story as the targeted list? - **F3**, the PLS-DA trap again, with 2,500 variables and only 2 components. - **F4**, one test per feature: 2,500 tests, so the correction matters even more; then the names proposed for the significant features, in `soybean_exudates_annotation.csv`, and how much you can trust them; and the volcano plot of the 2,500 features, as in D5. - **F5**, putting the three compartments and the two methods together, and **F6**, one limitation. ### Wrap-up (17:10-17:15) Click **Render**. Check that your PDF shows all the plots of blocks B to D, warm-up included. Blocks E and F follow the same logic as C and D. Submit `WORKBOOK.pdf` on Moodle **by Sunday evening (25 October)**. If you see errors in your PDF, or if it does not render, submit in addition your `WORKBOOK.qmd` and your `results/` folder (your saved plots). ::: callout-note ## How this is graded (10 points) | Criterion | Points | What it looks for | |---|---:|---| | Reproducible workflow | 2 | The workbook renders (or the code is complete and in order), inside a clean project. | | Working analysis | 3 | All the code of blocks A to F (warm-up included) ran, and all your answers use the right numbers. | | Required result | 2 | The leaf PCA (E2): coloured by water, shapes by heat, the variance of each axis in the axis titles, saved and inserted with a caption. | | Biological interpretation | 2 | Your answers on the biology, block by block, and especially E5 and F5, with your own numbers. | | Reflection | 1 | Your answers on the methods (the trap, the corrections, the putative names) and your limitation (F6). | Every question counts, except the exit questions. ::: ## Checkpoint You are done when all of the following are true: - you can explain why concentrations are log-transformed and scaled before a PCA; - you turned a long table into a matrix with `pivot_wider()`; - you can read a PCA: the share of variance, the scores and the loadings; - you can read a heatmap with clustering, and say what its trees are computed from; - you can explain why a PLS-DA that separates groups proves nothing until it is cross-validated, with the numbers of D2 and D3; - you can read a volcano plot, and say why the p-values were adjusted; - you can compare the leaf with the root, and say what is different; - you can say what targeted and untargeted metabolomics each see, and why a putative name from an untargeted search is not an identification; - you can say what the heatwave does to the exudates, and why it matters for the soil bacteria of Session 4. ## Optional challenge For your own interest, not graded. **A. The stem.** Redo the PCA and the volcano plot for the stem. Is it closer to the root or to the leaf? **B. A heatmap of the leaf.** Adapt C5 to the leaf. Does the first split of the top tree separate the water regimes, or the climates? **C. The heatwave, alone.** Run the PLS-DA of the root on `heat` instead of `water`, with its cross-validation. Is the heatwave as easy to predict as drought? Look back at C3 to explain why. **D. Drought in the exudates.** Redo F4 on `water` instead of `heat`. How many features change under drought, compared with the heatwave? Is it the same ranking of the two stresses as in the targeted list (F1)? ## Going further - [Principal component analysis](https://pure.psu.edu/en/publications/points-of-significance-principal-component-analysis/) and [Clustering](https://pure.psu.edu/en/publications/points-of-significance-clustering/), Points of Significance (Nature Methods): two pages each. - [mixOmics](https://mixomics.org/): tutorials on PCA, PLS-DA and their variants. - [Modern Statistics for Modern Biology](https://www.huber.embl.de/msmb/), chapters on multiple testing and multivariate analysis. ## Before the next session Session 6 scales up again: from 34 molecules to thousands of genetic markers, in a **GWAS** (genome-wide association study) on *Arabidopsis*. You will run thousands of the linear models of Session 4, correct for multiple testing as in today's volcano plot, and use a PCA once more, this time to correct for the genetic structure of the population.