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 :

In RStudio, click File > New Project... > New Directory > New Project.

Give the project a name, for example omics_practicals, and choose where to save it on your computer. Click Create Project.

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.

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

Copy the data file I sent you, TB_batch_corrected.RData, into the data-raw folder (we will use it in the second part).

Copy this practical (the .Rmd file) into the analysis folder.

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).

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.

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.

# Path to the root of the project
here()

## [1] "/home/ah3/Desktop/Work/Teaching/PHDS_omics_data_course"

# Check that the data file exists
file.exists(here("data-raw", "TB_batch_corrected.RData"))

## [1] TRUE

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.

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))

##    control    treated 
## 0.07462564 0.20862196

t.test(control, treated, var.equal = TRUE)$p.value

## [1] 0.7671784

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.

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?

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

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?

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?

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

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?

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

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.

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?

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"))
)

##           Decision
## Truth      Significant Not significant
##   Positive         100               0
##   Negative           9              91

# The numbers TP, FP, FN, TN and the three rates
results_mix <- evaluate_decision(significant, is_pos)
results_mix

##    TP FP FN TN power  FPR        FDP
## 1 100  9  0 91     1 0.09 0.08256881

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.

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?

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

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.

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

##    TP  FP FN   TN power        FPR       FDP
## 1 100 490  0 9410     1 0.04949495 0.8305085

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().

p_bonferroni <- p.adjust(p_1d, method = "bonferroni")

results_bonferroni <- evaluate_decision(p_bonferroni < 0.05, is_pos)
results_bonferroni

##   TP FP FN   TN power FPR FDP
## 1 58  0 42 9900  0.58   0   0

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.

p_bh <- p.adjust(p_1d, method = "BH")

results_bh <- evaluate_decision(p_bh < 0.05, is_pos)
results_bh

##   TP FP FN   TN power          FPR        FDP
## 1 97  7  3 9893  0.97 0.0007070707 0.06730769

Here are the three results side by side.

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)))

##               method  TP  FP FN   TN power   FPR   FDP
## 1      No correction 100 490  0 9410  1.00 0.049 0.831
## 2         Bonferroni  58   0 42 9900  0.58 0.000 0.000
## 3 Benjamini-Hochberg  97   7  3 9893  0.97 0.001 0.067

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?

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.

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 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.

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)

## [1] 19387   239

table(metadata_all$testGroup)

## 
##  Control ActiveTB 
##      202       37

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.

# 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)

##                    logFC     AveExpr        t      P.Value    adj.P.Val
## ENSG00000152766 5.138118  1.56404852 20.76697 1.313894e-55 2.547246e-51
## ENSG00000108387 3.507219  2.30138267 20.03362 2.188978e-53 2.121886e-49
## ENSG00000183347 3.936295 -0.25939566 19.57237 7.284027e-52 4.707181e-48
## ENSG00000108950 3.544596  0.38022211 17.90084 2.771553e-46 1.343303e-42
## ENSG00000197646 4.065395 -0.21780826 16.68325 3.602383e-42 1.396788e-38
## ENSG00000100336 3.790662 -0.04602701 16.52300 1.259232e-41 4.068788e-38
##                         B
## ENSG00000152766 115.10172
## ENSG00000108387 110.34229
## ENSG00000183347 105.80284
## ENSG00000108950  93.71105
## ENSG00000197646  84.32271
## ENSG00000100336  83.14660

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.

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
)

##                                genes tested 
##                                    19387.00 
##                          p-value below 0.05 
##                                    11918.00 
## expected by chance if no gene was different 
##                                      969.35

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?

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?

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).

# Your code here

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?

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.)

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?