# Task overview > In science, we collect data and use this data to assess the evidence for or against a null hypothesis. In the context of omics studies, we do not test one hypothesis, but thousands of simultaneous hypotheses. For example, in differential gene expression analysis, we test each gene for differential expression between treated and controls. Each individual test comes with a risk of making a false conclusion based on the data, and these risks add up as more tests are performed. In this practical, we will demonstrate the problems associated with **test multiplicity**, and show how we can correct for this. > In the first part, you will see a **demonstration with simulated data**. Because we generate the data ourselves, we know which variables are true positives (i.e. for which the exposure has a true effect between treated and control). This is useful, as it allows us to examine the rates at which we make errors. In this part, the code will mainly be given, and you will asked to interpret some of the results. > In the second part, you will perform a basic differential gene expression analysis on real data, following the example from a previous practical called `PHDS_omics_GMA_2026_solutions.Rmd`. We will use this example to show the difference between the results when using different choices of multiplicity correction. # Setup ## 0a) Set up your project ### Why do we need a project? > When you analyse data in R, you need to tell R *where* files are. For example, you might load data into your environment by defining a path `"C:/Users/me/Documents/data/file.xlsx"`. However, this is problematic, since it is specific to your own computer setup. The code will then break as soon as you move the files, share them with a colleague, or re-run the analysis on a different machine. > To avoid this, we will introduce some basic organisational skills for an R project. Firstly, we will contain the work in a single **R project folder** (the root), and refer to every file *relative to that folder*. This makes specifying paths for loading data and saving outputs more robust. > Using a project folder also helps to organise files in a clean way. > > - **raw data** : the original data files, which you should never modify; > - **analysis** : the code (such as this `.Rmd` file) that reads the raw data and produces results; > - **output** : everything the code creates (figures, tables, processed data), which should be able to be perfectly regenerated by re-running the code. ### Create the project folder > If you already created a project for a previous practical, you can reuse it : open it and go to "Add the files" below. Otherwise : > 1. In RStudio, click `File > New Project... > New Directory > New Project`. > 2. Give the project a name, for example `omics_practicals`, and choose where to save it on your computer. Click `Create Project`. > 3. RStudio restarts inside your new project. You can check this in the top-right corner of RStudio, which now shows the project name, and in the `Files` pane, which shows a single file called e.g. `omics_practicals.Rproj`. > Next, create three sub-folders inside the project. In the `Files` pane of RStudio, click `New Folder` and create `data-raw`, `analysis`, and `output`. Alternatively, run the following code once in the R console. ```{r, eval = F} dir.create("data-raw") dir.create("analysis") dir.create("output") ``` > Your project should now look like this: ``` omics_practicals/ ├── omics_practicals.Rproj <- marks the "root" of the project ├── data-raw/ <- the raw data files ├── analysis/ <- the practicals (.Rmd files) └── output/ <- figures and results that you save ``` ### Add the files > 1. Copy the data file I sent you, `TB_batch_corrected.RData`, into the `data-raw` folder (we will use it in the second part). > 2. Copy this practical (the `.Rmd` file) into the `analysis` folder. > 3. In the `Files` pane of RStudio, open the `analysis` folder and click on the `.Rmd` file to open it. > **Important:** always start your work by opening the project (double-click the `.Rproj` file, or use `File > Open Project`), and open the `.Rmd` file from *within* RStudio's `Files` pane. If you open the `.Rmd` file directly from your Downloads folder, it will not be part of the project and the paths below will not work. ### How do we refer to files? > By default, R Markdown runs code from the folder where the `.Rmd` file is saved (here, `analysis/`), which would make paths confusing. Instead, we will use the `here` package. The function `here()` always returns the path to the root of your project (the folder containing the `.Rproj` file), wherever the code is run from. We then build paths from this root. For example, `here("output", "my_plot.pdf")` points to a file in the `output` folder. The same code works on any computer and operating system, and you never have to type a full path. ## 0b) Install and load the required packages > If you ran the installation script I sent you, all the packages are already installed. Otherwise, the following chunk checks which of the packages needed for this practical you are missing and installs them : `dplyr`, `tidyr`, `ggplot2` and `here` (from CRAN), and `edgeR` and `limma` (from Bioconductor). ```{r, echo = F, results = 'hide', message = F} # This sets the options for the outputs when knitting the rmd file knitr::opts_chunk$set(warning = F, message = T, fig.width = 8, fig.height = 5) ``` ```{r, eval = F, results = 'hide', message = FALSE, warning = FALSE} cran_packages <- c("dplyr", "tidyr", "ggplot2", "here") # Install missing CRAN packages install.packages(setdiff(cran_packages, installed.packages()[, "Package"])) # ---- BioConductor Packages ---- bioc_packages <- c("edgeR", "limma") # Load BiocManager and install missing Bioconductor packages if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install(setdiff(bioc_packages, installed.packages()[, "Package"])) ``` > Now load the packages. ```{r, echo = F, results = 'hide', message = F} library(dplyr) library(tidyr) library(ggplot2) library(here) library(edgeR) library(limma) # Default look of all the ggplot2 figures theme_set(theme_minimal(base_size = 14)) # Set a seed for reproducibility set.seed(123) ``` ## 0c) Check your project setup > Run the following chunk to check that everything is in the right place. The first line should print the path to your project folder (the one containing the `.Rproj` file), and the data file should be found. ```{r} # Path to the root of the project here() # Check that the data file exists file.exists(here("data-raw", "TB_batch_corrected.RData")) ``` # 1) Demonstration on simulated data ## 1a) One gene with no true signal > Our setting is the following. We measure the expression of a gene in 10 controls and in 10 treated, and we want to know whether the gene is **differentially expressed** : is its average expression different between treated and controls? One basic way to answer this could be with a two-sample t-test, which gives a **p-value**. Conventionally, we call the gene "significant" if its p-value is below 0.05. > Since we generate the data, we choose whether there is a true difference between the generation of this gene between treated and controls. In this first example, we draw the expression of the controls from a normal distribution with mean 0 and standard deviation 1, and that of the treated from a normal distribution with mean `delta` and standard deviation 1. The number `delta` is the **true difference** in mean expression between treated and controls, in standard deviations. If `delta = 0`, the gene is truly not differentially expressed. Run the following chunk : it simulates one gene with no true difference. ```{r} n <- 10 # number of controls, and number of treated delta <- 0 # true difference between treated and controls : 0 means that the gene is not different control <- rnorm(n, mean = 0, sd = 1) # expression of the gene in the 10 controls treated <- rnorm(n, mean = delta, sd = 1) # expression of the gene in the 10 treated # Mean expression in each group, and p-value of the t-test c(control = mean(control), treated = mean(treated)) t.test(control, treated, var.equal = TRUE)$p.value ``` > **Try it** : run the chunk several times. The data are drawn at random, so the p-value will change in value every time, although the gene is never truly different. This is due to **sampling variability**. > To study what happens in general, we need to repeat the experiment many times. The following chunk defines a function, `simulate_pvalues()`, which repeats exactly the experiment above for many genes : for each element of `delta` (the true difference of one gene), it draws the expression of the controls and of the treated, runs the t-test, and returns the p-value. The result is a vector with one p-value per gene. ```{r} simulate_pvalues <- function(delta, n = 10) { sapply(delta, function(d) { control <- rnorm(n, mean = 0, sd = 1) treated <- rnorm(n, mean = d, sd = 1) t.test(control, treated, var.equal = TRUE)$p.value }) } ``` > Let's now simulate 100 genes, none of which is truly different (`delta = 0` for all genes), and look at their p-values. The dashed horizontal line shows what a perfectly flat distribution would look like, and the red line marks 0.05. **Try it** : now increase the number of genes (for example to 1,000, 10,000 or 100,000) and re-run the chunk. What happens to the shape of the histogram when the number of genes increases? What distribution does this resemble? Why do you think this is? ```{r} n_genes <- 100 # number of genes : change this value and re-run the chunk p_null <- simulate_pvalues(delta = rep(0, n_genes)) prop_signif_null <- mean(p_null < 0.05) plot_null <- ggplot(data.frame(p = p_null), aes(x = p)) + geom_histogram(aes(y = after_stat(density)), breaks = seq(0, 1, 0.05), colour = "white") + geom_hline(yintercept = 1, linetype = "dashed") + geom_vline(xintercept = 0.05, colour = "red", linetype = "dashed") + labs(x = "p-value", y = "Density", title = sprintf("p-values of %s genes with no true difference", format(n_genes, big.mark = ",")), subtitle = sprintf("%.1f%% of the genes have a p-value below 0.05", 100 * prop_signif_null)) ggsave(here("output", "pvalues_no_signal.pdf"), plot_null, width = 8, height = 5) plot_null ``` > **Answer:** With a small number of genes (100), the histogram is irregular : there are too few p-values to see a pattern, and the bars differ by chance. The more genes, the closer the histogram gets to a flat line (the dashed line) : it resembles a **uniform distribution** between 0 and 1. This is not a coincidence, it is how a p-value is defined. The p-value of a test is the probability, *if there is no true difference*, of observing a result at least as extreme as the one we observed. So, among many experiments with no true difference, 5% of them give a result among the 5% most extreme, which means that 5% have a p-value below 0.05, 10% have a p-value below 0.10, and so on : the proportion of p-values below any value x is exactly x. In other words, the cumulative distribution function of the p-values is F(x) = x, which is the cumulative distribution function of the uniform distribution on (0, 1). > What proportion of the genes with no true difference are declared significant at the 5% level? How does this proportion change when you increase the number of genes? > **Answer:** With the `r format(n_genes, big.mark = ",")` genes of the chunk above, `r sprintf("%.1f", 100 * prop_signif_null)`% of the genes are declared significant. As the number of genes increases, this proportion gets closer and closer to 5%. Since the p-values are uniform, the probability that a p-value falls below 0.05 is exactly 0.05. This probability of declaring significant a gene that is not different is the **false positive rate** (FPR) of the test : by choosing a threshold of 0.05, we choose to accept a 5% false positive rate. ## 1b) Genes with a true difference > So, we just learned about the false positive rate, i.e. the proportion of true negatives that we declare positive after a statistical test. For a **valid** statistical test, the false positive rate is expected to match the **claimed** rate (i.e. the significance level of the test) exactly. > Now, what happens when the gene *is* truly differentially expressed? The below chunk simulates 10,000 genes with a small difference in expression ($\delta = 0.5$). **Try it** : re-run the chunk with increasingly large values of $\delta$ (for example 1, 2 and 3). What happens to the histogram of p-values? How does this differ from what we previously saw when $\delta = 0$, and why? ```{r} n_genes <- 10000 # number of genes delta <- 0.5 # true difference of the genes : change this value and re-run the chunk p_signal <- simulate_pvalues(rep(delta, n_genes)) prop_signif_signal <- mean(p_signal < 0.05) plot_signal <- ggplot(data.frame(p = p_signal), aes(x = p)) + geom_histogram(breaks = seq(0, 1, 0.05), colour = "white") + geom_vline(xintercept = 0.05, colour = "red", linetype = "dashed") + labs(x = "p-value", y = "Number of genes", title = sprintf("p-values of 10,000 genes with a true difference of delta = %s", delta), subtitle = sprintf("%.1f%% of the genes have a p-value below 0.05", 100 * prop_signif_signal)) ggsave(here("output", "pvalues_signal.pdf"), plot_signal, width = 8, height = 5) plot_signal ``` > **Answer:** When `delta = 0.5`, the histogram is no longer flat : there are more small p-values than large ones. The larger `delta`, the more the p-values are concentrated near 0, and the more genes are below the red line (0.05). This differs from `delta = 0`, where all p-values were equally likely. When there is a true difference, the observed difference between the groups tends to be large compared to the variability within the groups, which is unlikely if there was no true difference, and hence gives a small p-value. The larger the true difference, the more often this happens. > Among the genes with a true difference, the proportion of genes declared significant (i.e. with a p-value below 0.05) is called the **power** of the test : the probability to detect a real difference. The following chunk computes it for a range of values of `delta`, with 10, 20 and 50 participants per group. What do you notice about the results of this experiment? ```{r} power_df <- expand.grid(delta = seq(0, 3, 0.25), n = c(10, 20, 50)) %>% rowwise() %>% mutate(proportion_significant = mean(simulate_pvalues(rep(delta, 10000), n = n) < 0.05)) %>% ungroup() %>% mutate(n = factor(n, labels = c("10 per group", "20 per group", "50 per group"))) plot_power <- ggplot(power_df, aes(x = delta, y = proportion_significant, colour = n)) + geom_hline(yintercept = 0.05, linetype = "dashed") + geom_line() + geom_point() + labs(x = "True difference (delta)", y = "Proportion of genes declared significant", colour = "Participants", title = "Proportion of genes declared significant, for different true differences") ggsave(here("output", "power_curve.pdf"), plot_power, width = 8, height = 5) plot_power ``` > **Answer:** There are two things to notice. (1) The power increases with the size of the true difference : from about 5% when `delta = 0` (where this proportion is the false positive rate) to `r sprintf("%.0f", 100 * power_df$proportion_significant[power_df$delta == 1 & power_df$n == "10 per group"])`% when `delta = 1`, and almost 100% when `delta = 3`, with 10 participants per group. The larger the difference, the easier it is to detect. (2) The power increases with the number of participants : for the same `delta`, the curve is higher with 50 participants per group than with 10 (for `delta = 1`, `r sprintf("%.0f", 100 * power_df$proportion_significant[power_df$delta == 1 & power_df$n == "50 per group"])`% against `r sprintf("%.0f", 100 * power_df$proportion_significant[power_df$delta == 1 & power_df$n == "10 per group"])`%), so a smaller difference is enough to be detected almost every time. This is because, with more participants, the mean of each group is estimated more precisely, so a given true difference stands out more clearly from the random variation, and the test is more likely to give a small p-value. ## 1c) A mix of genes with and without a true difference > In an omics study, we test thousands of genes, and most of them are not differentially expressed. We call the genes that are truly differentially expressed the **positives**, and the genes that are not the **negatives**. After testing, each gene is either declared significant or not. Comparing our decision with the truth, each gene falls into one of four categories : > > - **TP** (true positive) : a positive, declared significant; > - **FP** (false positive) : a negative, declared significant; > - **FN** (false negative) : a positive, not declared significant; > - **TN** (true negative) : a negative, not declared significant. > From these four numbers, we compute three rates : > > - the **false positive rate** (FPR) = FP / (FP + TN) : among the negatives, the proportion declared significant; > - the **power** = TP / (TP + FN) : among the positives, the proportion declared significant; > - the **false discovery proportion** (FDP) = FP / (FP + TP) : among the genes declared significant, the proportion which are negatives. > The following helper function counts TP, FP, FN and TN from our decisions (`significant`, `TRUE` if the gene is declared significant) and from the truth (`is_pos`, `TRUE` if the gene is a positive), and computes the three rates. ```{r} evaluate_decision <- function(significant, is_pos) { TP <- sum(significant & is_pos) FP <- sum(significant & !is_pos) FN <- sum(!significant & is_pos) TN <- sum(!significant & !is_pos) data.frame( TP = TP, FP = FP, FN = FN, TN = TN, power = TP / (TP + FN), FPR = FP / (FP + TN), FDP = ifelse(TP + FP == 0, 0, FP / (FP + TP)) # 0 if there is no discovery ) } ``` > Let's simulate a first small study : **100 positives** (with `delta = 3`, so the power is almost 100%, we expect to detect all 100 positives) and **100 negatives**. We test all the genes, and declare significant those with a p-value below 0.05. The chunk prints the contingency table of the truth against our decisions, and the numbers and rates computed by the helper function. > **Try it** : now increase the number of negatives (`n_neg`) to 1,000, and then to 10,000, keeping 100 positives, and re-run the chunk each time. What happens to the power, to the false positive rate and to the false discovery proportion? ```{r} n_pos <- 100 # number of positives : genes which are truly different n_neg <- 100 # number of negatives : genes which are not different delta_pos <- 3 # true difference of the positives # The truth for each gene, and the true difference of each gene is_pos <- rep(c(TRUE, FALSE), c(n_pos, n_neg)) delta <- ifelse(is_pos, delta_pos, 0) # Simulate and test all the genes, and declare significant those with a p-value below 0.05 p_mix <- simulate_pvalues(delta) significant <- p_mix < 0.05 # Contingency table : the truth against our decision table( Truth = factor(ifelse(is_pos, "Positive", "Negative"), levels = c("Positive", "Negative")), Decision = factor(ifelse(significant, "Significant", "Not significant"), levels = c("Significant", "Not significant")) ) # The numbers TP, FP, FN, TN and the three rates results_mix <- evaluate_decision(significant, is_pos) results_mix ``` > **Answer:** The power (almost 100%) and the false positive rate (about 5%) do not change. The false discovery proportion increases a lot (from about 30% with 1,000 negatives to more than 80% with 10,000 negatives). > How do you explain that the false discovery proportion increases, although the power and the false positive rate stay the same (up to sampling variability)? **Hint**: use the error definitions above to justify this. > **Answer:** The false positive rate is a proportion *of the negatives* : 5% of them are declared significant, whatever their number. So the number of false positives is about 5% of the number of negatives, and it grows with it : with 10,000 negatives, about 500 of them are declared significant. The number of true positives, on the other hand, cannot exceed the 100 positives. The false discovery proportion, FP / (FP + TP), is therefore more and more dominated by the false positives. A 5% false positive rate guarantees that 5% of the negatives are wrongly declared significant, not that 5% of the discoveries are false : the two answer different questions, and **controlling the false positive rate does not control the false discovery proportion**. > The following chunk runs the simulation for an increasing number of negatives (always with 100 positives), and plots the three quantities as a function of the proportion of negatives among all the genes. > Look at the three plots. What do they show? From which proportion of negatives are more than half of the discoveries false? ```{r} n_neg_values <- c(100, 200, 500, 1000, 2000, 5000, 10000) summary_1c <- lapply(n_neg_values, function(n_neg) { is_pos <- rep(c(TRUE, FALSE), c(n_pos, n_neg)) delta <- ifelse(is_pos, delta_pos, 0) evaluate_decision(simulate_pvalues(delta) < 0.05, is_pos) %>% select(power, FPR, FDP) %>% mutate(prop_neg = n_neg / (n_pos + n_neg)) }) %>% bind_rows() plot_power_1c <- ggplot(summary_1c, aes(x = prop_neg, y = power)) + geom_line() + geom_point() + ylim(0, 1) + labs(x = "Proportion of negatives among all the genes", y = "Power", title = "Power") plot_fpr_1c <- ggplot(summary_1c, aes(x = prop_neg, y = FPR)) + geom_line() + geom_point() + ylim(0, 1) + labs(x = "Proportion of negatives among all the genes", y = "False positive rate", title = "False positive rate") plot_fdp_1c <- ggplot(summary_1c, aes(x = prop_neg, y = FDP)) + geom_line() + geom_point() + ylim(0, 1) + labs(x = "Proportion of negatives among all the genes", y = "False discovery proportion", title = "False discovery proportion") ggsave(here("output", "power_by_proportion_negatives.pdf"), plot_power_1c, width = 6, height = 4) ggsave(here("output", "fpr_by_proportion_negatives.pdf"), plot_fpr_1c, width = 6, height = 4) ggsave(here("output", "fdp_by_proportion_negatives.pdf"), plot_fdp_1c, width = 6, height = 4) plot_power_1c plot_fpr_1c plot_fdp_1c ``` > **Answer:** The power and the false positive rate stay roughly constant (at almost 100% and about 5%, apart from some random variation) whatever the proportion of negatives, whereas the false discovery proportion increases steadily towards 100% as the proportion of negatives increases. With 100 positives, more than half of the discoveries are false once there are more than about 2,000 negatives (a proportion of negatives above about 95%), since 5% of 2,000 is 100 false positives for about 100 true positives. On the grid used here, the false discovery proportion is above 50% from a proportion of negatives of `r sprintf("%.0f", 100 * min(summary_1c$prop_neg[summary_1c$FDP > 0.5]))`%. In omics studies, we can often expect that a large number of variables are true negatives. For example, in transcriptomics, a large proportion of our genes have functions related to routine cellular processes which are not expected to be impacted by exposure. Therefore, in these cases, we can expect a large part of the genes declared significant at the 5% level to be false discoveries. ## 1d) Correcting for multiple testing : Bonferroni and Benjamini-Hochberg > We have just seen that testing each gene at the 5% level does not control the proportion of false discoveries. Two classic corrections address this problem, and they control two different error rates, defined from the numbers of false positives (FP) and true positives (TP) : > > - the **family-wise error rate** (FWER) : the probability of making *at least one* false positive, P(FP >= 1); > - the **false discovery rate** (FDR) : the *average* false discovery proportion, i.e. the average of FP / (FP + TP) over many repeated studies. > Note that these two error rates describe the behaviour of a procedure over many repeated studies. In a single study, we only observe the numbers FP and TP, and therefore the false discovery proportion of that study, and whether or not it contains at least one false positive. > We now simulate a study with **100 positives and 9,900 negatives** (10,000 genes), and look at what happens when we apply no correction (each p-value is compared to 0.05), the Bonferroni correction, and the Benjamini-Hochberg (BH) correction, one by one. Run the following chunk to simulate the p-values, and to look at the results with no correction. ```{r} n_pos <- 100 n_neg <- 9900 delta_pos <- 3 is_pos <- rep(c(TRUE, FALSE), c(n_pos, n_neg)) p_1d <- simulate_pvalues(ifelse(is_pos, delta_pos, 0)) # No correction : compare each p-value to 0.05 results_none <- evaluate_decision(p_1d < 0.05, is_pos) results_none ``` > **The Bonferroni correction** controls the FWER. It declares significant the genes whose p-value is below 0.05 / m, where m is the number of tests (here 0.05 / 10,000 = 0.000005). This is the same as multiplying each p-value by m (a "Bonferroni-adjusted p-value") and comparing it to 0.05. In R, the adjusted p-values are given by `p.adjust()`. ```{r} p_bonferroni <- p.adjust(p_1d, method = "bonferroni") results_bonferroni <- evaluate_decision(p_bonferroni < 0.05, is_pos) results_bonferroni ``` > **The Benjamini-Hochberg (BH) correction** controls the FDR. It ranks the p-values from the smallest to the largest, p(1) <= p(2) <= ... <= p(m). It finds the largest rank k such that p(k) <= (k / m) x 0.05, and declares significant the k genes with the smallest p-values. The threshold (k / m) x 0.05 increases with the rank, so for the first genes it is as strict as Bonferroni, but it becomes much more permissive as more discoveries are made. ```{r} p_bh <- p.adjust(p_1d, method = "BH") results_bh <- evaluate_decision(p_bh < 0.05, is_pos) results_bh ``` > Here are the three results side by side. ```{r} bind_rows( list("No correction" = results_none, "Bonferroni" = results_bonferroni, "Benjamini-Hochberg" = results_bh), .id = "method" ) %>% mutate(across(c(power, FPR, FDP), ~ round(.x, 3))) ``` > Compare the three results. What happens to the numbers of false positives (FP), true positives (TP) and false negatives (FN) with each correction, and to the three rates? What does each correction gain, and what does it cost? Which error rate does each of them control? > **Answer:** With no correction, there are many false positives (FP = `r results_none$FP`, so a false discovery proportion of `r sprintf("%.0f", 100 * results_none$FDP)`%) but all (or nearly all) the positives are found (TP = `r results_none$TP`, FN = `r results_none$FN`). The **Bonferroni** correction removes (nearly) all the false positives (FP = `r results_bonferroni$FP`) and controls the FWER, the probability of making at least one false positive. The price is a loss of power : the threshold 0.05 / m is so strict that some of the positives are missed (TP = `r results_bonferroni$TP` and FN = `r results_bonferroni$FN`, a power of `r sprintf("%.0f", 100 * results_bonferroni$power)`%). The **Benjamini-Hochberg** correction accepts a few false positives (FP = `r results_bh$FP`) in order to find many more positives than Bonferroni (TP = `r results_bh$TP`, FN = `r results_bh$FN`, a power of `r sprintf("%.0f", 100 * results_bh$power)`%). It controls the FDR : its false discovery proportion (`r sprintf("%.0f", 100 * results_bh$FDP)`% here) is, on average over many studies, at most 5%. Both corrections also lower the false positive rate, to almost 0. If you re-run the simulation, the numbers change, as each run is a new study, but the pattern stays the same. > Which correction would you choose in each of the following situations ? (a) A confirmatory study, in which the few genes found significant will be used to take a clinical decision. (b) An exploratory study of 20,000 genes, to find candidate genes to be followed up in other experiments. > **Answer:** (a) A correction which controls the FWER (Bonferroni), since a single false positive is costly : we want to be nearly sure that none of the genes is wrong. (b) A correction which controls the FDR (Benjamini-Hochberg), since the candidates will be validated afterwards : we accept a small proportion of false leads in exchange for finding more of the true genes. # 2) Application to real data : differential gene expression in tuberculosis > We now come back to real data. We use the transcriptomic data from [Singhania et al, 2018](https://www.nature.com/articles/s41467-018-04579-w) that you analysed in the genomics data analysis practical : blood samples from healthy controls, patients with latent tuberculosis infection (LTBI) and patients with active tuberculosis, collected in three cohorts. As in that practical, we compare the patients with **active tuberculosis** to the others (the controls and the LTBI patients grouped together), and we look for the genes that are **differentially expressed** between the two groups. > This is the same question as in the simulations of the first part, but this time we do not know the truth : for each of the thousands of genes, we do not know whether it is really differentially expressed. ## 2a) Differential expression analysis > The data file contains the data after the pre-processing steps of the previous practical : the raw counts after filtering and batch correction with ComBat-Seq (`counts_combat`, with one row per gene and one column per sample), and the information on the samples (`metadata_all`), including the variable `testGroup`, which groups the controls and the LTBI patients into one control group, and the active tuberculosis patients into the other group. Load the data. ```{r} load(here("data-raw", "TB_batch_corrected.RData")) # Number of genes (rows) and samples (columns), and number of samples in each group dim(counts_combat) table(metadata_all$testGroup) ``` > We perform the differential expression analysis with the voom-limma pipeline, as in the previous practical. For each gene, it tests whether its mean expression differs between the active tuberculosis patients and the control group, and returns a p-value. ```{r} # Design matrix : the coefficient "testGroupActiveTB" compares ActiveTB to control design <- model.matrix(~ testGroup, data = metadata_all) # TMM normalisation factors, then the mean-variance relationship with voom() dge <- DGEList(counts = counts_combat, group = metadata_all$testGroup) %>% calcNormFactors(method = "TMM") voom_data <- voom(dge, design) # Fit the linear models, then compute the empirical Bayes moderated tests fit <- lmFit(voom_data, design) fit <- eBayes(fit, robust = TRUE) # One row per gene, with the log2 fold-change (logFC) and the raw p-value (P.Value) results_dge <- topTable(fit, coef = 2, n = Inf) head(results_dge) ``` ## 2b) The p-values > Let's look at the raw p-values of all the genes. The following chunk plots their histogram (the red line marks 0.05), and gives the number of genes tested, the number of genes with a raw p-value below 0.05, and the number of such genes that we would expect by chance if no gene were differentially expressed. ```{r} n_tested <- nrow(results_dge) n_raw <- sum(results_dge$P.Value < 0.05) plot_pvalues_dge <- ggplot(results_dge, aes(x = P.Value)) + geom_histogram(breaks = seq(0, 1, 0.05), colour = "white") + geom_vline(xintercept = 0.05, colour = "red", linetype = "dashed") + labs(x = "p-value", y = "Number of genes", title = "Raw p-values of the differential expression analysis") ggsave(here("output", "pvalues_dge.pdf"), plot_pvalues_dge, width = 8, height = 5) plot_pvalues_dge c( "genes tested" = n_tested, "p-value below 0.05" = n_raw, "expected by chance if no gene was different" = 0.05 * n_tested ) ``` > Compare this histogram with those of the simulations of the first part. What do you notice? What do the flat part of the histogram and the peak near 0 represent? > **Answer:** The histogram looks like the mix of genes with and without a true difference that we simulated in the first part : there is a peak of small p-values, which comes from the genes that are truly differentially expressed (like the positives), on top of a flat part, which comes from the genes that are not differentially expressed (like the negatives, whose p-values are uniform). > How many genes have a raw p-value below 0.05? How many would we expect if none of the genes was differentially expressed? Looking at the flat part of the histogram (hint : the p-values above 0.5 come almost only from the genes that are not different), what proportion of the genes seems not to be differentially expressed? Using what you learned in 1c, what proportion of the genes declared significant at the 5% level (without correction) could be false positives? > **Answer:** `r n_raw` of the `r n_tested` genes have a raw p-value below 0.05, while we would only expect `r round(0.05 * n_tested)` if no gene was different. The p-values above 0.5 are almost all from genes that are not different, and they make up `r sprintf("%.0f", 100 * mean(results_dge$P.Value > 0.5))`% of the p-values, which covers half of the interval from 0 to 1. So the genes that are not different are about 2 x `r sprintf("%.0f", 100 * mean(results_dge$P.Value > 0.5))`% = `r sprintf("%.0f", 100 * min(1, 2 * mean(results_dge$P.Value > 0.5)))`% of all the genes. Among them, 5% have a p-value below 0.05, i.e. about `r round(0.05 * min(1, 2 * mean(results_dge$P.Value > 0.5)) * n_tested)` genes, so around `r sprintf("%.0f", 100 * 0.05 * min(1, 2 * mean(results_dge$P.Value > 0.5)) * n_tested / n_raw)`% of the `r n_raw` genes declared significant at the 5% level without correction may be false positives. This is the same phenomenon as in 1c : a large number of genes that are not different produces many false positives, even with a 5% threshold for each test. ## 2c) Correcting the p-values > Now apply the corrections. Using the function `p.adjust()`, compute the adjusted p-values of the genes with the Bonferroni, Holm, Benjamini-Hochberg (`method = "BH"`) and Benjamini-Yekutieli (`method = "BY"`) corrections, and store them in new columns of `results_dge`. For each correction, count the number of genes declared significant at the 5% level (adjusted p-value below 0.05), and compare it with the number of genes without correction (raw p-value below 0.05). ```{r} results_dge <- results_dge %>% mutate( p_bonferroni = p.adjust(P.Value, method = "bonferroni"), p_holm = p.adjust(P.Value, method = "holm"), p_bh = p.adjust(P.Value, method = "BH"), p_by = p.adjust(P.Value, method = "BY") ) n_discoveries <- c( "No correction" = sum(results_dge$P.Value < 0.05), "Bonferroni" = sum(results_dge$p_bonferroni < 0.05), "Holm" = sum(results_dge$p_holm < 0.05), "Benjamini-Hochberg" = sum(results_dge$p_bh < 0.05), "Benjamini-Yekutieli" = sum(results_dge$p_by < 0.05) ) n_discoveries ``` > Which corrections are the strictest and which are the most permissive? How do you explain the differences between Bonferroni and Holm, and between BH and BY? Among the genes declared significant with the BH correction, about how many do you expect to be false positives? > **Answer:** Without correction, `r n_discoveries["No correction"]` genes are declared significant, and the corrections reduce this number. The strictest are the two corrections that control the FWER, Bonferroni (`r n_discoveries["Bonferroni"]` genes) and Holm (`r n_discoveries["Holm"]` genes), which give very similar results (Holm is never stricter than Bonferroni, and slightly less strict). The corrections that control the FDR declare many more genes significant : Benjamini-Hochberg (`r n_discoveries["Benjamini-Hochberg"]` genes) and Benjamini-Yekutieli (`r n_discoveries["Benjamini-Yekutieli"]` genes). BY is more conservative than BH, because it is valid whatever the dependence between the tests (genes are often correlated) and divides the BH threshold by a factor c(m) = 1 + 1/2 + ... + 1/m, about `r round(sum(1 / seq_len(n_tested)), 1)` here. With the BH correction, the FDR is controlled at 5%, so among the `r n_discoveries["Benjamini-Hochberg"]` genes declared significant, we expect, on average, at most 5%, i.e. about `r round(0.05 * n_discoveries["Benjamini-Hochberg"])` false positives. In contrast, we expect several hundreds of false positives without correction (see 2b). ## 2d) How do the corrections work? > Each correction compares the sorted p-values $p_{(1)} \le p_{(2)} \le \dots \le p_{(m)}$ (where $m$ is the number of genes tested, and $i$ is the rank of a p-value) to a threshold that depends on the rank $i$ and on the significance level $\alpha = 0.05$ : > > - Bonferroni : $p_{(i)} \le \dfrac{\alpha}{m}$ > - Holm : $p_{(i)} \le \dfrac{\alpha}{m - i + 1}$ > - Benjamini-Hochberg : $p_{(i)} \le \dfrac{i}{m}\,\alpha$ > - Benjamini-Yekutieli : $p_{(i)} \le \dfrac{i}{m}\,\dfrac{\alpha}{c(m)}$, with $c(m) = \sum_{j=1}^{m} \dfrac{1}{j}$ > The following chunk plots the sorted p-values against their rank (black points), on a log scale. The coloured lines are the thresholds of the four corrections. The dashed vertical lines show the number of genes declared significant by each correction : each correction declares significant all the genes up to a given rank. (We only show the first ranks.) ```{r} alpha <- 0.05 rank <- seq_len(n_tested) # The sorted p-values pvalues_sorted <- data.frame(rank = rank, p = sort(results_dge$P.Value)) # The thresholds of the four corrections, for each rank c_m <- sum(1 / rank) thresholds <- data.frame( rank = rank, "Bonferroni" = alpha / n_tested, "Holm" = alpha / (n_tested - rank + 1), "Benjamini-Hochberg" = rank / n_tested * alpha, "Benjamini-Yekutieli" = rank / n_tested * alpha / c_m, check.names = FALSE ) %>% pivot_longer(-rank, names_to = "correction", values_to = "threshold") # The number of genes declared significant by each correction n_rejected <- data.frame( correction = c("Bonferroni", "Holm", "Benjamini-Hochberg", "Benjamini-Yekutieli"), n = c(sum(p.adjust(results_dge$P.Value, "bonferroni") < alpha), sum(p.adjust(results_dge$P.Value, "holm") < alpha), sum(p.adjust(results_dge$P.Value, "BH") < alpha), sum(p.adjust(results_dge$P.Value, "BY") < alpha)) ) plot_thresholds <- ggplot() + geom_point(data = pvalues_sorted, aes(x = rank, y = p), size = 0.5) + geom_line(data = thresholds, aes(x = rank, y = threshold, colour = correction), linewidth = 0.8) + geom_vline(data = n_rejected, aes(xintercept = n, colour = correction), linetype = "dashed") + scale_y_log10() + coord_cartesian(xlim = c(0, 1.5 * max(n_rejected$n, 100)), ylim = c(1e-8, 1)) + labs(x = "Rank of the p-value", y = "p-value (log scale)", colour = "Correction", title = "Sorted p-values and thresholds of the corrections") ggsave(here("output", "pvalues_thresholds.pdf"), plot_thresholds, width = 9, height = 6) plot_thresholds ``` > Describe how the threshold changes with the rank for each correction. Which is the strictest for the first ranks? Which correction declares significant the most genes? > **Answer:** The Bonferroni threshold is the same for all ranks (a flat line, $\alpha / m$). The Holm threshold starts at the same value, and increases with the rank, slowly at first and then very quickly for the last ranks. The Benjamini-Hochberg and Benjamini-Yekutieli thresholds increase linearly with the rank (straight lines), BY being lower than BH by a factor $c(m)$. For the first ranks, the strictest are Bonferroni and Holm for the first rank, and BY (a threshold of $\alpha / (m\,c(m))$ for the first rank, below the Bonferroni threshold). As the rank increases, the BH threshold becomes much higher than the Bonferroni and Holm thresholds, so BH declares significant the most genes (`r n_rejected$n[n_rejected$correction == "Benjamini-Hochberg"]` genes).