Försöksdjursvetenskap (teoretisk) och avancerad databearbetning
Lab 7 — Sample exploration and differential expression: PCA and DESeq2
Birong Zhang
21 October 2026
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:
- 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.
- DESeq2 — Identify differentially expressed genes. We compare COVID-19 with healthy for every gene, correct for the fact that we test ~15,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 …) must 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
# 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 <- readRDS(file.path(data_dir, "gse161731_counts.rds")) # raw counts
meta <- readRDS(file.path(data_dir, "gse161731_meta.rds")) # clean patient table
ann <- readRDS(file.path(data_dir, "gse161731_gene_annotation.rds")) # gene annotation
dim(counts)## [1] 35447 177##
## healthy Bacterial CoV_other COVID_19 Influenza
## 19 24 40 77 17## [1] TRUE1 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.4745891.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## [1] 17980 1771.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.31810par(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")
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)
ggplot(pca_df, aes(x = PC1, y = PC2, colour = cohort)) +
geom_point(size = 2.5, alpha = 0.85) +
scale_colour_brewer(palette = "Set2") +
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).
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) +
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## 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.13545180Interpret. In blood, the top PCs are usually a mix of the infection response (interferon genes) and which cells are in the tube (neutrophil vs lymphocyte genes). Look at the gene names above: which of the two do you see?
1.4 Exercise 1
1.4.1 Question
- Redo the PCA with the 2000 most variable genes. What % does PC1 explain now?
- Draw PC1 as a box plot per cohort. Which cohorts separate?
- 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## [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.
- Redo the PCA with the 1000 most variable genes and report the % explained by PC1 and PC2.
- Draw PC2 as a box plot per cohort, using the colour
palette
"Set1". Which cohorts separate on PC2? - Is S100A8 among the 20 genes with the largest PC1 loadings? In one or two sentences: is PC1 mainly about interferon (virus) or about neutrophils (cell composition / bacteria)?
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_f, not CPM.
# Keep the two groups we want to compare.
meta_de <- meta[meta$cohort %in% c("healthy", "COVID_19"), ]
meta_de$cohort <- droplevels(meta_de$cohort) # forget the 3 unused groups
counts_de <- counts_f[, rownames(meta_de)] # the same samples, same order
table(meta_de$cohort)##
## healthy COVID_19
## 19 77## [1] TRUEThe design formula tells DESeq2 what explains the
differences between samples. Here it is simply ~ cohort:
“gene expression depends on which group the patient is in”. (Adding age
and sex, ~ age + gender + cohort, is shown in the Extra
section.)
dds <- DESeqDataSetFromMatrix(
countData = counts_de, # raw counts — never CPM
colData = meta_de, # the patient table, same order as the columns
design = ~ cohort # the comparison
)
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 <- results(dds, contrast = c("cohort", "COVID_19", "healthy"), alpha = 0.05)
summary(res)##
## out of 17980 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up) : 1334, 7.4%
## LFC < 0 (down) : 2155, 12%
## outliers [1] : 0, 0%
## low counts [2] : 3486, 19%
## (mean count < 19)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?resultsWhat 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 just tested about 15,000 genes. With a cut-off of p < 0.05, 1 in 20 genes with no real change will still look significant by chance — that could be ~750 false discoveries!
The adjusted p-value (padj) corrects
for this. If you keep all genes with padj < 0.05, you
expect only about 5 % of your list to be false
discoveries. Always use padj, never the raw
pvalue, to call genes significant.
## [1] 5427## [1] 3489Why 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.
## [1] 3486## 19.38776%
## 18.551872.3 Reading and visualising the result
2.3.1 The results table
res_df <- as.data.frame(res)
res_df$gene <- rownames(res_df) # keep the gene name as a column
res_df <- res_df[order(res_df$padj), ] # most significant genes first
head(res_df[, 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_df[res_df$gene %in% check, c("gene", "log2FoldChange", "padj")]2.3.2 Volcano plot
Each dot is a gene. Left–right: how much it changes. Up: how significant it is. The interesting genes are in the top corners.
vol <- res_df[!is.na(res_df$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
## 324 13927 243# the 20 most significant genes get a label
top_lab <- rbind(
head(vol[vol$log2FoldChange > 0, ][order(vol[vol$log2FoldChange > 0, ]$padj), ], 10),
head(vol[vol$log2FoldChange < 0, ][order(vol[vol$log2FoldChange < 0, ]$padj), ], 10)
)
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.
best <- res_df$gene[1] # the most significant gene
plotCounts(dds, gene = best, intgroup = "cohort", main = best)
2.3.4 Heatmap of the top genes
top50 <- head(res_df$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
- How many genes have
padj < 0.05, and how many of those also change more than 2-fold (|log2FC| > 1)? - Redraw the volcano plot, labelling only the top 10
genes and colouring up-regulated genes
"red"and down-regulated genes"blue". - Make a
plotCounts()figure for IFI27. Does the direction fit a viral infection?
2.4.2 Solution
## [1] 3489## [1] 567# 2 Changed: head(vol, 10) and the two colours.
top_lab10 <- rbind(
head(vol[vol$log2FoldChange > 0, ][order(vol[vol$log2FoldChange > 0, ]$padj), ], 5),
head(vol[vol$log2FoldChange < 0, ][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.
- How many genes have
padj < 0.01, and how many of those also change more than 4-fold (|log2FC| > 2)? - Redraw the volcano plot, labelling only the top 5
genes and colouring up-regulated genes
"darkorange"and down-regulated genes"purple". - Make a
plotCounts()figure for CD3D (a T-cell gene). Is it higher or lower in COVID-19? In one sentence, link this to lymphopenia (fewer lymphocytes in the blood), which is common in severe COVID-19.
3 Save the results
write.csv(res_df, file.path(res_dir, "DE_COVID19_vs_healthy.csv"), row.names = FALSE)
saveRDS(dds, file.path(res_dir, "dds_covid_vs_healthy.rds"))
sessionInfo()## R version 4.5.0 (2025-04-11)
## Platform: aarch64-apple-darwin20
## Running under: macOS 26.7
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: Europe/Stockholm
## tzcode source: internal
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] RColorBrewer_1.1-3 pheatmap_1.0.13 ggrepel_0.9.8
## [4] dplyr_1.2.1 ggplot2_4.0.3 edgeR_4.6.3
## [7] limma_3.64.3 DESeq2_1.48.2 SummarizedExperiment_1.40.0
## [10] Biobase_2.70.0 MatrixGenerics_1.22.0 matrixStats_1.5.0
## [13] GenomicRanges_1.62.1 Seqinfo_1.0.0 IRanges_2.44.0
## [16] S4Vectors_0.48.1 BiocGenerics_0.56.0 generics_0.1.4
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.6 xfun_0.60 bslib_0.12.0 lattice_0.23-1
## [5] vctrs_0.7.3 tools_4.5.0 parallel_4.5.0 tibble_3.3.1
## [9] pkgconfig_2.0.3 Matrix_1.7-6 S7_0.2.2 lifecycle_1.0.5
## [13] GenomeInfoDbData_1.2.14 compiler_4.5.0 farver_2.1.2 statmod_1.5.2
## [17] codetools_0.2-20 GenomeInfoDb_1.44.3 htmltools_0.5.9 sass_0.4.10
## [21] yaml_2.3.12 pillar_1.11.1 jquerylib_0.1.4 BiocParallel_1.44.0
## [25] DelayedArray_0.36.1 cachem_1.1.0 abind_1.4-8 tidyselect_1.2.1
## [29] locfit_1.5-9.12 digest_0.6.39 labeling_0.4.3 fastmap_1.2.0
## [33] grid_4.5.0 cli_3.6.6 SparseArray_1.10.10 magrittr_2.0.5
## [37] S4Arrays_1.10.1 withr_3.0.3 scales_1.4.0 UCSC.utils_1.4.0
## [41] rmarkdown_2.32 XVector_0.50.0 httr_1.4.9 otel_0.2.0
## [45] evaluate_1.0.5 knitr_1.51 viridisLite_0.4.3 rlang_1.3.0
## [49] Rcpp_1.1.2 glue_1.8.1 rstudioapi_0.19.0 jsonlite_2.0.0
## [53] R6_2.6.1Take-home messages
- Raw counts depend on sequencing depth; CPM makes samples comparable for plots.
- 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).
- DESeq2 needs raw counts and a clear comparison: which group is the reference.
- Testing thousands of genes creates false positives by chance — use padj.
- 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.
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 sets3.1 Adjusting for age and sex
In Lab 6 we saw that COVID-19 patients and healthy controls may differ in age. We can ask DESeq2 to take age and sex into account.
dds_adj <- DESeqDataSetFromMatrix(countData = counts_de, colData = meta_de,
design = ~ age + gender + cohort) # variable of interest last
dds_adj <- DESeq(dds_adj)
res_adj <- results(dds_adj, contrast = c("cohort", "COVID_19", "healthy"), alpha = 0.05)
sum(res$padj < 0.05, na.rm = TRUE) # simple model
sum(res_adj$padj < 0.05, na.rm = TRUE) # adjusted model
a <- rownames(res)[which(res$padj < 0.05)]
b <- rownames(res_adj)[which(res_adj$padj < 0.05)]
length(intersect(a, b)) # significant in both models3.2 ORA — over-representation analysis
ORA needs three things:
- the gene list (significant genes, up and down separately);
- the gene sets (pathways);
- 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_df$gene[!is.na(res_df$pvalue)] # every gene that was tested
genes_up <- res_df$gene[!is.na(res_df$padj) & res_df$padj < 0.05 & res_df$log2FoldChange > 1]
genes_down <- res_df$gene[!is.na(res_df$padj) & res_df$padj < 0.05 & res_df$log2FoldChange < -1]
length(background); length(genes_up); length(genes_down)3.2.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")3.2.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 <- enrichKEGG(gene = map_up$ENTREZID, universe = map_bg$ENTREZID,
organism = "hsa", pAdjustMethod = "BH", pvalueCutoff = 0.05)
ora_kegg <- setReadable(ora_kegg, OrgDb = org.Hs.eg.db, keyType = "ENTREZID") # IDs back to symbols
head(as.data.frame(ora_kegg)[, c("Description", "p.adjust", "Count")], 10)
barplot(ora_kegg, showCategory = 10, title = "ORA, KEGG — up in COVID-19")
as.data.frame(ora_kegg)$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.3 Try it yourself
3.2.3.1 Question
- Re-run the Hallmark ORA with a looser list (padj < 0.05, any fold change up). How many pathways are significant now?
- Re-run it with the wrong background (all genes in
logcpm). Do the p-values get smaller or larger? - How many genes are in
HALLMARK_INTERFERON_ALPHA_RESPONSE, and how many of them are in your up-regulated list?
3.2.3.2 Solution
# 1
genes_up_loose <- res_df$gene[!is.na(res_df$padj) & res_df$padj < 0.05 & res_df$log2FoldChange > 0]
ora_loose <- enricher(genes_up_loose, universe = background, TERM2GENE = h_sets)
nrow(as.data.frame(ora_loose)); nrow(as.data.frame(ora_h_up))
# 2
ora_wrongbg <- enricher(genes_up, universe = rownames(logcpm), TERM2GENE = h_sets)
head(as.data.frame(ora_wrongbg)[, c("ID", "p.adjust")], 5)
head(as.data.frame(ora_h_up)[, c("ID", "p.adjust")], 5)
# 3
ifn_set <- h_sets$gene_symbol[h_sets$gs_name == "HALLMARK_INTERFERON_ALPHA_RESPONSE"]
length(ifn_set)
length(intersect(ifn_set, genes_up))3.3 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.3.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.
3.3.2 Try it yourself
3.3.2.1 Question
- Which Hallmark pathway has the most negative NES, and what might that mean for the blood of COVID-19 patients?
- Draw the enrichment plot for
HALLMARK_INFLAMMATORY_RESPONSEand report its NES and adjusted p-value. - Rank by
log2FoldChangeinstead ofstat. Do the top three pathways change?
3.3.2.2 Solution
# 1
gsea_df[which.min(gsea_df$NES), c("ID", "NES", "p.adjust")]
# 2
gseaplot2(gsea_h, geneSetID = "HALLMARK_INFLAMMATORY_RESPONSE",
title = "HALLMARK_INFLAMMATORY_RESPONSE", pvalue_table = TRUE)
gsea_df[gsea_df$ID == "HALLMARK_INFLAMMATORY_RESPONSE", c("NES", "p.adjust")]
# 3
rank_lfc <- sort(setNames(rank_df$log2FoldChange, rank_df$gene), decreasing = TRUE)
gsea_lfc <- GSEA(rank_lfc, TERM2GENE = h_sets, minGSSize = 15,
pvalueCutoff = 0.05, eps = 0, seed = TRUE)
head(as.data.frame(gsea_lfc)[, c("ID", "NES", "p.adjust")], 3)
head(gsea_df[, c("ID", "NES", "p.adjust")], 3)
# `stat` is usually better: a big fold change measured on 4 reads should not
# outrank a moderate one measured on 4000.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.