Before we start

0.1 Introduction

In Lab 6 we downloaded the GSE161731 blood study, gave the genes readable names, and cleaned the patient table. Today we ask the biological question: which genes are switched on or off in the blood of COVID-19 patients compared with healthy people?

We do it in two steps:

  1. PCA — Explore sample variation. First we make samples comparable (Counts Per Million: CPM), then we draw a “map” of all samples (Principal Component Analysis: PCA) to see which samples resemble each other, and whether the biggest difference in the data is the disease.
  2. DESeq2 — Identify differentially expressed genes. We compare COVID-19 with healthy for every gene, correct for the fact that we test ~18,000 genes at once, and show the result as a volcano plot, single-gene plots and a heatmap.

Our positive control from Lab 6 still holds: if the analysis is right, the interferon genes (IFI27, ISG15, RSAD2) should come out as higher in COVID-19.

At the end there is an Extra section on pathway analysis (ORA and GSEA). It is not part of the assignments — it is there for those who want to go further.

How each part works (same as Lab 6): an Exercise with questions 1–3, a Solution tab, and then an Assignment that is almost the same as the exercise with one small change.

Part Topic Ends with
1 PCA: what does the data look like before we model it? Exercise 1 → Assignment 15
2 DESeq2: which genes differ between COVID-19 and healthy? Exercise 2 → Assignment 16
Extra Pathway analysis: ORA and GSEA (optional) Try-it-yourself questions

0.2 Packages

Tip: If you see Update all/some/none? [a/s/n] in the RStudio Console, type a and press Enter to update all outdated packages and continue installation (or n to skip updates).

Always check that your previous command has finished before running another: the > prompt means R is ready, while + means R is waiting for you to complete an unfinished command.

# Run once.
BiocManager::install(c("DESeq2", "edgeR"))
install.packages(c("ggplot2", "dplyr", "ggrepel", "pheatmap", "RColorBrewer"))
library(DESeq2)         # differential expression
library(edgeR)          # cpm() and filterByExpr()
library(ggplot2)        # plots
library(dplyr)
library(ggrepel)        # gene labels that do not overlap
library(pheatmap)       # heatmaps
library(RColorBrewer)   # colour palettes

theme_set(theme_bw(base_size = 12))
data_dir <- "data"; res_dir <- "results"
dir.create(res_dir, showWarnings = FALSE)

0.3 Loading the objects from Lab 6

counts <- read.csv(file.path(data_dir, "gse161731_counts.csv"), row.names = 1)           # raw counts
meta   <- read.csv(file.path(data_dir, "gse161731_meta.csv"), row.names = 1)             # clean patient table

dim(counts)
## [1] 35447   177
table(meta$cohort)
## 
## Bacterial CoV_other  COVID_19   healthy Influenza 
##        24        40        77        19        17
# counts <- counts[,rownames(meta)]
all(colnames(counts) == rownames(meta))   # the order check, one more time
## [1] TRUE

1 PCA: what does the data look like before we model it?

1.1 CPM first: make the samples comparable

Imagine two blood samples with exactly the same biology, but one was sequenced twice as deeply. Every gene in that sample gets about twice as many reads — not because of biology, but because of the machine.

CPM (counts per million) fixes this: for each sample, divide each gene’s count by the sample’s total reads, and multiply by one million. The question changes from “how many reads?” to “how many reads out of every million?”

lib_size <- colSums(counts)          # total reads per sample
range(lib_size) / 1e6                # smallest and largest, in millions
## [1]   1.501329 128.474589
# CPM by hand for one gene in the first 5 samples:
counts["IFI27", 1:5] / lib_size[1:5] * 1e6

1.1.1 Remove genes that are (almost) never expressed

Many genes are silent in blood. They carry no information, so we remove them. filterByExpr() keeps a gene if it has a reasonable number of reads in at least as many samples as the smallest group.

keep_gene <- filterByExpr(counts, group = meta$cohort)
table(keep_gene)                     # FALSE = removed, TRUE = kept
## keep_gene
## FALSE  TRUE 
## 17467 17980
counts_f <- counts[keep_gene, ]      # the filtered raw counts (used for DESeq2 in Part 2)
dim(counts_f)
## [1] 17980   177

1.1.2 log2 CPM

A few genes have enormous counts and most have small ones. Taking log2 squeezes this range so that “twice as much” is always one step up, whether a gene goes from 10 to 20 or from 1000 to 2000.

logcpm <- cpm(counts_f, log = TRUE, prior.count = 1)   # log2 CPM (a small number is added so log2(0) is not a problem)
logcpm[1:5, 1:4]
##        SSS94189 SSS97389 SSSDU09.03S19498 SSS94478
## ACTB   13.47764 12.97710         14.03524 13.20341
## MT-CO1 13.06796 11.85181         13.04154 11.64574
## HLA-B  12.85560 12.66163         11.04469 12.57538
## HLA-C  12.56977 12.40191         10.91248 12.44266
## CSF3R  12.97123 13.55698         12.14039 13.31810
par(mfrow = c(1, 2))                                    # two plots side by side
hist(cpm(counts_f)[, 1], breaks = 50, main = "CPM", xlab = "CPM")
hist(logcpm[, 1],        breaks = 50, main = "log2 CPM", xlab = "log2 CPM")

par(mfrow = c(1, 1))

1.2 PCA: a map of the samples

Each sample has thousands of gene values — impossible to plot directly. PCA summarises them into a few new axes:

  • PC1 is the direction in which the samples differ the most;
  • PC2 is the second biggest difference, and so on.

Samples close together on the PCA plot have similar gene expression. The key question: is the biggest difference in my data the one I care about (the disease)?

# 1. Use the 500 genes that vary most between samples (they carry the signal).
gene_var  <- apply(logcpm, 1, var)                          # variance of each gene
top_genes <- names(sort(gene_var, decreasing = TRUE))[1:500]
mat       <- logcpm[top_genes, ]

# 2. Run the PCA. prcomp() wants samples in rows, so we flip the table with t().
pca <- prcomp(t(mat), center = TRUE, scale. = FALSE)

# 3. How much of the total difference does each PC capture (in %)?
var_pct <- round(100 * pca$sdev^2 / sum(pca$sdev^2), 1)
var_pct[1:6]
## [1] 13.4  9.6  6.3  5.1  3.8  2.5
# Scree plot: the % captured by each PC.
barplot(var_pct[1:10], names.arg = paste0("PC", 1:10),
        ylab = "% variance explained", main = "Scree plot")

# Put the sample positions (pca$x) together with the patient table.
pca_df <- data.frame(PC1 = pca$x[, 1], PC2 = pca$x[, 2], meta)
pca_df$cohort <- factor(pca_df$cohort, levels = c("healthy", "Bacterial", "CoV_other", "COVID_19", "Influenza"))

ggplot(pca_df, aes(x = PC1, y = PC2, colour = cohort)) +
  geom_point(size = 2.5, alpha = 0.85) +
  scale_colour_brewer(palette = "Set2") +
  stat_ellipse(
    aes(group = cohort),
    type = "norm",
    level = 0.95,
    linewidth = 0.8
  )+
  labs(x = paste0("PC1 (", var_pct[1], "%)"),
       y = paste0("PC2 (", var_pct[2], "%)"),
       title = "PCA of the 500 most variable genes")

What you should see. Bacterial infection gives a strong, distinct response (many neutrophils), so bacterial samples tend to separate. COVID-19, influenza and the other coronaviruses all trigger a similar interferon response, so they overlap. Healthy and mild viral samples can also overlap. PCA shows only the biggest differences; the finer ones need DESeq2 (Part 2).

stat_ellipse() adds a separate ellipse for each cohort to show the distribution of samples. - aes(group = cohort): Draws a separate ellipse for each group. - type = "norm": Assumes a multivariate normal distribution. - level = 0.95: Sets the ellipse level to 95%. - linewidth = 0.8: Controls the thickness of the ellipse outline.

Interpretation: Ellipses help visualize within-group variability and overlap between cohorts in PCA space. They do not >indicate statistical significance.

1.3 Which variable drives which PC?

Colour the same map by other variables — and test it, rather than guessing.

ggplot(pca_df, aes(PC1, PC2, colour = age)) +
  geom_point(size = 2.5) +
  scale_colour_viridis_c() +
  labs(title = "Same PCA, coloured by age")

ggplot(pca_df, aes(PC1, PC2, colour = gender)) +
  geom_point(size = 2.5) +
  scale_colour_brewer(palette = "Set3") +
    stat_ellipse(
    aes(group = gender),
    type = "norm",
    level = 0.95,
    linewidth = 0.8
  )+
  labs(title = "Same PCA, coloured by sex")

# Which genes push samples along PC1? These are the "loadings".
load1 <- sort(pca$rotation[, 1])
head(load1, 10)     # genes pulling samples to the left
##       CD177        MMP8       TMT1B     ADAMTS2      SLC51A     PCOLCE2      FAM20A        OLAH 
## -0.13814574 -0.11996857 -0.11617768 -0.10800190 -0.10497790 -0.10033004 -0.09908607 -0.09696268 
##      NECAB1   TNFAIP8L3 
## -0.08863570 -0.08840498
tail(load1, 10)     # genes pulling samples to the right
##    TRBV6-3      LRRN3      MXRA8      TRDV2     CACNG6      OLIG2        NOG    SIGLEC8     ALOX15 
## 0.08476489 0.08787174 0.08807177 0.08928849 0.08942153 0.09840933 0.10181153 0.11960078 0.13181815 
##     PRSS33 
## 0.13545180

Interpretation: PC1 appears to reflect differences in immune cell composition. For example, CD177 and MMP8 are associated with neutrophils, while TRBV6-3 and TRDV2 are related to T cells. This suggests that differences in blood cell populations may contribute to the distinct gene expression profiles between COVID-19 and bacterial infection.

1.4 Exercise 1

1.4.1 Question

  1. Redo the PCA with the 2000 most variable genes. What % does PC1 explain now?
  2. Draw PC1 as a box plot per cohort. Which cohorts separate?
  3. Is IFI27 among the 20 genes with the largest PC1 loadings?

1.4.2 Solution

# 1  Only the number 500 changed to 2000.
top2k <- names(sort(gene_var, decreasing = TRUE))[1:2000]
pca2  <- prcomp(t(logcpm[top2k, ]), center = TRUE, scale. = FALSE)
round(100 * pca2$sdev^2 / sum(pca2$sdev^2), 1)[1:3]
## [1] 16.2  7.6  3.6
# 2
ggplot(pca_df, aes(x = cohort, y = PC1, fill = cohort)) +
  geom_boxplot(alpha = 0.7, outlier.shape = NA) +
  geom_jitter(width = 0.15, size = 0.8) +
  scale_fill_brewer(palette = "Set2") +
  theme(legend.position = "none")

# 3  abs(): ignore the sign, we want the biggest pushes in either direction.
top_load <- names(sort(abs(pca$rotation[, 1]), decreasing = TRUE))[1:20]
"IFI27" %in% top_load
## [1] FALSE
top_load
##  [1] "CD177"     "PRSS33"    "ALOX15"    "MMP8"      "SIGLEC8"   "TMT1B"     "ADAMTS2"   "SLC51A"   
##  [9] "NOG"       "PCOLCE2"   "FAM20A"    "OLIG2"     "OLAH"      "CACNG6"    "TRDV2"     "NECAB1"   
## [17] "TNFAIP8L3" "MXRA8"     "LRRN3"     "KLF14"

1.5 Assignment 15 — PCA

Hint: start from the code of Exercise 1 and change one small thing.

  1. Redo the PCA with the 1000 most variable genes and report the % explained by PC1 and PC2.
  2. Draw PC2 as a box plot per cohort. Which cohorts separate on PC2?
  3. Is S100A8 among the 50 genes with the largest PC1 loadings?

2 DESeq2: which genes differ between COVID-19 and healthy?

2.1 The comparison, and the design formula

We compare two groups: COVID-19 and healthy. DESeq2 always works on the raw counts — it does its own normalisation inside, so we give it counts with raw integer read counts, not CPM.

# Keep the two groups we want to compare.
meta_de   <- meta[meta$cohort %in% c("healthy", "COVID_19"), ]
counts_de <- counts_f[, rownames(meta_de)]               # the same samples, same order

table(meta_de$cohort)
## 
## COVID_19  healthy 
##       77       19
all(rownames(meta_de) == colnames(counts_de)) # always check
## [1] TRUE
# Age (number) vs group -> Wilcoxon test
tapply(meta_de$age, meta_de$cohort, median)   # median age in each group
## COVID_19  healthy 
##       33       18
wilcox.test(age ~ cohort, data = meta_de) # p-value = 1.762e-11, age is significant between healthy and COVID_19
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  age by cohort
## W = 1461.5, p-value = 1.762e-11
## alternative hypothesis: true location shift is not equal to 0
# Sex (category) vs group -> Fisher 
tab2 <- table(meta_de$gender, meta_de$cohort)
tab2
##         
##          COVID_19 healthy
##   Female       32       8
##   Male         45      11
fisher.test(tab2) # p-value = 1, gender is NOT significant between healthy and COVID_19
## 
##  Fisher's Exact Test for Count Data
## 
## data:  tab2
## p-value = 1
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.3163508 3.1428759
## sample estimates:
## odds ratio 
##  0.9779987

The design formula tells DESeq2 which factors to consider when comparing gene expression between samples.

  • ~ cohort: Tests for differences in gene expression between cohorts, such as COVID-19 and bacterial infection.
  • ~ cohort + age + batch: Also accounts for differences in age and technical batch effects, helping us estimate the cohort effect more fairly. Including age and batch helps reduce confounding, so observed expression differences are less likely to be explained by these factors instead of the disease group.

Important: Only include variables when the data support it. For example, if cohort and batch are completely confounded (each cohort was processed in a separate batch), their effects cannot be reliably separated. Also, DESeq2 requires the variables in the design formula to be present in the sample metadata.

dds <- DESeqDataSetFromMatrix(
  countData = counts_de,   # raw counts — never CPM or normalised data
  colData   = meta_de,     # the patient table, same order as the columns
  design    = ~ cohort+age+batch     # the comparison: Adjusting for age and batch
)

dds <- dds[rowSums(counts(dds)) >= 10, ]   # drop genes with fewer than 10 reads in total
dds
## class: DESeqDataSet 
## dim: 17980 96 
## metadata(1): version
## assays(1): counts
## rownames(17980): ACTB MT-CO1 ... SNORD62A MIR5192
## rowData names(0):
## colnames(96): SSSDU18.02S0011619 SSSDU18.02S0011620 ... SSSDU09.02S0000156
##   SSSDU09.02S0000153
## colData names(8): subject_id age ... hospitalized batch
# DESeq() does everything in one line:
#   1. corrects for sequencing depth (like CPM, but more robust)
#   2. estimates how variable each gene is between patients of the same group
#   3. tests every gene: is it different between COVID-19 and healthy?
dds <- DESeq(dds)

# Pull out the comparison: c(variable, group of interest, reference group)
res_dds <- results(dds, contrast = c("cohort", "COVID_19", "healthy"), alpha = 0.05)
summary(res_dds)
## 
## out of 17980 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 822, 4.6%
## LFC < 0 (down)     : 1723, 9.6%
## outliers [1]       : 159, 0.88%
## low counts [2]     : 4790, 27%
## (mean count < 45)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
res <- as.data.frame(res_dds)
head(res )

What the columns mean:

Column Meaning in plain words
baseMean average expression of the gene over all samples
log2FoldChange how much the gene changes: +1 = 2× higher in COVID-19, −1 = 2× lower, 0 = no change
lfcSE how uncertain that change is
stat change divided by uncertainty
pvalue the p-value of this one gene
padj the adjusted p-value — the one to use (see 2.2)

2.2 Multiple testing

We tested around 18,000 genes simultaneously. If we use a raw p-value cut-off of 0.05, about 5% of genes with no real change may appear significant by chance — potentially around 900 false positives.

This is why we need multiple-testing correction: the more genes we test, the greater the chance of finding false positives just by chance. DESeq2 adjusts p-values using the Benjamini–Hochberg method to control the false discovery rate (FDR).

The adjusted p-value (padj) corrects for this. Using padj < 0.05 means we aim to keep the expected proportion of false discoveries among the genes we call significant at around 5%, under the method’s assumptions. It does not guarantee that exactly 5% of our selected genes are false positives.

sum(res$pvalue < 0.05, na.rm = TRUE)   # "significant" by raw p-value
## [1] 4457
sum(res$padj   < 0.05, na.rm = TRUE)   # significant after correction (fewer, but trustworthy)
## [1] 2545

Why do some genes have padj = NA? DESeq2 does not waste tests on genes with too few reads to ever reach significance. They get padj = NA. That is information, not an error.

sum(is.na(res$padj))              # genes without an adjusted p-value
## [1] 4949
metadata(res_dds)$filterThreshold     # the expression level below which genes were skipped
## 27.14286% 
##  44.53633

2.3 Reading and visualising the result

2.3.1 The results table

res <- res[order(res$padj), ]     # most significant genes first
res$gene <- rownames(res)
head(res[, c("gene", "baseMean", "log2FoldChange", "padj")], 15)
# The positive control: where are the known interferon genes?
check <- c("IFI27", "IFI44L", "ISG15", "RSAD2", "OAS1", "MX1", "CD3D", "S100A8")
res[res$gene %in% check, c("gene", "log2FoldChange", "padj")]

2.3.2 Volcano plot

Each dot represents one gene in a volcano plot.

  • X-axis (log₂ fold change): Shows how much gene expression changes between the two groups. Positive values mean higher expression in one group; negative values mean higher expression in the other.
  • Y-axis (−log₁₀ adjusted p-value): Shows statistical significance. The higher a dot is, the stronger the evidence that the gene is differentially expressed.
  • Top-left and top-right corners: These contain genes with statistically significant changes and large expression differences, making them especially interesting candidates for further investigation.

Genes near the centre have smaller expression changes, while genes near the bottom have weaker statistical evidence. Remember to check both the fold change and adjusted p-value when identifying significant genes.

vol <- res[!is.na(res$padj), ]            # drop genes without padj

# A "status" column, only to colour the dots.
vol$status <- "not significant"
vol$status[vol$padj < 0.05 & vol$log2FoldChange >  1] <- "up in COVID-19"
vol$status[vol$padj < 0.05 & vol$log2FoldChange < -1] <- "down in COVID-19"
table(vol$status)
## 
## down in COVID-19  not significant   up in COVID-19 
##              410            12579               42
# the 20 most significant genes get a label
top_lab <- rbind(
  head(vol[vol$log2FoldChange > 1, ][order(vol[vol$log2FoldChange > 0, ]$padj), ], 10),
  head(vol[vol$log2FoldChange < -1, ][order(vol[vol$log2FoldChange < 0, ]$padj), ], 10)
)   

top_lab
ggplot(vol, aes(x = log2FoldChange, y = -log10(padj), colour = status)) +
  geom_point(size = 1, alpha = 0.6) +
  scale_colour_manual(values = c("up in COVID-19"   = "firebrick",
                                 "down in COVID-19" = "steelblue",
                                 "not significant"  = "grey80")) +
  geom_vline(xintercept = c(-1, 1), linetype = 2, colour = "grey40") +  # 2x change lines
  geom_hline(yintercept = -log10(0.05), linetype = 2, colour = "grey40") + # padj = 0.05 line
  geom_text_repel(data = top_lab, aes(label = gene), size = 3, max.overlaps = 20) +
  labs(x = "log2 fold change (COVID-19 / healthy)", y = "-log10 adjusted p",
       colour = NULL, title = "COVID-19 vs healthy, whole blood")

2.3.3 One gene at a time

A single-gene plot is the best way to convince yourself (and a reviewer) that a change is real: one dot per patient.

plotCounts() shows the normalised read counts of a selected gene across individual samples, grouped by cohort. It helps us visualise differences in gene expression and variation between samples.

best <- res$gene[1]    # the most significant gene
plotCounts(dds, gene = best, intgroup = "cohort", main = best)

2.3.4 Heatmap of the top genes

pheatmap() visualises the expression levels of multiple genes across samples using colours to represent relative expression. It helps us identify patterns of gene expression and similarities or differences between samples.

top50 <- head(res$gene, 50)                 # the 50 most significant genes
mat50 <- logcpm[top50, colnames(counts_de)]    # their log2 CPM in these samples

pheatmap(mat50,
         scale          = "row",   # compare each gene with itself: red = above its average
         annotation_col = data.frame(cohort = meta_de$cohort, row.names = rownames(meta_de)),
         show_colnames  = FALSE,
         fontsize_row   = 7,
         cluster_rows   = TRUE,    # group genes with similar patterns
         cluster_cols   = FALSE,   # keep healthy and COVID-19 side by side
         color          = colorRampPalette(rev(brewer.pal(9, "RdBu")))(100),
         main           = "Top 50 differentially expressed genes")

2.4 Exercise 2

2.4.1 Question

  1. How many genes have padj < 0.05, and how many of those also change more than 2-fold (|log2FC| > 1)?
  2. Redraw the volcano plot, labelling only the top 10 genes and colouring up-regulated genes "red" and down-regulated genes "blue".
  3. Make a plotCounts() figure for IFI27. Does the direction fit a viral infection?

2.4.2 Solution

# 1  abs() = size of the change, ignoring direction
sum(res$padj < 0.05, na.rm = TRUE)
## [1] 2545
sum(res$padj < 0.05 & abs(res$log2FoldChange) > 1, na.rm = TRUE)
## [1] 452
# 2  Changed: head(vol, 10) and the two colours.
top_lab10 <- rbind(
  head(vol[vol$log2FoldChange > 1, ][order(vol[vol$log2FoldChange > 0, ]$padj), ], 5),
  head(vol[vol$log2FoldChange < -1, ][order(vol[vol$log2FoldChange < 0, ]$padj), ], 5)
)   

ggplot(vol, aes(x = log2FoldChange, y = -log10(padj), colour = status)) +
  geom_point(size = 1, alpha = 0.6) +
  scale_colour_manual(values = c("up in COVID-19"   = "red",
                                 "down in COVID-19" = "blue",
                                 "not significant"  = "grey80")) +
  geom_vline(xintercept = c(-1, 1), linetype = 2, colour = "grey40") +
  geom_hline(yintercept = -log10(0.05), linetype = 2, colour = "grey40") +
  geom_text_repel(data = top_lab10, aes(label = gene), size = 3) +
  labs(x = "log2 fold change (COVID-19 / healthy)", y = "-log10 adjusted p",
       colour = NULL, title = "COVID-19 vs healthy, whole blood")

# 3  IFI27 is higher in COVID-19: a typical interferon response to a virus.
plotCounts(dds, gene = "IFI27", intgroup = "cohort", main = "IFI27")

2.5 Assignment 16 — Differential expression

Hint: start from the code of Exercise 2 and change one small thing.

  1. How many genes have padj < 0.01, and how many of those also change more than 4-fold (|log2FC| > 2)?
  2. Redraw the volcano plot, labelling only the top 30 genes and colouring up-regulated genes "darkorange" and down-regulated genes "purple".
  3. Make a plotCounts() figure for CD3D (a T-cell gene).

3 Save the results

write.csv(res, file.path(res_dir, "DE_COVID19_vs_healthy.csv"), row.names = FALSE)

Take-home messages

  1. Raw counts depend on sequencing depth; CPM makes samples comparable for plots.
  2. Look before you test: PCA shows whether the biggest difference in your data is the biology you care about, or something else (cell composition, age, batch).
  3. DESeq2 needs raw counts and a clear comparison: which group is the reference.
  4. Testing thousands of genes creates false positives by chance — use padj.
  5. Check a positive control (here: interferon genes up in COVID-19) before you believe anything new.

Extra — Pathway analysis: ORA and GSEA (optional)

This section is not part of the assignments. It is here if you want to go one step further: from a list of genes to the biological processes behind them.

A list of 1000 genes is hard to read. Pathway analysis asks: are genes of a known process (e.g. “interferon response”) over-represented among my changed genes?

  • ORA (over-representation analysis): take the significant genes, and count how many belong to each pathway compared with chance.
  • GSEA (gene set enrichment analysis): rank all genes from most up to most down, and ask whether a pathway’s genes are bunched at the top or the bottom.

By default the code below is not run when you knit (it needs extra packages and an internet connection). You can still run the chunks one by one, or set run_extra <- TRUE below.

run_extra <- FALSE   # change to TRUE to run the Extra section when knitting
BiocManager::install(c("clusterProfiler", "enrichplot", "org.Hs.eg.db"))
install.packages("msigdbr")
library(clusterProfiler)   # ORA and GSEA
library(enrichplot)        # dotplot(), gseaplot2()
library(org.Hs.eg.db)      # gene ID conversion
library(msigdbr)           # MSigDB Hallmark gene sets

3.1 ORA — over-representation analysis

ORA needs three things:

  1. the gene list (significant genes, up and down separately);
  2. the gene sets (pathways);
  3. the background — all genes that were actually tested, not all human genes (many genes are never expressed in blood, and including them makes everything look falsely significant).
background <- res$gene[!is.na(res$pvalue)]   # every gene that was tested

genes_up   <- res$gene[!is.na(res$padj) & res$padj < 0.05 & res$log2FoldChange >  1]
genes_down <- res$gene[!is.na(res$padj) & res$padj < 0.05 & res$log2FoldChange < -1]
length(background); length(genes_up); length(genes_down)

3.1.1 ORA against MSigDB Hallmark

The 50 Hallmark gene sets are short, curated summaries of well-known biological states — a good first choice.

hallmark <- msigdbr(species = "Homo sapiens", collection = "H")
# hallmark <- msigdbr(species = "Homo sapiens", category = "H")   # older msigdbr versions

h_sets <- as.data.frame(hallmark[, c("gs_name", "gene_symbol")])   # pathway name, gene
length(unique(h_sets$gs_name))
ora_h_up <- enricher(gene = genes_up, universe = background, TERM2GENE = h_sets,
                     pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.2)
head(as.data.frame(ora_h_up)[, c("ID", "GeneRatio", "BgRatio", "p.adjust", "Count")], 10)

# Dot size = number of our genes in the pathway; colour = adjusted p-value.
dotplot(ora_h_up, showCategory = 10, title = "ORA, Hallmark — up in COVID-19")
ora_h_down <- enricher(gene = genes_down, universe = background, TERM2GENE = h_sets,
                       pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.2)
head(as.data.frame(ora_h_down)[, c("ID", "GeneRatio", "BgRatio", "p.adjust", "Count")], 10)

3.1.2 ORA against KEGG

KEGG uses Entrez gene IDs, so we translate the symbols first.

map_up <- bitr(genes_up,   fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)
map_bg <- bitr(background, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)

ora_kegg_up <- enrichKEGG(gene = map_up$ENTREZID, universe = map_bg$ENTREZID,
                       organism = "hsa", pAdjustMethod = "BH", pvalueCutoff = 0.05)
ora_kegg_up <- setReadable(ora_kegg_up, OrgDb = org.Hs.eg.db, keyType = "ENTREZID")  # IDs back to symbols
head(as.data.frame(ora_kegg_up)[, c("Description", "p.adjust", "Count")], 10)

barplot(ora_kegg_up, showCategory = 10, title = "ORA, KEGG — up in COVID-19")

as.data.frame(ora_kegg_up)$geneID[1]   # which of our genes drove the top pathway?
map_down <- bitr(genes_down,   fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)
map_bg <- bitr(background, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)

ora_kegg_down <- enrichKEGG(gene = map_down$ENTREZID, universe = map_bg$ENTREZID,
                       organism = "hsa", pAdjustMethod = "BH", pvalueCutoff = 0.05)
ora_kegg_down <- setReadable(ora_kegg_down, OrgDb = org.Hs.eg.db, keyType = "ENTREZID")  # IDs back to symbols
head(as.data.frame(ora_kegg_down)[, c("Description", "p.adjust", "Count")], 10)

barplot(ora_kegg_down, showCategory = 10, title = "ORA, KEGG — down in COVID-19")

as.data.frame(ora_kegg_down)$geneID[1]   # which of our genes drove the top pathway?

Pitfalls. (1) The result depends on your cut-offs. (2) A gene missing from the databases can never be found. (3) KEGG contains a “Coronavirus disease – COVID-19” pathway built partly from studies like this one; finding it is reassuring, not a discovery.

3.2 GSEA — gene set enrichment analysis

ORA throws away every gene below the cut-off. GSEA keeps all genes, ranked from most up to most down in COVID-19. A pathway where 200 genes each go up a little — none significant on its own — can still be detected. That is typical for immune programmes in blood.

We rank by stat (change ÷ uncertainty): a big change measured precisely ranks high; a big change measured on very few reads does not.

rank_df <- as.data.frame(res)
rank_df$gene <- rownames(rank_df)
rank_df <- rank_df[!is.na(rank_df$stat), ]   # GSEA cannot use missing values

ranked <- rank_df$stat
names(ranked) <- rank_df$gene
ranked <- sort(ranked, decreasing = TRUE)    # most up first

head(ranked, 8); tail(ranked, 8)
gsea_h <- GSEA(geneList = ranked, TERM2GENE = h_sets,
               minGSSize = 15, maxGSSize = 500,
               pvalueCutoff = 0.05, pAdjustMethod = "BH",
               eps = 0, seed = TRUE)

gsea_df <- as.data.frame(gsea_h)
head(gsea_df[, c("ID", "setSize", "NES", "p.adjust")], 10)

NES (normalised enrichment score) is the number to report. Positive = the pathway’s genes are mostly up in COVID-19; negative = mostly down.

# KEGG needs Entrez IDs, so translate the names of the ranked list.
map_all <- bitr(names(ranked), fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)
ranked_entrez <- ranked[map_all$SYMBOL]
names(ranked_entrez) <- map_all$ENTREZID
ranked_entrez <- sort(ranked_entrez[!duplicated(names(ranked_entrez))], decreasing = TRUE)

gsea_kegg <- gseKEGG(geneList = ranked_entrez, organism = "hsa",
                     minGSSize = 15, pvalueCutoff = 0.05, eps = 0, seed = TRUE)
gsea_kegg <- setReadable(gsea_kegg, OrgDb = org.Hs.eg.db, keyType = "ENTREZID")

dotplot(gsea_kegg, x = "NES", showCategory = 10, split = ".sign") +
  facet_grid(. ~ .sign) +
  ggtitle("GSEA, KEGG — COVID-19 vs healthy")

3.2.1 The enrichment plot

# Top: running score (a peak on the left = pathway genes are at the "up" end).
# Middle: where the pathway's genes sit in the ranked list.
gseaplot2(gsea_h, geneSetID = "HALLMARK_INTERFERON_ALPHA_RESPONSE",
          title = "HALLMARK_INTERFERON_ALPHA_RESPONSE", pvalue_table = TRUE)

# Three pathways together
sel <- head(order(gsea_df$p.adjust), 3)
gseaplot2(gsea_h, geneSetID = sel, title = "Three most significant Hallmark pathways")

# The "leading edge": the genes that actually drive the enrichment.
le <- gsea_df$core_enrichment[gsea_df$ID == "HALLMARK_INTERFERON_ALPHA_RESPONSE"]
strsplit(le, "/")[[1]][1:20]

ORA vs GSEA. ORA asks “is my list enriched?”; GSEA asks “is this pathway shifted in my ranking?”. They usually agree on strong signals (interferon) and differ on weak, coordinated ones, where GSEA is more sensitive.

write.csv(as.data.frame(ora_h_up), file.path(res_dir, "ORA_hallmark_up.csv"), row.names = FALSE)
write.csv(as.data.frame(ora_kegg_up), file.path(res_dir, "ora_kegg_up_up.csv"),     row.names = FALSE)
write.csv(gsea_df,                 file.path(res_dir, "GSEA_hallmark.csv"),   row.names = FALSE)

Further reading

  • Love MI, Huber W, Anders S (2014) Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology 15:550.
  • McClain MT et al. (2021) Nat Commun 12:1079 — the GSE161731 paper.
  • Extra: Subramanian A et al. (2005) Gene set enrichment analysis. PNAS 102:15545; Wu T et al. (2021) clusterProfiler 4.0. The Innovation 2:100141.