How today is organised

Four hours, five blocks, the same shape as Lab 4. Each block ends with a small exercise; the solution is in the next tab. The five exercises are the practice version of the five assignment tasks in section.

	Block	Time	Ends with
1	PCA and sample–sample distances: quality control before modelling	45 min	Exercise 1
2	DESeq2 step by step: design, size factors, dispersion, testing	80 min	Exercise 2
	break	15 min
3	Reading and visualising a DE result; covariates and thresholds	30 min	Exercise 3
4	Over-representation analysis (ORA): KEGG and MSigDB Hallmark	35 min	Exercise 4
5	GSEA: Hallmark and KEGG, dot plots and enrichment plots	35 min	Exercise 5

0.1 Introduction

This lab builds on the count matrix and metadata prepared in Lab 4 and moves from data handling to statistical inference. We begin with principal component analysis and sample-to-sample distances as exploratory quality control, then work through differential expression with DESeq2 and the interpretation of the results table. The differentially expressed genes are visualised with volcano, per-gene count and heatmap plots. The lab closes with pathway-level interpretation using KEGG and the MSigDB Hallmark collection, covering both over-representation analysis and gene set enrichment analysis, with attention to the choice of background set and ranking metric.

0.2 Learning objectives

By the end of Lab 5 you should be able to:

LO1 — Explore structure. Run a PCA on log2 CPM values, decide how many components matter, and identify which clinical variable drives each component.
LO2 — Model. Build a DESeq2 object with the correct design formula.
LO3 — Interpret. Read a results table correctly (log2FC, p, padj, NA values), and produce the standard volcano / heatmap figures.
LO4 — ORA. Run over-representation analysis on a significant gene list against KEGG and Hallmark, with the correct background, and read the dot plot.
LO5 — GSEA. Rank all genes, run GSEA, and interpret the normalised enrichment score and the running-score plot.

0.3 Assignment

Assignment 1 — PCA (LO1)

Produce a PCA of the samples coloured by cohort, and report the percentage of variance explained by PC1 and PC2.
Test whether PC1 is associated with cohort, and whether it is associated with age.
In two sentences: is the main axis of variation biological, technical, or driven by blood cell composition — and how did you decide?

Assignment 2 — DESeq2 model (LO2)

Build a DESeqDataSet for COVID-19 vs healthy with healthy as the reference level, and say in one sentence why the reference level matters.
Refit the model including age and sex as covariates and report how many genes change status (significant / not significant) compared with the simple model.
Compare the significant gene sets between the two models.

Assignment 3 — Reading the result (LO3)

Report the number of genes with padj < 0.05 and |log2FC| > 1, and list the top ten by adjusted p-value.
Visualize and interpret key genes using plotCounts() function.
Explain in two sentences why some genes have padj = NA, and check where three interferon genes of your choice sit in the ranking.

Assignment 4 — ORA (LO4)

Run ORA on the up-regulated genes against KEGG and against MSigDB Hallmark, using all tested genes as the background. Show both dot plots (top 10).
Repeat for the down-regulated genes and describe how the two answers differ.
In two sentences: name one pathway you expected to see and one you did not, and give a possible biological explanation for the unexpected one.

Assignment 5 — GSEA (LO5)

Build the ranked gene list, run GSEA against Hallmark, and show the dot plot of the top 10 pathways by adjusted p-value.
Show an enrichment plot for HALLMARK_INTERFERON_ALPHA_RESPONSE and a second plot combining three pathways.
In two or three sentences: what does a positive versus a negative NES mean here, and why can GSEA find a pathway that ORA misses?

0.4 Packages

BiocManager::install(c("DESeq2", "apeglm", "clusterProfiler", "enrichplot",
                       "org.Hs.eg.db", "DOSE"))
install.packages(c("ggplot2", "dplyr", "tidyr", "ggrepel", "pheatmap",
                   "RColorBrewer", "msigdbr"))

library(DESeq2)            # differential expression
library(apeglm)            # log fold change shrinkage
library(clusterProfiler)   # ORA and GSEA
library(enrichplot)        # dotplot(), gseaplot2()
library(org.Hs.eg.db)      # gene ID conversion
library(msigdbr)           # MSigDB gene sets (Hallmark)
library(ggplot2)
library(dplyr)
library(tidyr)
library(ggrepel)           # non-overlapping labels on the volcano plot
library(pheatmap)
library(RColorBrewer)

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

0.5 Loading the objects from Lab 4

counts <- readRDS(file.path(data_dir, "gse161731_counts_filtered.rds"))  # raw counts, filtered
meta   <- readRDS(file.path(data_dir, "gse161731_meta.rds"))             # clean sample table
logcpm <- readRDS(file.path(data_dir, "gse161731_logcpm.rds"))           # log2 CPM
ann    <- readRDS(file.path(data_dir, "gse161731_gene_annotation.rds"))  # symbol / ensembl

dim(counts)

## [1] 17549   137

table(meta$cohort)

## 
## Bacterial CoV_other  COVID_19   healthy Influenza 
##        24        35        45        16        17

all(colnames(counts) == rownames(meta))   # the alignment check, one more time

## [1] TRUE

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

Principal component analysis rotates the samples into new axes ordered by how much variance they explain. PC1 is the single direction along which the samples differ most. It answers the first question of every RNA-seq project: is the biggest difference in my data the difference I care about?

1.1 Doing it step by step

# 1. Use the most variable genes. Genes that barely vary add noise, not signal.
gene_var  <- apply(logcpm, 1, var)                       # variance of each row (gene)
top_genes <- names(sort(gene_var, decreasing = TRUE))[1:500]
mat       <- logcpm[top_genes, ]
dim(mat)

## [1] 500 137

# 2. prcomp() expects SAMPLES in rows, so transpose the matrix.
#    center = TRUE  : subtract the mean of each gene (always do this)
#    scale. = FALSE : keep the genes' relative variances (usual choice for log CPM)
pca <- prcomp(t(mat), center = TRUE, scale. = FALSE)

# 3. Percentage of variance explained by each component.
var_pct <- round(100 * pca$sdev^2 / sum(pca$sdev^2), 1)
var_pct[1:6]

## [1] 16.0  9.3  6.5  5.2  3.6  3.0

# A scree plot: how many components carry real structure?
barplot(var_pct[1:10], names.arg = paste0("PC", 1:10),
        ylab = "% variance explained", main = "Scree plot")

# pca$x holds the sample coordinates. Put the first two PCs in a data.frame with
# the metadata so ggplot can colour by any clinical variable.
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")

This pattern is actually expected, not an error. Bacterial infections trigger a much stronger and more distinct innate immune response (neutrophil activation, myeloid gene upregulation) than viral infections, so bacterial samples separate clearly along the top variable genes. In contrast, COVID-19, other coronaviruses, and influenza all converge on a similar interferon-stimulated gene (ISG) response, so they overlap heavily with each other. Many “healthy” samples in GSE161731 may also be mild or early/convalescent-phase individuals whose transcriptional signature isn’t strongly different from mild viral infection, which is why they cluster together too. Note that PC1 and PC2 here only explain ~24% of total variance combined, so finer distinctions between viral subtypes and healthy controls are likely captured in higher-order PCs rather than this 2D projection — you may need supervised analysis (e.g., differential expression, ISG-focused gene panels, or classification models) to separate them more clearly.

1.2 Which variable drives which component?

Colour the same plot by other variables — and test it, do not guess.

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

# Continuous variable vs a PC -> correlation
cor.test(pca_df$PC1, pca_df$age, method = "spearman")

## 
##  Spearman's rank correlation rho
## 
## data:  pca_df$PC1 and pca_df$age
## S = 584344, p-value = 1.259e-05
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.3635814

# Categorical variable vs a PC -> Kruskal-Wallis (more than two groups)
kruskal.test(PC1 ~ cohort, data = pca_df)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  PC1 by cohort
## Kruskal-Wallis chi-squared = 68.275, df = 4, p-value = 5.248e-14

kruskal.test(PC2 ~ cohort, data = pca_df)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  PC2 by cohort
## Kruskal-Wallis chi-squared = 1.7648, df = 4, p-value = 0.7789

# Which genes push the samples along PC1? The largest loadings.
load1 <- sort(pca$rotation[, 1])
head(load1, 10)     # strongest negative

##       CD177        MMP8      SLC51A     ADAMTS2     PCOLCE2       TMT1B        OLAH       KLF14 
## -0.13322090 -0.11964174 -0.10967879 -0.10555907 -0.10446258 -0.10253619 -0.09778483 -0.09504801 
##      NECAB1      FAM20A 
## -0.09312594 -0.09258384

tail(load1, 10)     # strongest positive

##     CACNG6    TRBV6-3      LRRN3      TRDV2     PTGDR2      OLIG2        NOG    SIGLEC8     ALOX15 
## 0.08345015 0.08353969 0.08437038 0.08633682 0.08880176 0.08917472 0.09510913 0.11383173 0.12398797 
##     PRSS33 
## 0.12773563

Interpret. In whole blood the top PCs are usually a mix of the infection response (interferon genes) and the cell composition of the tube (neutrophil vs lymphocyte genes). Look at the loadings above and decide which one you are seeing.

1.3 Sample–sample distances as a complement

# dist() computes the Euclidean distance between the ROWS, so transpose again.
d   <- dist(t(mat))
dmat <- as.matrix(d)

pheatmap(dmat,
          cluster_rows = FALSE,
         cluster_cols = FALSE,
         clustering_distance_rows = d,     # cluster using the same distances
         clustering_distance_cols = d,
         annotation_col = data.frame(cohort = meta$cohort, row.names = rownames(meta)),
         show_rownames = FALSE, 
         show_colnames = FALSE,
         color = colorRampPalette(rev(brewer.pal(9, "RdBu")))(100),
         main = "Sample-to-sample distances")

1.4 Exercise 1

1.4.1 Question

Redo the PCA using the 2000 most variable genes. Does PC1 explain more or less variance, and does the picture change?
Plot PC1 as a box plot per cohort. Which cohorts separate?
Is IFI27 among the top 20 PC1 loadings?

1.4.2 Solution

# 1
top2k  <- names(sort(gene_var, decreasing = TRUE))[1:2000]
pca2   <- prcomp(t(logcpm[top2k, ]), center = TRUE, scale. = FALSE)
round(100 * pca2$sdev[1:3]^2 / sum(pca2$sdev^2), 1)

## [1] 16.9  8.6  4.3

# 2
ggplot(pca_df, aes(cohort, PC1, fill = cohort)) +
  geom_boxplot(alpha = .7, outlier.shape = NA) + geom_jitter(width = .15, size = .8) +
  theme(legend.position = "none", axis.text.x = element_text(angle = 45, hjust = 1))

# 3  Rank genes by the absolute size of their PC1 loading.
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"   "SLC51A"    "ADAMTS2"   "PCOLCE2"  
##  [9] "TMT1B"     "OLAH"      "NOG"       "KLF14"     "NECAB1"    "FAM20A"    "OLIG2"     "TNFAIP8L3"
## [17] "PTGDR2"    "CYP19A1"   "AP3B2"     "TRDV2"

2 Block 2 — DESeq2, step by step

2.1 The comparison, and the design formula

# Keep the two cohorts we want to compare.
meta <- meta %>% arrange(desc(meta$cohort))
keep_s   <- meta$cohort %in% c("healthy", "COVID_19")
meta_de  <- meta[keep_s, ]
counts_de <- counts[, rownames(meta_de)]     # DESeq2 wants a matrix of integers

# droplevels() removes the three cohorts we are not using from the factor.
meta_de$cohort <- droplevels(meta_de$cohort)

# relevel(): make "healthy" the reference. All log2 fold changes will then be
# "COVID-19 relative to healthy" — positive = higher in COVID-19.
meta_de$cohort <- relevel(meta_de$cohort, ref = "healthy")

table(meta_de$cohort)

## 
##  healthy COVID_19 
##       16       45

dim(counts_de)

## [1] 17549    61

all(rownames(meta_de) == colnames(counts_de))

## [1] TRUE

A design formula lists the variables the model must account for. The variable of interest goes last.

~ cohort — the simple comparison.
~ age + gender + cohort — the same comparison, adjusted for age and sex. Use this when a covariate differs between groups (as age does here; see Lab 4).

dds <- DESeqDataSetFromMatrix(
  countData = counts_de,     # raw counts — never CPM or TPM
  colData   = meta_de,       # sample table, rows in the same order as the columns
  design    = ~ cohort       # we start simple; we add covariates later
)
dds

## class: DESeqDataSet 
## dim: 17549 61 
## metadata(1): version
## assays(1): counts
## rownames(17549): ACTB MT-CO1 ... EIF2AP1 RN7SKP134
## rowData names(0):
## colnames(61): SDU09.02S0000103 SDU09.02S0000113 ... SDU18.02S0011631 SDU18.02S0011633
## colData names(8): subject_id age ... hospitalized batch

# A light extra filter: keep genes with at least 10 reads in total.
# (Lab 4 already filtered; this only removes leftovers.)
dds <- dds[rowSums(counts(dds)) >= 10, ]
nrow(dds)

## [1] 17549

# Everything above can be run in one line — `dds <- DESeq(dds)` — which does size factors, dispersion and the test together
dds <- DESeq(dds) # 1. dds <- estimateSizeFactors(dds); 2. dds <- estimateDispersions(dds); 3. dds <- nbinomWaldTest(dds) 
resultsNames(dds)

## [1] "Intercept"                  "cohort_COVID_19_vs_healthy"

# results(): pull out the comparison. contrast = c(variable, numerator, denominator)
res <- results(dds, contrast = c("cohort", "COVID_19", "healthy"), alpha = 0.05)
summary(res)

## 
## out of 17549 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 489, 2.8%
## LFC < 0 (down)     : 443, 2.5%
## outliers [1]       : 0, 0%
## low counts [2]     : 3403, 19%
## (mean count < 24)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

head(as.data.frame(res))

# What the columns mean:
# baseMean        mean of the normalised counts across all samples
# log2FoldChange  effect size; +1 = twice as high in COVID-19
# lfcSE           standard error of that estimate
# stat            log2FoldChange / lfcSE, the Wald statistic (used for GSEA later)
# pvalue          raw p-value
# padj            Benjamini-Hochberg adjusted p-value (FDR)
colnames(as.data.frame(res))

## [1] "baseMean"       "log2FoldChange" "lfcSE"          "stat"           "pvalue"        
## [6] "padj"

2.2 Multiple testing

# We tested ~15,000 genes. At p < 0.05 we would expect ~750 false positives by chance.
sum(res$pvalue < 0.05, na.rm = TRUE) #4098

## [1] 3130

sum(res$padj   < 0.05, na.rm = TRUE) #1766

## [1] 932

hist(res$pvalue, breaks = 50, col = "grey", main = "p-value histogram", xlab = "p")

Read the histogram. A flat background with a spike near zero is healthy. A hill in the middle or a spike near 1 means the model is misspecified.

2.3 Exercise 2

2.3.1 Question

How many genes are significant at padj < 0.05, and how many of those also have |log2FC| > 1?
Refit the model as ~ age + gender + cohort and report the number of significant genes. Did it go up or down?
How much do the two models’ gene lists overlap?

2.3.2 Solution

# 1
sum(res$padj < 0.05, na.rm = TRUE)

## [1] 932

sum(res$padj < 0.05 & abs(res$log2FoldChange) > 1, na.rm = TRUE)

## [1] 311

# 2
# dds_adj <- dds
# design(dds_adj) <- ~ age + gender + cohort   # replace the design…
dds_adj <- DESeqDataSetFromMatrix(
  countData = counts_de,   
  colData   = meta_de,       
  design    = ~ cohort + age + gender    # add covariates
)
dds_adj

## class: DESeqDataSet 
## dim: 17549 61 
## metadata(1): version
## assays(1): counts
## rownames(17549): ACTB MT-CO1 ... EIF2AP1 RN7SKP134
## rowData names(0):
## colnames(61): SDU09.02S0000103 SDU09.02S0000113 ... SDU18.02S0011631 SDU18.02S0011633
## colData names(8): subject_id age ... hospitalized batch

dds_adj <- DESeq(dds_adj)                    # …and refit everything
res_adj <- results(dds_adj, contrast = c("cohort", "COVID_19", "healthy"), alpha = 0.05)
sum(res_adj$padj < 0.05, na.rm = TRUE)

## [1] 549

# 3
a <- rownames(res)[which(res$padj < 0.05)]
b <- rownames(res_adj)[which(res_adj$padj < 0.05)]
length(intersect(a, b)); length(setdiff(a, b)); length(setdiff(b, a))

## [1] 314

## [1] 618

## [1] 235

3 Block 3 — Reading and visualising the result

3.1 The table

res_df <- as.data.frame(res)
head(res_df)

res_df$gene <- rownames(res_df)          # keep the symbol as a column
res_df <- res_df[order(res_df$padj), ]   # best genes first

head(res_df[, c("gene", "baseMean", "log2FoldChange", "pvalue", "padj")], 15)

# Why are some p/padj values NA? 
# all-zero counts, or an extreme Cook's-distance outlier
sum(is.na(res_df$pvalue))

## [1] 0

# independent filtering: low-count genes are removed before multiple-testing correction, which increases power for the genes that remain
sum(is.na(res_df$padj) & !is.na(res_df$pvalue))

## [1] 3403

metadata(res)$filterThreshold      # the baseMean cut DESeq2 chose automatically

## 19.38776% 
##  24.16068

# The positive control: where do known interferon genes sit?
check <- c("IFI27", "IFI44L", "ISG15", "RSAD2", "OAS1", "MX1", "SIGLEC1", "CD3D", "IL7R", "S100A8")
res_df[res_df$gene %in% check, c("gene", "log2FoldChange", "padj")]

3.2 Volcano plot

vol <- res_df[!is.na(res_df$padj), ]                    # drop genes with no padj

# A status column, used only for colouring.
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 
##               99            13835              212

top_lab <- head(vol[order(vol$padj), ], 20)             # label the 15 best genes

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") +
  geom_hline(yintercept = -log10(0.05), linetype = 2, colour = "grey40") +
  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")

3.3 The top genes, one by one and as a heatmap

# plotCounts() draws the normalised counts of one gene, split by group.
best <- rownames(res_df)[1]
plotCounts(dds, gene = best, intgroup = "cohort", main = best) # DESeq2 normalized count = raw count/size factor

top50 <- head(res_df$gene[!is.na(res_df$padj)], 50)   # 50 best genes
mat50 <- logcpm[top50, colnames(counts_de)]           # their log2 CPM, these samples

pheatmap(mat50,
         scale          = "row",
         annotation_col = data.frame(cohort = meta_de$cohort, row.names = rownames(meta_de)),
         show_colnames  = FALSE,
         fontsize_row   = 7,
         cluster_rows      = TRUE,    
         cluster_cols      = FALSE,      
         color          = colorRampPalette(rev(brewer.pal(9, "RdBu")))(100),
         main           = "Top 50 differentially expressed genes")

3.4 Exercise 3

3.4.1 Question

List the ten most up-regulated and the ten most down-regulated genes with padj < 0.05, sorted by log2FC.
Make a plotCounts() figure for IFI27 and for CD3D. Do the directions match the biology of viral infection and lymphopenia?
How many genes have padj = NA, and what is the baseMean threshold DESeq2 used?

3.4.2 Solution

# 1
sig <- res_df[!is.na(res_df$padj) & res_df$padj < 0.05, ]
head(sig[order(-sig$log2FoldChange), c("gene", "log2FoldChange", "padj")], 10)

head(sig[order( sig$log2FoldChange), c("gene", "log2FoldChange", "padj")], 10)

# 2
plotCounts(dds, gene = "IFI27", intgroup = "cohort", main = "IFI27")

plotCounts(dds, gene = "CD3D",  intgroup = "cohort", main = "CD3D")

# 3
sum(is.na(res$padj))

## [1] 3403

metadata(res)$filterThreshold

## 19.38776% 
##  24.16068

4 Block 4 — Over-representation analysis (ORA)

4.1 The idea, and the one thing people get wrong

ORA asks: among my 800 significant genes, are there more interferon genes than you would expect by chance? It is a Fisher test on a 2×2 table, so it needs three inputs:

the gene list (significant genes, usually split into up and down),
the gene sets (KEGG pathways, Hallmark collections),
the background — and this must be every gene that was actually tested, not “all human genes”. Using the whole genome as background inflates the significance of everything, because our matrix already excludes genes not expressed in blood.

# Background: every gene with a p-value, i.e. every gene the model actually tested.
background <- res_df$gene[!is.na(res_df$pvalue)]
length(background)

## [1] 17549

# The two gene lists.
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(genes_up); length(genes_down)

## [1] 212

## [1] 99

4.2 ORA against MSigDB Hallmark

The 50 Hallmark sets are curated, non-redundant summaries of well-defined biological states — the best first choice for a clinical transcriptomics study.

# msigdbr() downloads the sets as a data.frame, one row per gene-set membership.
# Newer versions use `collection = "H"`; older versions use `category = "H"`.
hallmark <- msigdbr(species = "Homo sapiens", collection = "H")
# hallmark <- msigdbr(species = "Homo sapiens", category = "H")   # older msigdbr

# clusterProfiler wants a two-column table: set name, then gene.
h_sets <- as.data.frame(hallmark[, c("gs_name", "gene_symbol")])
head(h_sets)

length(unique(h_sets$gs_name))    # 50 Hallmark sets

## [1] 50

ora_h_up <- enricher(
  gene          = genes_up,       # our up-regulated genes
  universe      = background,     # the background
  TERM2GENE     = h_sets,         # the gene 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 genes in the overlap; colour = adjusted p-value;
# x position = GeneRatio, the fraction of our list that falls in that set.
dotplot(ora_h_up, showCategory = 10,
        title = "ORA, Hallmark — up in COVID-19") +
  theme(axis.text.y = element_text(size = 9))

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)

# dotplot(ora_h_down, showCategory = 10,
#         title = "ORA, Hallmark — down in COVID-19") +
#   theme(axis.text.y = element_text(size = 9))

4.3 ORA against KEGG

KEGG pathways are indexed by Entrez IDs, so we must convert our symbols first. Some symbols will not convert — report how many.

# bitr() = "biological id translator". It returns a two-column data.frame.
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)

nrow(map_up); length(genes_up)     # how many of our genes could be converted

## [1] 212

## [1] 212

head(map_up)

ora_kegg <- enrichKEGG(
  gene          = map_up$ENTREZID,
  universe      = map_bg$ENTREZID,
  organism      = "hsa",           # human
  pAdjustMethod = "BH",
  pvalueCutoff  = 0.05
)

# setReadable() turns the Entrez IDs in the result back into gene symbols.
ora_kegg <- setReadable(ora_kegg, OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
head(as.data.frame(ora_kegg)[, c("ID", "Description", "p.adjust", "Count")], 10)

barplot(ora_kegg, showCategory = 10, title = "ORA, KEGG — up in COVID-19") +
  theme(axis.text.y = element_text(size = 9))

# Which of our genes drove the top pathway? Always look.
as.data.frame(ora_kegg)$geneID[1]

## [1] "CCNB2/CDC20/BUB1/PKMYT1/CCNB1/CDT1/CDC25A/CDCA5/PLK1/CCNA1/E2F2/CCNA2/ORC1/CDK1/CDC45/CDC6/TICRR/ESPL1/TRIP13"

Pitfalls to state out loud. (1) The result depends entirely on your significance cut-off — change padj or the fold-change threshold and the pathways change. (2) A gene absent from every database cannot be found, however important it is. (3) KEGG contains disease pathways (“Coronavirus disease – COVID-19”) that are partly built from studies like this one; finding them is reassuring, not a discovery.

4.4 Exercise 4

4.4.1 Question

Re-run the Hallmark ORA using a looser gene list (padj < 0.05, no fold-change cut-off). How many pathways are significant now, compared with before?
Re-run it with the wrong background (all genes in logcpm instead of only the tested ones). 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?

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

## [1] 5

## [1] 5

loose_terms <- as.data.frame(ora_loose)$ID
strict_terms <- as.data.frame(ora_h_up)$ID

setdiff(loose_terms, strict_terms)

## character(0)

setdiff(strict_terms, loose_terms)

## character(0)

intersect(loose_terms, strict_terms)

## [1] "HALLMARK_G2M_CHECKPOINT"            "HALLMARK_INTERFERON_ALPHA_RESPONSE"
## [3] "HALLMARK_E2F_TARGETS"               "HALLMARK_INTERFERON_GAMMA_RESPONSE"
## [5] "HALLMARK_MITOTIC_SPINDLE"

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

wrongbg_df <- as.data.frame(ora_wrongbg)
correct_df <- as.data.frame(ora_h_up)

setdiff(wrongbg_df$ID[wrongbg_df$p.adjust < 0.05], correct_df$ID[correct_df$p.adjust < 0.05])

## character(0)

setdiff(correct_df$ID[correct_df$p.adjust < 0.05], wrongbg_df$ID[wrongbg_df$p.adjust < 0.05])

## character(0)

# 3
ifn_set <- h_sets$gene_symbol[h_sets$gs_name == "HALLMARK_INTERFERON_ALPHA_RESPONSE"]
length(ifn_set)

## [1] 97

length(intersect(ifn_set, genes_up))

## [1] 18

5 Block 5 — GSEA

5.1 Why rank instead of cut

ORA throws away everything below the cut-off. GSEA keeps all genes, ranks them from most up to most down, and asks whether the members of a pathway are clustered at one end of that ranking. So a pathway in which 200 genes each move a little — never reaching significance on their own — can still come out clearly. This is exactly the situation for immune programmes in blood.

5.2 Building the ranked list

Rank by the Wald statistic (stat): it combines the size of the effect and its precision, and it is signed.

rank_df <- as.data.frame(res)              # use the UNshrunken result: it has `stat`
rank_df$gene <- rownames(rank_df)
rank_df <- rank_df[!is.na(rank_df$stat), ] # GSEA cannot use missing values

ranked <- rank_df$stat                     # the values…
names(ranked) <- rank_df$gene              # …with the gene symbols as names
ranked <- sort(ranked, decreasing = TRUE)  # GSEA requires a sorted vector

head(ranked, 8)    # most up in COVID-19

##    IFI27     TPX2    CCNB2     RRM2     TYMS IGHV1-24    MKI67     OTOF 
## 9.486781 7.371511 7.056436 6.905009 6.555583 6.498146 6.435114 6.407280

tail(ranked, 8)    # most down in COVID-19

##   RABGGTB    RPS27A      RPS6     RPL11      RPL6    EEF1A1    RPL7P9     NEAT1 
## -5.498438 -5.504737 -5.533187 -5.576179 -5.586265 -6.040015 -6.118251 -6.222182

length(ranked)

## [1] 17549

5.3 GSEA against Hallmark

gsea_h <- GSEA(
  geneList      = ranked,
  TERM2GENE     = h_sets,
  minGSSize     = 15,        # ignore very small sets: unstable
  maxGSSize     = 500,       # ignore very large sets: uninformative
  pvalueCutoff  = 0.05,
  pAdjustMethod = "BH",
  eps           = 0,         # allow very small p-values instead of rounding to 0
  seed          = TRUE       # reproducible permutations
)

gsea_df <- as.data.frame(gsea_h)
nrow(gsea_df)

## [1] 32

head(gsea_df[, c("ID", "setSize", "NES", "pvalue", "p.adjust")], 10)

NES, the normalised enrichment score, is the number to report. Positive = the pathway sits at the top of the ranking (up in COVID-19); negative = at the bottom (down in COVID-19). It is normalised for set size, so NES values are comparable between pathways.

5.4 GSEA against KEGG

# KEGG needs Entrez IDs, so translate the names of the ranked vector.
map_all <- bitr(names(ranked), fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)

ranked_entrez <- ranked[map_all$SYMBOL]        # keep the genes that converted
names(ranked_entrez) <- map_all$ENTREZID       # rename them to Entrez
ranked_entrez <- ranked_entrez[!duplicated(names(ranked_entrez))]
ranked_entrez <- sort(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")

head(as.data.frame(gsea_kegg)[, c("Description", "setSize", "NES", "p.adjust")], 10)

dotplot(gsea_kegg,  x = "NES", showCategory = 10, split = ".sign") +
  facet_grid(. ~ .sign) +
  theme(axis.text.y = element_text(size = 8)) +
  ggtitle("GSEA, KEGG — COVID-19 vs healthy")

5.5 The enrichment plot

# One pathway. Top panel: the running enrichment score. Middle: where the pathway's
# genes sit in the ranking. Bottom: the ranking metric itself.
gseaplot2(gsea_h,
          geneSetID = "HALLMARK_INTERFERON_ALPHA_RESPONSE",
          title     = "HALLMARK_INTERFERON_ALPHA_RESPONSE",
          pvalue_table = TRUE)

# Several pathways in one panel: use the row numbers of the result table.
sel <- head(order(gsea_df$p.adjust), 3)      # the three best pathways
gsea_df$ID[sel]

## [1] "HALLMARK_G2M_CHECKPOINT"            "HALLMARK_INTERFERON_ALPHA_RESPONSE"
## [3] "HALLMARK_INTERFERON_GAMMA_RESPONSE"

gseaplot2(gsea_h, geneSetID = sel, pvalue_table = FALSE,
          title = "Three most significant Hallmark pathways")

# The "leading edge": the genes that actually produced the enrichment.
le <- gsea_df$core_enrichment[gsea_df$ID == "HALLMARK_INTERFERON_ALPHA_RESPONSE"]
le

## [1] "IFI27/LY6E/USP18/EPSTI1/OAS1/TRIM26/IFI44L/RSAD2/OASL/CMPK2/ISG15/MX1/IFITM3/CNP/TRIM14/LAP3/HERC6/HELZ2/IFIT3/MOV10/PLSCR1/RTP4/LGALS3BP/TMEM140/IRF7/CMTR1/IFI44/CSF1/NUB1/UBE2L6/ADAR/DHX58/IFITM1/BATF2/SELL/IFI35/EIF2AK2/IFIH1/TRAFD1/TAP1/TRIM21/RNF31/LAMP3/DDX60/PARP12/BST2/CCRL2/TRIM25/HLA-C/IFIT2/WARS1/IRF2/PARP9/TDRD7/STAT2/TRIM5/SAMD9L/GMPR"

strsplit(le, "/")[[1]][1:20]     # first 20 leading-edge genes

##  [1] "IFI27"  "LY6E"   "USP18"  "EPSTI1" "OAS1"   "TRIM26" "IFI44L" "RSAD2"  "OASL"   "CMPK2" 
## [11] "ISG15"  "MX1"    "IFITM3" "CNP"    "TRIM14" "LAP3"   "HERC6"  "HELZ2"  "IFIT3"  "MOV10"

ORA vs GSEA. ORA answers “is my list enriched?”; GSEA answers “is this pathway shifted in my ranking?”. They usually agree on the strong signals (interferon) and disagree on the weak, coordinated ones — where GSEA is more sensitive. Report which one you used, with the exact cut-off or ranking metric.

5.6 Exercise 5

5.6.1 Question

Which Hallmark pathway has the most negative NES, and what does that mean biologically for COVID-19 blood?
Draw the enrichment plot for HALLMARK_INFLAMMATORY_RESPONSE, and report its NES and adjusted p-value.
Re-run the Hallmark GSEA ranking by log2FoldChange instead of stat. Do the top three pathways change? Which ranking do you trust more, and why?

5.6.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 <- rank_df$log2FoldChange
names(rank_lfc) <- rank_df$gene
rank_lfc <- sort(rank_lfc, 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 preferable: a large fold change measured on 4 reads should not
# outrank a moderate one measured on 4000.

6 Exporting and reporting

# Full results tables, so that anyone can check your numbers.
write.csv(res_df,                 file.path(res_dir, "DE_COVID19_vs_healthy.csv"), row.names = FALSE)
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),file.path(res_dir, "ORA_kegg_up.csv"),           row.names = FALSE)
write.csv(gsea_df,                file.path(res_dir, "GSEA_hallmark.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.6.2
## 
## 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] tidyr_1.3.2                 dplyr_1.2.1                 ggplot2_4.0.3              
##  [7] msigdbr_26.1.1              org.Hs.eg.db_3.21.0         AnnotationDbi_1.72.0       
## [10] enrichplot_1.28.4           clusterProfiler_4.16.0      apeglm_1.30.0              
## [13] DESeq2_1.48.2               SummarizedExperiment_1.40.0 Biobase_2.70.0             
## [16] MatrixGenerics_1.22.0       matrixStats_1.5.0           GenomicRanges_1.62.1       
## [19] Seqinfo_1.0.0               IRanges_2.44.0              S4Vectors_0.48.1           
## [22] BiocGenerics_0.56.0         generics_0.1.4             
## 
## loaded via a namespace (and not attached):
##   [1] rstudioapi_0.19.0       jsonlite_2.0.0          magrittr_2.0.5         
##   [4] ggtangle_0.1.2          farver_2.1.2            rmarkdown_2.32         
##   [7] fs_2.1.0                vctrs_0.7.3             memoise_2.0.1          
##  [10] ggtree_3.16.3           htmltools_0.5.9         S4Arrays_1.10.1        
##  [13] SparseArray_1.10.10     gridGraphics_0.5-1      sass_0.4.10            
##  [16] bslib_0.12.0            plyr_1.8.9              cachem_1.1.0           
##  [19] igraph_2.3.3            lifecycle_1.0.5         pkgconfig_2.0.3        
##  [22] Matrix_1.7-6            R6_2.6.1                fastmap_1.2.0          
##  [25] gson_0.2.1              GenomeInfoDbData_1.2.14 digest_0.6.39          
##  [28] numDeriv_2016.8-1.1     aplot_0.3.1             colorspace_2.1-3       
##  [31] patchwork_1.3.2         RSQLite_3.53.3          labeling_0.4.3         
##  [34] httr_1.4.9              abind_1.4-8             compiler_4.5.0         
##  [37] bit64_4.8.6             withr_3.0.3             S7_0.2.2               
##  [40] BiocParallel_1.44.0     DBI_1.3.0               R.utils_2.13.0         
##  [43] MASS_7.3-66             rappdirs_0.3.4          DelayedArray_0.36.1    
##  [46] tools_4.5.0             otel_0.2.0              ape_5.8-1              
##  [49] R.oo_1.27.1             glue_1.8.1              nlme_3.1-171           
##  [52] GOSemSim_2.34.0         grid_4.5.0              reshape2_1.4.5         
##  [55] fgsea_1.34.2            gtable_0.3.6            R.methodsS3_1.8.2      
##  [58] data.table_1.18.6.1     XVector_0.50.0          pillar_1.11.1          
##  [61] stringr_1.6.0           yulab.utils_0.2.5       emdbook_1.3.14         
##  [64] splines_4.5.0           treeio_1.32.0           lattice_0.23-1         
##  [67] bit_4.6.0               tidyselect_1.2.1        GO.db_3.21.0           
##  [70] locfit_1.5-9.12         Biostrings_2.78.0       knitr_1.51             
##  [73] gridExtra_2.3.1         xfun_0.60               stringi_1.8.9          
##  [76] UCSC.utils_1.4.0        lazyeval_0.2.3          ggfun_0.2.1            
##  [79] yaml_2.3.12             evaluate_1.0.5          codetools_0.2-20       
##  [82] bbmle_1.0.26            tibble_3.3.1            qvalue_2.40.0          
##  [85] ggplotify_0.1.3         cli_3.6.6               jquerylib_0.1.4        
##  [88] Rcpp_1.1.2              GenomeInfoDb_1.44.3     coda_0.19-4.1          
##  [91] png_0.1-9               bdsmatrix_1.3-7         parallel_4.5.0         
##  [94] assertthat_0.2.1        blob_1.3.0              DOSE_4.2.0             
##  [97] viridisLite_0.4.3       mvtnorm_1.4-2           tidytree_0.4.8         
## [100] scales_1.4.0            purrr_1.2.2             crayon_1.5.3           
## [103] rlang_1.3.0             cowplot_1.2.0           fastmatch_1.1-8        
## [106] KEGGREST_1.50.0

A reporting checklist

Before a figure leaves your computer, you should be able to state:

The accession and the exact samples used (n per group, after which exclusions).
The design formula, the reference level, and the contrast.
The filtering rule, and the number of genes tested.
The significance rule: padj threshold, fold-change threshold, and whether the fold-change threshold was tested or applied afterwards.
Whether fold changes are shrunken, and with which method.
For pathways: the database and its version, the background set, the ranking metric, and the multiple-testing method.
The software versions (sessionInfo()).

Take-home messages

Look at the structure (PCA) before you model. If PC1 is a technical variable, fix that first.
DESeq2 needs raw counts; normalisation is inside the model, not before it.
The design formula is the analysis. Putting age or sex in it is a scientific decision, not a formality.
padj = NA is information, not a bug.
ORA needs a correct background; GSEA needs a well-chosen ranking metric.
Always check a positive control — here, the interferon response — before you believe anything novel in the list.

Försöksdjursvetenskap (teoretisk) och avancerad databearbetning

Lab 5 — PCA, differential expression with DESeq2, and pathway analysis (ORA and GSEA)

Birong Zhang

21 October 2026