How today is organised

Four blocks, four hours. Each block ends with a short exercise you do on your own before we move on.

	Block	Time	Exercise
1	Design matrices and DESeq2	75 min	Ex. 1 (design), Ex. 2 (diagnostics)
	break	15 min
2	edgeR	45 min	Ex. 3 (normalisation & test choice)
3	Comparing DESeq2 and edgeR	40 min	Ex. 4 (who disagrees, and why)
	break	10 min
4	Pathway analysis: ORA and GSEA	65 min	Ex. 5 (background & ranking)

0.1 Learning objectives

Translate a biological question into a design matrix and read the coefficients back out.
Run DESeq2 step by step — normalisation, dispersion, testing, shrinkage — and interpret every diagnostic plot.
Run the same contrast in edgeR, and explain how TMM and quasi-likelihood differ from median-of-ratios and the Wald test.
Quantify the agreement between two DE pipelines and characterise the genes they disagree about.
Perform over-representation analysis and GSEA correctly, with particular attention to the background set.

0.2 Packages

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

library(SummarizedExperiment)
library(DESeq2)
library(edgeR)
library(ggplot2)
library(dplyr)
library(tidyr)
library(ggrepel)
library(pheatmap)
library(RColorBrewer)
library(matrixStats)
library(gridExtra)
library(apeglm)
library(UpSetR)

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

0.3 Loading yesterday’s object

se <- readRDS(file.path(data_dir, "se_raw.rds"))
se

## class: SummarizedExperiment 
## dim: 63152 216 
## metadata(0):
## assays(1): counts
## rownames(63152): ENSG00000000003 ENSG00000000005 ... ENSG00000272545 TC%
## rowData names(8): ensembl_gene_id external_gene_name ... gene_strand symbol
## colnames(216): L400T L401T ... L887T L890T
## colData names(14): geo_accession title ... ps_who sex

table(colData(se)$histology, colData(se)$sex)

##        
##         female male
##   LUAD      68   38
##   LUSC      24   43
##   other     10   14

head(rowData(se))

## DataFrame with 6 rows and 8 columns
##                 ensembl_gene_id external_gene_name   gene_biotype chromosome_name start_position
##                     <character>        <character>    <character>     <character>      <integer>
## ENSG00000000003 ENSG00000000003             TSPAN6 protein_coding               X       99883667
## ENSG00000000005 ENSG00000000005               TNMD protein_coding               X       99839799
## ENSG00000000419 ENSG00000000419               DPM1 protein_coding              20       49551404
## ENSG00000000457 ENSG00000000457              SCYL3 protein_coding               1      169821804
## ENSG00000000460 ENSG00000000460           C1orf112 protein_coding               1      169631245
## ENSG00000000938 ENSG00000000938                FGR protein_coding               1       27938575
##                 end_position gene_strand      symbol
##                    <integer> <character> <character>
## ENSG00000000003     99894988           -      TSPAN6
## ENSG00000000005     99854882           +        TNMD
## ENSG00000000419     49575092           -        DPM1
## ENSG00000000457    169863408           -       SCYL3
## ENSG00000000460    169823221           +    C1orf112
## ENSG00000000938     27961788           -         FGR

Today’s question: which genes distinguish lung adenocarcinoma (LUAD) from lung squamous cell carcinoma (LUSC)? We verified those labels against marker genes yesterday, so we know they mean what they say.

keep_samples <- colData(se)$histology %in% c("LUAD", "LUSC")
se2 <- se[, keep_samples]
colData(se2)$histology <- relevel(droplevels(factor(colData(se2)$histology)), ref = "LUAD")

age_var <- grep("^age", colnames(colData(se2)), value = TRUE)[1]
colData(se2)$age <- as.numeric(colData(se2)[[age_var]])
colData(se2)$age[is.na(colData(se2)$age)] <- median(colData(se2)$age, na.rm = TRUE)

dim(se2); table(colData(se2)$histology)

## [1] 63152   173

## 
## LUAD LUSC 
##  106   67

# If the room is slow, subsample to 20 + 20 and everything below still runs in seconds.
# set.seed(1)
# idx <- c(sample(which(colData(se2)$histology == "LUAD"), 20),
#          sample(which(colData(se2)$histology == "LUSC"), 20))
# se2 <- se2[, idx]

1 Block 1 — Design matrices and DESeq2

1.1 What the formula actually means

~ sex + histology says: model each gene’s log-expression as a baseline, plus a shift for sex, plus a shift for histology — and report the histology effect after adjusting for sex.

cd <- as.data.frame(colData(se2))
design <- model.matrix(~ sex + histology, data = cd)
head(design)

##       (Intercept) sexmale histologyLUSC
## L400T           1       1             0
## L401T           1       1             0
## L404T           1       0             0
## L406T           1       1             1
## L413T           1       0             0
## L420T           1       0             1

colnames(design)

## [1] "(Intercept)"   "sexmale"       "histologyLUSC"

(Intercept) — mean of the reference group (female, LUAD)
sexmale — male vs female shift
histologyLUSC — the coefficient we want

qr(design)$rank == ncol(design)     # full rank? if FALSE, two variables are aliased

## [1] TRUE

table(cd$sex, cd$histology)         # covariate balance across the contrast

##         
##          LUAD LUSC
##   female   68   24
##   male     38   43

Discuss. If every LUSC sample were male and every LUAD sample female, sex and histology would be perfectly confounded, the design would be rank-deficient, and no statistical method could separate them. Partial imbalance is fine — that is exactly what adjustment handles. Check this cross-tab before running anything.

~ condition                # simplest two-group comparison
~ batch + condition        # adjust for a known technical batch
~ patient + condition      # paired design (tumour/normal from the same patient)
~ condition * time         # main effects plus interaction
~ 0 + group                # cell means; build contrasts explicitly with makeContrasts()

1.2 Filtering

counts <- assay(se2, "counts")
dge0 <- DGEList(counts, group = cd$histology)
keep <- filterByExpr(dge0, design = design)
table(keep)

## keep
## FALSE  TRUE 
## 39098 24054

se_f <- se2[keep, ]
dim(se_f)

## [1] 24054   173

Filtering before testing is not cheating, provided the filter is independent of the group labels under the null. filterByExpr uses overall expression level and group sizes, never group means, so it qualifies. Dropping roughly half the rows buys real power back at the multiple-testing step.

1.3 Building the DESeq2 object

dds <- DESeqDataSetFromMatrix(
  countData = assay(se_f, "counts"),
  colData   = colData(se_f),
  design    = ~ sex + histology
)
rowData(dds) <- rowData(se_f)
dds

## class: DESeqDataSet 
## dim: 24054 173 
## metadata(1): version
## assays(1): counts
## rownames(24054): ENSG00000000003 ENSG00000000419 ... ENSG00000272540 TC%
## rowData names(8): ensembl_gene_id external_gene_name ... gene_strand symbol
## colnames(173): L400T L401T ... L887T L890T
## colData names(14): geo_accession title ... ps_who sex

1.4 Step 1 — size factors (median of ratios)

DESeq2 does not normalise by total counts. It uses the median of ratios to a pseudo-reference sample, which is robust to a few very highly expressed genes dominating a library.

dds <- estimateSizeFactors(dds)
head(sizeFactors(dds))

##     L400T     L401T     L404T     L406T     L413T     L420T 
## 0.9519501 1.0155416 1.0852169 1.1794930 0.9551240 0.9085463

class(dds)

## [1] "DESeqDataSet"
## attr(,"package")
## [1] "DESeq2"

rowRanges(dds)

## GRangesList object of length 24054:
## $ENSG00000000003
## GRanges object with 0 ranges and 0 metadata columns:
##    seqnames    ranges strand
##       <Rle> <IRanges>  <Rle>
##   -------
##   seqinfo: no sequences
## 
## $ENSG00000000419
## GRanges object with 0 ranges and 0 metadata columns:
##    seqnames    ranges strand
##       <Rle> <IRanges>  <Rle>
##   -------
##   seqinfo: no sequences
## 
## $ENSG00000000457
## GRanges object with 0 ranges and 0 metadata columns:
##    seqnames    ranges strand
##       <Rle> <IRanges>  <Rle>
##   -------
##   seqinfo: no sequences
## 
## ...
## <24051 more elements>

colnames(mcols(rowRanges(dds)))

## [1] "ensembl_gene_id"    "external_gene_name" "gene_biotype"       "chromosome_name"   
## [5] "start_position"     "end_position"       "gene_strand"        "symbol"

# By hand, to see there is no magic in it:
cts <- counts(dds)
nz  <- cts[rowSums(cts == 0) == 0, , drop = FALSE]
loggeo <- rowMeans(log(nz))
manual_sf <- apply(nz, 2, function(x) exp(median(log(x) - loggeo)))

par(mfrow = c(1, 2))
plot(sizeFactors(dds), manual_sf, pch = 16, cex = 0.6,
     xlab = "DESeq2 size factor", ylab = "manual median-of-ratios"); abline(0, 1, col = "red")
plot(colSums(cts) / 1e6, sizeFactors(dds), pch = 16, cex = 0.6,
     xlab = "library size (millions)", ylab = "size factor", main = "Size factor vs depth")

par(mfrow = c(1, 1))

If the right-hand panel were a perfect line, CPM and DESeq2 normalisation would agree exactly. The scatter around it is the composition effect: samples where a handful of genes soak up a large share of the reads.

1.5 Step 2 — dispersion

Counts are not Poisson; biological replicates are overdispersed. DESeq2 fits a mean–dispersion trend across all genes and shrinks each gene’s noisy estimate toward it. This is what makes a 3-vs-3 experiment analysable at all.

dds <- estimateDispersions(dds)
plotDispEsts(dds, main = "Dispersion estimates")

Reading the plot:

black — per-gene maximum-likelihood estimates, very noisy at low counts
red — the fitted mean–dispersion trend
blue — final shrunken values, used for testing
blue circles — dispersion outliers, deliberately not shrunk (shrinking them would create false positives)

A healthy plot decreases with the mean and hugs the trend. A cloud sitting far above the line usually means an unmodelled batch or an outlier sample.

1.6 Step 3 — testing

dds <- nbinomWaldTest(dds)
resultsNames(dds)

## [1] "Intercept"              "sex_male_vs_female"     "histology_LUSC_vs_LUAD"

res <- results(dds, name = "histology_LUSC_vs_LUAD", alpha = 0.05)
# res = results(dds, contrast=c("histology", "LUSC","LUAD"), alpha=0.05)
summary(res)

## 
## out of 24054 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 6501, 27%
## LFC < 0 (down)     : 6220, 26%
## outliers [1]       : 1053, 4.4%
## low counts [2]     : 0, 0%
## (mean count < 5)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

res_df <- as.data.frame(res)
res_df$ensembl <- rownames(res_df)
res_df$symbol  <- rowData(dds)$symbol[match(rownames(res_df), rownames(dds))]
res_df <- res_df[order(res_df$padj), ]
head(res_df[, c("symbol", "baseMean", "log2FoldChange", "pvalue", "padj")], 15)

Sanity check. Are the top genes what yesterday’s marker plots predicted? KRT5, TP63, DSG3 on one side and NAPSA, NKX2-1, SFTPC on the other means the pipeline is behaving. Always build a positive control into an analysis.

1.7 Independent filtering and `padj = NA`

table(is.na(res$padj))

## 
## FALSE  TRUE 
## 23001  1053

metadata(res)$filterThreshold

##       0% 
## 4.742495

plot(metadata(res)$filterNumRej, type = "b",
     xlab = "quantile of baseMean filtered", ylab = "number of rejections",
     main = "Independent filtering optimisation")

padj = NA means the gene was not tested, not that it was non-significant. Two causes: independent filtering judged it too lowly expressed to ever reach significance, or Cook’s distance flagged it as driven by a count outlier.

1.8 Log fold change shrinkage

A raw log2FoldChange for a low-count gene is unreliable — 3 counts vs 0 gives a huge fold change on no evidence. Shrinkage pulls estimates toward zero in proportion to their uncertainty.

res_shr <- lfcShrink(dds, coef = "histology_LUSC_vs_LUAD", type = "apeglm")

par(mfrow = c(1, 2))
plotMA(res,     ylim = c(-8, 8), main = "MA plot: raw LFC")
plotMA(res_shr, ylim = c(-8, 8), main = "MA plot: apeglm-shrunken LFC")

par(mfrow = c(1, 1))

Rule. Unshrunken results for p-values and significance ranking; shrunken LFCs for plotting, for ranking by effect size, and as the GSEA input. Never filter on abs(log2FoldChange) > 1 using unshrunken values.

1.9 Multiple testing

par(mfrow = c(1, 2))
hist(res$pvalue, breaks = 50, col = "grey80", xlab = "raw p-value",
     main = "p-value distribution")
hist(res$pvalue[res$baseMean > 50], breaks = 50, col = "grey80", xlab = "raw p-value",
     main = "well-expressed genes only")

par(mfrow = c(1, 1))

What the shape tells you:

flat with a spike at 0 — the expected picture: uniform nulls plus real signal
completely flat — no differential expression, or badly underpowered
hill-shaped — usually a misspecified model: unmodelled batch, or correlated samples
peak near 1 — over-conservative variance estimates, sometimes from over-filtering

c(raw_p05    = sum(res$pvalue < 0.05, na.rm = TRUE),
  BH_FDR05   = sum(res$padj   < 0.05, na.rm = TRUE),
  bonferroni = sum(p.adjust(res$pvalue, "bonferroni") < 0.05, na.rm = TRUE))

##    raw_p05   BH_FDR05 bonferroni 
##      13640      12721       5346

BH controls the expected proportion of false discoveries among the genes you call significant. Bonferroni controls the probability of any false positive. For discovery-stage transcriptomics BH is the right trade-off.

1.10 Testing against a fold-change threshold

An FDR-significant gene with a 1.05-fold change is real and irrelevant. Rather than post-hoc filtering on effect size (which invalidates the FDR), test the composite null directly:

res_thr <- results(dds, name = "histology_LUSC_vs_LUAD", lfcThreshold = 1, alpha = 0.05)
sum(res_thr$padj < 0.05, na.rm = TRUE)

## [1] 1422

1.11 Visualising the DESeq2 result

vol <- as.data.frame(res_shr)
vol$symbol   <- rowData(dds)$symbol[match(rownames(vol), rownames(dds))]
vol$padj_raw <- res$padj[match(rownames(vol), rownames(res))]
vol <- vol[!is.na(vol$padj_raw), ]
vol$class <- with(vol, ifelse(padj_raw < 0.05 & log2FoldChange >  1, "up in LUSC",
                       ifelse(padj_raw < 0.05 & log2FoldChange < -1, "up in LUAD", "ns")))
table(vol$class)

## 
##         ns up in LUAD up in LUSC 
##      20159       1312       1530

lab <- vol %>% dplyr::filter(class != "ns") %>% group_by(class) %>%
  slice_max(abs(log2FoldChange) * -log10(padj_raw + 1e-300), n = 12)

ggplot(vol, aes(log2FoldChange, -log10(padj_raw), colour = class)) +
  geom_point(size = 0.7, alpha = 0.6) +
  geom_vline(xintercept = c(-1, 1), linetype = 2, colour = "grey40") +
  geom_hline(yintercept = -log10(0.05), linetype = 2, colour = "grey40") +
  geom_text_repel(data = lab, aes(label = symbol), size = 3, max.overlaps = 30,
                  show.legend = FALSE) +
  scale_colour_manual(values = c(ns = "grey75", "up in LUSC" = "#D7301F",
                                 "up in LUAD" = "#2B8CBE")) +
  labs(x = "shrunken log2 fold change (LUSC / LUAD)", y = "-log10 adjusted p",
       colour = NULL, title = "LUSC vs LUAD, DESeq2")

vsd <- vst(dds, blind = FALSE)   # blind = FALSE: use the design when fitting the trend
plotPCA(vsd, intgroup = c("histology", "sex")) +
  aes(colour = histology, shape = sex) +
  labs(
    title = "PCA on variance-stabilised counts",
    colour = "Histology",
    shape = "Sex"
  )

top_genes <- rownames(res_df)[1:40]
mat <- assay(vsd)[top_genes, ]
rownames(mat) <- rowData(dds)$symbol[match(top_genes, rownames(dds))]
mat_z <- t(scale(t(mat)))

pheatmap(mat_z,
         annotation_col = data.frame(histology = colData(vsd)$histology,
                                     sex = colData(vsd)$sex,
                                     row.names = colnames(mat_z)),
         annotation_colors = list(histology = c(LUAD = "#2B8CBE", LUSC = "#D7301F"),
                                  sex = c(female = "#8DA0CB", male = "#E78AC3")),
         show_colnames = FALSE, fontsize_row = 8, clustering_method = "ward.D2",
         breaks = seq(-3, 3, length.out = 101),
         color = colorRampPalette(rev(brewer.pal(11, "RdBu")))(100),
         main = "Top 40 DE genes (vst, row z-scored)")

Careful. This heatmap will always look convincing, because the genes were chosen because they separate the groups. It illustrates a result; it is not evidence for one. Evidence would be the same genes separating an independent cohort.

1.12 Exercise 1 — does the covariate matter?

~10 minutes.

1.12.1 Question

Refit DESeq2 without sex in the design (~ histology only).
How many genes change their significance call at FDR 5%?
Which chromosomes are those genes on? Explain what you see.
Does the histologyLUSC fold change itself change much? Plot the two against each other.

1.12.2 Solution

dds_nosex <- DESeqDataSetFromMatrix(assay(se_f, "counts"), colData(se_f), ~ histology)
dds_nosex <- DESeq(dds_nosex, quiet = TRUE)
res_nosex <- results(dds_nosex, name = "histology_LUSC_vs_LUAD")

ok <- !is.na(res$padj) & !is.na(res_nosex$padj)
changed <- rownames(res)[ok & xor(res$padj < 0.05, res_nosex$padj < 0.05)]
length(changed)

## [1] 1121

sort(table(rowData(dds)$chromosome_name[match(changed, rownames(dds))]),decreasing = TRUE)

## 
##   1  19   2   X   7  12  17   3  11  14   6  10   5   9   4  16   8  20  15  18  22  13  21   Y  MT 
## 121  83  72  66  65  60  54  53  50  50  49  48  48  44  41  39  37  32  31  23  23  14  14   3   1

all_chr <- table(rowData(dds)$chromosome_name)

changed_chr <- table(
  rowData(dds)$chromosome_name[
    match(changed, rownames(dds))
  ]
)

df<-data.frame(
  chromosome = names(changed_chr),
  changed = as.integer(changed_chr),
  total = as.integer(all_chr[names(changed_chr)]),
  proportion = as.integer(changed_chr) /
               as.integer(all_chr[names(changed_chr)])
)

df <- df %>% arrange(desc(proportion))
head(df,5)

plot(res$log2FoldChange, res_nosex$log2FoldChange, pch = ".",
     xlab = "with sex in design", ylab = "without sex",
     main = "log2 fold change, LUSC vs LUAD"); abline(0, 1, col = "red")

X and Y should be over-represented among the genes whose call changes: without a sex term, sex-linked variation goes into the residual, inflating the variance estimate for those genes. The fold changes barely move — adding a balanced-ish covariate mostly changes the precision, not the estimate.

1.13 Exercise 2 — reading the diagnostics

~10 minutes.

1.13.1 Question

Take the gene with the largest raw log2FoldChange (unshrunken) among genes with baseMean < 20. Plot its counts with plotCounts(). Would you believe it?
What is its shrunken LFC? By what factor did it shrink?
Pick any gene with padj = NA and work out which of the two causes applies.
Re-plot the p-value histogram using only genes with baseMean < 10. Why does it look different?

1.13.2 Solution

low <- res[which(res$baseMean < 20), ]
g <- rownames(low)[which.max(abs(low$log2FoldChange))]
rowData(dds)$symbol[match(g, rownames(dds))]

## [1] "KPRP"

par(mfrow = c(1, 2))
plotCounts(dds, gene = g, intgroup = "histology", main = "raw counts")
hist(res$pvalue[res$baseMean < 10], breaks = 50, col = "grey80",
     main = "p-values, baseMean < 10", xlab = "raw p-value")

par(mfrow = c(1, 1))

c(raw = res[g, "log2FoldChange"], shrunk = res_shr[g, "log2FoldChange"])

##                    raw shrunk.ENSG00000203786 
##               7.797589               9.806862

# padj = NA: filtered, or Cook's outlier?
na_genes <- rownames(res)[is.na(res$padj)]
data.frame(baseMean = res[na_genes[1:5], "baseMean"],
           below_filter = res[na_genes[1:5], "baseMean"] < metadata(res)$filterThreshold,
           pvalue_is_NA = is.na(res[na_genes[1:5], "pvalue"]))   # TRUE => Cook's outlier

The low-count gene is usually driven by two or three samples; shrinkage typically pulls it several-fold toward zero. The low-expression p-value histogram is concentrated near 1 because discrete counts give few attainable p-values and the test has almost no power — which is precisely why independent filtering removes these genes before FDR correction.

2 Block 2 — edgeR

Same generalised linear model family, different implementation choices at every step. Running edgeR is not redundancy — it is the cheapest robustness check you can buy.

Step	DESeq2	edgeR
Normalisation	median of ratios (size factors)	TMM (normalisation factors × library size)
Dispersion	gene-wise, shrunk to a fitted trend (empirical Bayes)	common → trended → tagwise, shrunk toward the trend
Test	Wald on the coefficient (or LRT)	quasi-likelihood F test (`glmQLFTest`) or LRT (`glmLRT`)
Extra uncertainty	dispersion outlier handling, Cook’s distance	QL step models the uncertainty in the dispersion itself
Filtering	independent filtering inside `results()`	`filterByExpr()`, done up front by the user

2.1 The object and TMM normalisation

y <- DGEList(counts = assay(se_f, "counts"), group = cd$histology)
y <- calcNormFactors(y, method = "TMM")
head(y$samples)

TMM picks a reference sample, trims the extreme log-ratios and extreme abundances, and takes a weighted mean of what remains. The result is one factor per sample that corrects for composition — the same problem DESeq2’s median-of-ratios addresses, by a different route.

eff_lib <- y$samples$lib.size * y$samples$norm.factors
plot(sizeFactors(dds), eff_lib / mean(eff_lib), pch = 16, cex = 0.6,
     xlab = "DESeq2 size factor", ylab = "edgeR effective library size (scaled)",
     main = sprintf("r = %.3f", cor(sizeFactors(dds), eff_lib)))

These two normalisations almost always agree closely. When they do not, look for samples dominated by a few transcripts — haemoglobin, immunoglobulin, mitochondrial RNA.

2.2 Dispersion

y <- estimateDisp(y, design, robust = TRUE)
c(common_dispersion = y$common.dispersion,
  BCV = sqrt(y$common.dispersion))

## common_dispersion               BCV 
##         0.5288045         0.7271894

plotBCV(y)

The biological coefficient of variation is the square root of the dispersion: a BCV of 0.4 means expression varies by about 40% between biological replicates, independent of depth. Typical values: 0.4 for human patient samples, 0.1 for inbred mouse lines, ~0.01 for technical replicates. This one number tells you more about your power than the sample size does.

2.3 Quasi-likelihood test

fit <- glmQLFit(y, design, robust = TRUE)
plotQLDisp(fit)

qlf <- glmQLFTest(fit, coef = "histologyLUSC")
topTags(qlf, n = 10)

## Coefficient:  histologyLUSC 
##                    logFC     logCPM         F        PValue           FDR
## ENSG00000154227 9.088341  3.2264918 2114.9909 5.570459e-150 1.339918e-145
## ENSG00000187054 9.117010  3.3504197  922.5067  7.225084e-99  8.689608e-95
## ENSG00000260581 5.554466  0.5900447  674.4819  2.189223e-76  1.755319e-72
## ENSG00000224984 6.401246 -0.3760469  660.2547  4.297141e-76  2.584086e-72
## ENSG00000197641 9.419372  5.4282822  826.3343  1.281234e-73  6.163759e-70
## ENSG00000188508 7.361681  1.3450667  554.3536  1.196431e-72  4.796491e-69
## ENSG00000186081 9.366847  9.8620191  953.9183  1.628544e-72  5.596143e-69
## ENSG00000125998 9.093402  1.9542627  673.5159  3.104054e-70  9.333114e-67
## ENSG00000169594 7.031271  4.3652460  812.3110  2.524177e-67  6.746284e-64
## ENSG00000188373 8.251756  1.3966540  560.6972  7.234132e-67  1.740098e-63

res_edger <- topTags(qlf, n = Inf, sort.by = "none")$table
res_edger$symbol <- rowData(se_f)$symbol
sum(res_edger$FDR < 0.05)

## [1] 13520

The QL F-test differs from a likelihood ratio test in that it accounts for the uncertainty in the dispersion estimate itself, not just the dispersion value. It is more conservative and better calibrated, and it is the recommended default for designed experiments. glmLRT() is faster and slightly more powerful when you have many replicates; use it when n is large or dispersion is well estimated.

fit_lrt <- glmFit(y, design)
lrt <- glmLRT(fit_lrt, coef = "histologyLUSC")
c(QL = sum(res_edger$FDR < 0.05),
  LRT = sum(topTags(lrt, n = Inf)$table$FDR < 0.05))

##    QL   LRT 
## 13520 13174

2.4 Fold-change-thresholded testing

edgeR’s analogue of lfcThreshold, and a genuinely better tool than filtering on fold change afterwards:

tr <- glmTreat(fit, coef = "histologyLUSC", lfc = 1)
sum(topTags(tr, n = Inf)$table$FDR < 0.05)

## [1] 1981

2.5 edgeR’s own MD plot

plotMD(qlf, main = "edgeR: mean-difference plot")

2.6 Exercise 3 — normalisation and test choice

~10 minutes.

2.6.1 Question

Rebuild y without calcNormFactors() (i.e. leave norm.factors = 1). Re-estimate dispersion and re-test. How many DE genes do you get, and how does the top-10 list change?
Compare the number of DE genes from glmQLFTest, glmLRT and glmTreat(lfc = 1). Order them from most to least conservative and explain the ordering.
What is the BCV of this dataset? What would you expect it to be for a cell-line experiment, and what does that imply for the number of replicates you would need?

2.6.2 Solution

# 1
y_nonorm <- DGEList(assay(se_f, "counts"), group = cd$histology)
y_nonorm <- estimateDisp(y_nonorm, design, robust = TRUE)
qlf_nonorm <- glmQLFTest(glmQLFit(y_nonorm, design, robust = TRUE), coef = "histologyLUSC")
tt_nonorm <- topTags(qlf_nonorm, n = Inf, sort.by = "none")$table

c(TMM = sum(res_edger$FDR < 0.05), no_norm = sum(tt_nonorm$FDR < 0.05))

##     TMM no_norm 
##   13520   13400

data.frame(TMM     = res_edger$symbol[order(res_edger$FDR)][1:10],
           no_norm = rowData(se_f)$symbol[order(tt_nonorm$FDR)][1:10])

# 2
c(glmLRT = sum(topTags(lrt, n = Inf)$table$FDR < 0.05),
  glmQLFTest = sum(res_edger$FDR < 0.05),
  glmTreat_lfc1 = sum(topTags(tr, n = Inf)$table$FDR < 0.05))

##        glmLRT    glmQLFTest glmTreat_lfc1 
##         13174         13520          1981

# 3
sqrt(y$common.dispersion)

## [1] 0.7271894

glmLRT > glmQLFTest > glmTreat in number of calls: the LRT ignores dispersion uncertainty, QL accounts for it, and glmTreat additionally demands a minimum effect size. Without TMM, the gene list shifts modestly — composition effects are moderate here — but the ranking of borderline genes changes, and in a dataset with a dominant transcript the effect would be dramatic. A BCV around 0.1 (cell lines) instead of 0.4 (patients) means you need roughly an order of magnitude fewer replicates for the same power.

3 Block 3 — Comparing DESeq2 and edgeR

We now have two result tables for the same contrast on the same counts with the same design. How much do they agree, and where do they part company?

stopifnot(identical(rownames(res), rownames(res_edger)))

cmp <- data.frame(
  gene       = rownames(res),
  symbol     = rowData(se_f)$symbol,
  baseMean   = res$baseMean,
  lfc_deseq  = res$log2FoldChange,
  lfc_edger  = res_edger$logFC,
  p_deseq    = res$pvalue,
  p_edger    = res_edger$PValue,
  padj_deseq = res$padj,
  fdr_edger  = res_edger$FDR,
  stringsAsFactors = FALSE
)
cmp <- cmp[!is.na(cmp$padj_deseq) & !is.na(cmp$fdr_edger), ]
nrow(cmp)

## [1] 23001

3.1 Correlation of effect sizes

r_p <- cor(cmp$lfc_deseq, cmp$lfc_edger)
r_s <- cor(cmp$lfc_deseq, cmp$lfc_edger, method = "spearman")
c(pearson = r_p, spearman = r_s)

##   pearson  spearman 
## 0.9992235 0.9996272

ggplot(cmp, aes(lfc_deseq, lfc_edger)) +
  geom_point(size = 0.4, alpha = 0.3) +
  geom_abline(slope = 1, intercept = 0, colour = "red") +
  geom_smooth(method = "lm", se = FALSE, linetype = 2, colour = "blue") +
  labs(x = "DESeq2 log2FC", y = "edgeR log2FC",
       title = sprintf("Fold-change agreement (Pearson r = %.4f)", r_p))

Fold changes are near-identical, because both fit essentially the same negative binomial GLM. The estimates come from the model; the disagreement lives in the variance estimation and therefore in the p-values.

ggplot(cmp, aes(-log10(p_deseq), -log10(p_edger))) +
  geom_point(size = 0.4, alpha = 0.3) +
  geom_abline(slope = 1, intercept = 0, colour = "red") +
  coord_cartesian(xlim = c(0, 50), ylim = c(0, 50)) +
  labs(x = "-log10 p, DESeq2", y = "-log10 p, edgeR",
       title = sprintf("p-value agreement (Spearman rho = %.3f)",
                       cor(cmp$p_deseq, cmp$p_edger, method = "spearman")))

Points systematically below the diagonal mean edgeR’s QL test is the more conservative of the two at that end of the range.

3.2 Overlap of the significant sets

sig <- list(
  DESeq2 = cmp$gene[cmp$padj_deseq < 0.05],
  edgeR  = cmp$gene[cmp$fdr_edger  < 0.05]
)
sapply(sig, length)

## DESeq2  edgeR 
##  12721  12900

both   <- intersect(sig$DESeq2, sig$edgeR)
only_d <- setdiff(sig$DESeq2, sig$edgeR)
only_e <- setdiff(sig$edgeR,  sig$DESeq2)

c(both = length(both), DESeq2_only = length(only_d), edgeR_only = length(only_e),
  jaccard = length(both) / length(union(sig$DESeq2, sig$edgeR)))

##         both  DESeq2_only   edgeR_only      jaccard 
## 1.258200e+04 1.390000e+02 3.180000e+02 9.649513e-01

3.3 Who are the disagreements?

This is the part worth doing carefully. If method-specific genes were scattered randomly through the ranking, you would have a problem. They should instead sit right at the significance boundary. The genes detected by both DESeq2 and edgeR tend to have higher expression and larger effect sizes, whereas genes detected by only one method tend to have smaller fold changes. In particular, edgeR-only genes show substantially lower expression, suggesting that differences between the methods are more pronounced for lowly expressed genes.

cmp$call <- with(cmp, ifelse(gene %in% both, "both",
                      ifelse(gene %in% only_d, "DESeq2 only",
                      ifelse(gene %in% only_e, "edgeR only", "neither"))))
table(cmp$call)

## 
##        both DESeq2 only  edgeR only     neither 
##       12582         139         318        9962

sub <- cmp[cmp$call != "neither", ]

tapply(log10(sub$baseMean), sub$call, mean)

##        both DESeq2 only  edgeR only 
##    2.644470    2.887800    1.809384

tapply(log10(sub$baseMean), sub$call, median)

##        both DESeq2 only  edgeR only 
##    2.814546    2.960342    1.314448

tapply(abs(sub$lfc_deseq), sub$call, mean, na.rm = TRUE)

##        both DESeq2 only  edgeR only 
##   0.8534352   0.3027026   0.4092657

tapply(abs(sub$lfc_deseq), sub$call, median, na.rm = TRUE)

##        both DESeq2 only  edgeR only 
##   0.5548459   0.2260567   0.3561775

g1 <- ggplot(sub, aes(call, log10(baseMean), fill = call)) +
  geom_boxplot(outlier.size = 0.4) + labs(x = NULL, y = "log10 mean expression") +
  theme(legend.position = "none", axis.text.x = element_text(angle = 20, hjust = 1))
g2 <- ggplot(sub, aes(call, abs(lfc_deseq), fill = call)) +
  geom_boxplot(outlier.size = 0.4) + coord_cartesian(ylim = c(0, 3)) +
  labs(x = NULL, y = "|log2 fold change|") +
  theme(legend.position = "none", axis.text.x = element_text(angle = 20, hjust = 1))
grid.arrange(g1, g2, nrow = 1)

# Where do the method-specific genes sit in the *other* method's ranking?
cmp$rank_deseq <- rank(cmp$p_deseq)
cmp$rank_edger <- rank(cmp$p_edger)

summary(cmp$rank_edger[cmp$call == "DESeq2 only"])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   12901   13001   13118   13193   13257   16805

summary(cmp$rank_deseq[cmp$call == "edgeR only"])

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   12723   12878   13073   13317   13455   18694

nrow(cmp)

## [1] 23001

3.4 Concordance at the top

A cleaner summary than a single Venn diagram: for every list length n, what fraction of the top n genes is shared?

ord_d <- cmp$gene[order(cmp$p_deseq)]
ord_e <- cmp$gene[order(cmp$p_edger)]

ns <- seq(50, 3000, by = 50)
cat_df <- data.frame(
  n = ns,
  concordance = sapply(ns, function(n) length(intersect(ord_d[1:n], ord_e[1:n])) / n)
)

ggplot(cat_df, aes(n, concordance)) +
  geom_line(linewidth = 0.8, colour = "steelblue") +
  geom_hline(yintercept = 1, linetype = 2, colour = "grey50") +
  coord_cartesian(ylim = c(0, 1.02)) +
  labs(x = "top n genes by p-value", y = "fraction shared",
       title = "Concordance at the top: DESeq2 vs edgeR")

3.5 What to do with two lists

consensus <- cmp[cmp$padj_deseq < 0.05 & cmp$fdr_edger < 0.05 &
                   abs(cmp$lfc_deseq) > 1, ]
nrow(consensus)

## [1] 3059

head(consensus[order(consensus$padj_deseq), c("symbol", "baseMean", "lfc_deseq",
                                              "padj_deseq", "fdr_edger")], 10)

Practical guidance:

Pick one method in advance and report it as your primary analysis. Choosing the tool after seeing which gives more hits is p-hacking with extra steps.
Use the second method as a robustness statement: “of the N genes called by DESeq2, M (x%) were also significant in edgeR.”
If a headline conclusion depends on which of the two you ran, you do not have a result — you have a borderline gene, and it needs validation.
For downstream pathway analysis, use a single method’s full ranking. Do not intersect first: intersecting truncates the ranking and biases the gene set input.

Optional: a third opinion from limma-voom

library(limma)
v <- voom(y, design, plot = TRUE)
fit_v <- eBayes(lmFit(v, design))
res_voom <- topTable(fit_v, coef = "histologyLUSC", number = Inf, sort.by = "none")
sum(res_voom$adj.P.Val < 0.05)

voom estimates a precision weight per observation from the mean–variance trend and then uses limma’s linear-model machinery. It is fast, handles complex designs gracefully, and is the usual choice when n is large.

3.6 Exercise 4 — characterise the disagreement

~15 minutes.

3.6.1 Question

How many genes are called by DESeq2 only, and how many by edgeR only? What is the Jaccard index of the two significant sets?
Take the 5 most significant “DESeq2 only” genes. What are their edgeR FDR values? Are they far from 0.05, or just over the line?
Compare the baseMean distribution of “both” vs the method-specific genes. Which end of the expression range do the disagreements come from?
Re-do the overlap at FDR 0.01 and FDR 0.10. Does agreement improve or worsen as the threshold gets stricter, and why?

3.6.2 Solution

# 1
c(DESeq2_only = length(only_d), edgeR_only = length(only_e),
  jaccard = length(both) / length(union(sig$DESeq2, sig$edgeR)))

## DESeq2_only  edgeR_only     jaccard 
## 139.0000000 318.0000000   0.9649513

# 2
d_only <- cmp[cmp$call == "DESeq2 only", ]
head(d_only[order(d_only$padj_deseq), c("symbol", "baseMean", "lfc_deseq",
                                        "padj_deseq", "fdr_edger")], 5)

# 3
tapply(log10(sub$baseMean), sub$call, median)

##        both DESeq2 only  edgeR only 
##    2.814546    2.960342    1.314448

# 4
thresholds <- c(0.01, 0.05, 0.10)
do.call(rbind, lapply(thresholds, function(a) {
  sd_ <- cmp$gene[cmp$padj_deseq < a]; se_ <- cmp$gene[cmp$fdr_edger < a]
  data.frame(alpha = a, DESeq2 = length(sd_), edgeR = length(se_),
             shared = length(intersect(sd_, se_)),
             jaccard = round(length(intersect(sd_, se_)) / length(union(sd_, se_)), 3))
}))

The method-specific genes have edgeR FDRs just above 0.05 — they are boundary cases, not contradictions. They also sit at lower baseMean than the shared genes, which is where variance estimation is hardest and the two empirical Bayes schemes diverge most. Agreement is usually highest at strict thresholds (the strongest genes are unambiguous) and degrades as you loosen the cutoff and admit more borderline genes.

4 Block 4 — Pathway analysis

A list of 4,000 significant genes is not a result. Pathway analysis asks whether that list is enriched for genes sharing an annotation.

4.1 Two families of method

	Over-representation (ORA)	Gene set enrichment (GSEA)
Input	a list of significant genes	all genes, ranked
Question	is set S over-represented in my list relative to background?	are set S members concentrated at one end of the ranking?
Test	hypergeometric / Fisher	running enrichment statistic, permutation p-value
Weakness	discards everything below the cutoff; results depend on the cutoff	needs a defensible ranking statistic
Use when	you have a clear, well-populated gene list	effects are subtle but coordinated

4.2 Preparing the inputs

library(clusterProfiler)
library(org.Hs.eg.db)
library(enrichplot)

res_df$lfc_shrunk <- res_shr$log2FoldChange[match(res_df$ensembl, rownames(res_shr))]

sig_up   <- res_df$ensembl[which(res_df$padj < 0.05 & res_df$lfc_shrunk >  1)]
sig_down <- res_df$ensembl[which(res_df$padj < 0.05 & res_df$lfc_shrunk < -1)]

# THE BACKGROUND: every gene that was actually tested. Not "all human genes".
universe <- res_df$ensembl[!is.na(res_df$padj)]

c(up = length(sig_up), down = length(sig_down), universe = length(universe))

##       up     down universe 
##     1530     1312    23001

The single most common mistake in published enrichment analyses is using the whole genome as background when only ~15,000 genes were expressed and tested. Gene sets full of genes that could never have appeared in your foreground then come out “enriched” for purely arithmetic reasons. Always pass universe.

4.3 Over-representation analysis: GO

ego_up <- enrichGO(gene = sig_up, universe = universe, OrgDb = org.Hs.eg.db,
                   keyType = "ENSEMBL", ont = "BP", pAdjustMethod = "BH",
                   pvalueCutoff = 0.05, qvalueCutoff = 0.10, readable = TRUE)

head(as.data.frame(ego_up)[, c("Description", "GeneRatio", "BgRatio", "p.adjust", "Count")], 15)

# GO is a DAG: parent and child terms report nearly the same genes. Collapse them.
ego_up_s <- simplify(ego_up, cutoff = 0.7, by = "p.adjust", select_fun = min)
c(before = nrow(as.data.frame(ego_up)), after = nrow(as.data.frame(ego_up_s)))

## before  after 
##    409    160

dotplot(ego_up_s, showCategory = 15) + labs(title = "GO BP: up in LUSC")

Reading a dot plot: GeneRatio on the x-axis is the fraction of your list in the term, dot size is the absolute count, dot colour is the adjusted p-value. A tiny term with 3 genes and a beautiful p-value is usually less interesting than a broad term with 60.

ego_down <- enrichGO(sig_down, universe = universe, OrgDb = org.Hs.eg.db,
                     keyType = "ENSEMBL", ont = "BP", readable = TRUE)
dotplot(simplify(ego_down), showCategory = 15) + labs(title = "GO BP: up in LUAD")

Always run ORA separately for up and down. Merging them lets opposite effects cancel inside a term and produces bland, uninterpretable results.

fc <- setNames(res_df$lfc_shrunk, res_df$ensembl)
cnetplot(setReadable(ego_up_s, org.Hs.eg.db, keyType = "ENSEMBL"),
         showCategory = 5, foldChange = fc)+
    scale_colour_gradient2(
    low = "blue",
    mid = "white",
    high = "red",
    midpoint = 0
  )

4.4 Over-representation analysis: KEGG

KEGG needs Entrez IDs — a live demonstration of why Day 1’s identifier discipline matters.

to_entrez <- function(ens) {
  x <- AnnotationDbi::mapIds(org.Hs.eg.db, keys = ens, column = "ENTREZID",
                             keytype = "ENSEMBL", multiVals = "first")
  unique(na.omit(unname(x)))
}

kk <- enrichKEGG(gene = to_entrez(sig_up), universe = to_entrez(universe),
                 organism = "hsa", pvalueCutoff = 0.05)
head(as.data.frame(kk)[, c("Description", "GeneRatio", "p.adjust", "Count")], 10)

c(ensembl_in = length(sig_up), entrez_out = length(to_entrez(sig_up)))  # loss in translation

## ensembl_in entrez_out 
##       1530       1223

4.5 GSEA

GSEA needs a ranking of every gene, not a cutoff. The ranking statistic is a real choice:

shrunken log2 fold change — weights effect size (a good default)
sign(LFC) * -log10(p) — weights significance
the Wald statistic res$stat — a compromise between the two

rank_stat <- res_df$lfc_shrunk
names(rank_stat) <- res_df$symbol
rank_stat <- rank_stat[!is.na(rank_stat) & names(rank_stat) != "" &
                         !duplicated(names(rank_stat))]
rank_stat <- sort(rank_stat, decreasing = TRUE)
length(rank_stat); head(rank_stat, 3); tail(rank_stat, 3)

## [1] 23993

##          KPRP RP11-416L21.1         KRT77 
##      9.806862      9.582829      9.511635

##         PAEP RP11-148B3.2          FGB 
##    -6.535603    -6.569795    -6.677779

library(msigdbr)

# msigdbr renamed `category` to `collection` in recent versions — support both.
get_sets <- function(coll, sub = NULL) {
  tryCatch(msigdbr(species = "Homo sapiens", collection = coll, subcollection = sub),
           error = function(e) msigdbr(species = "Homo sapiens", category = coll,
                                       subcategory = sub))
}

hallmark <- get_sets("H")
t2g <- hallmark[, c("gs_name", "gene_symbol")]
length(unique(t2g$gs_name))

## [1] 50

set.seed(42)
gsea_h <- GSEA(rank_stat, TERM2GENE = t2g, pvalueCutoff = 0.25,
               eps = 0, seed = TRUE, verbose = FALSE)
gs <- as.data.frame(gsea_h)
head(gs[order(-abs(gs$NES)), c("Description", "setSize", "NES", "pvalue", "p.adjust")], 15)

gsea_df <- as.data.frame(gsea_h)

gsea_df <- as.data.frame(gsea_h) %>%
  mutate(sign = ifelse(NES > 0, "Activated", "Suppressed"))

ggplot(
  gsea_df,
  aes(
    x = NES,
    y = reorder(Description, NES),
    colour = NES,
    size = -log10(p.adjust)
  )
) +
  geom_point() +
  facet_grid(. ~ sign) +
    scale_colour_gradient2(
    low = "blue",
    mid = "white",
    high = "red",
    midpoint = 0
  ) +
  labs(
    title = "Hallmark gene sets, LUSC vs LUAD",
    x = "Normalized Enrichment Score (NES)",
    y = NULL,
    colour = "NES",
    size = "-log10(adjusted p-value)"
  ) +
  theme_bw()

top_set <- gs$ID[which.max(abs(gs$NES))]
gseaplot2(gsea_h, geneSetID = top_set, title = top_set)

Reading the running-score plot: the curve steps up each time a set member is hit while walking down the ranked list. A peak far to the left, with a dense band of ticks under it, means the set is concentrated among genes up in LUSC. NES is the enrichment score normalised for set size — that is the number you compare across sets, never the raw ES.

ridgeplot(gsea_h, showCategory = 15) +
  labs(
    x = "Gene-level shrunken log2 fold change",
    colour = "NES"
  )

4.6 Pitfalls worth stating out loud

Background set — as above; this alone changes conclusions.
Redundancy — GO returns dozens of near-identical terms. Use simplify(), or report at a fixed GO level.
Gene sets overlap heavily, so the FDR across sets is optimistic. Treat enrichment p-values as a ranking device, not as evidence of the same grade as gene-level FDR.
Set size — sets with < 10 or > 500 genes are rarely informative; both enrichGO and GSEA expose minGSSize/maxGSSize.
Annotation bias — well-studied genes carry more annotations, so “immune response” is enriched in almost everything.
Version everything — MSigDB and GO change between releases; record the version alongside the result.

4.7 Exercise 5 — background and ranking

~15 minutes.

4.7.1 Question

Re-run enrichGO on sig_up without the universe argument. Compare the top 8 terms with and without. What appears, and why?
Re-run GSEA using sign(log2FC) * -log10(pvalue) as the ranking statistic. Do the top hallmark sets change? Merge the two NES tables and correlate them.
Run GSEA on the edgeR result instead of DESeq2 (rank by res_edger$logFC). How much do the enriched sets overlap? What does that tell you about how much the DE-tool choice matters at the pathway level?
Swap Hallmark for C2 CP:REACTOME. Which collection is easier to interpret, and why?

4.7.2 Solution

# 1
ego_nobg <- enrichGO(sig_up, OrgDb = org.Hs.eg.db, keyType = "ENSEMBL",
                     ont = "BP", readable = TRUE)
data.frame(with_bg    = head(as.data.frame(ego_up)$Description, 8),
           without_bg = head(as.data.frame(ego_nobg)$Description, 8))

# 2
r2 <- sign(res_df$lfc_shrunk) * -log10(res_df$pvalue)
names(r2) <- res_df$symbol
r2 <- r2[
  is.finite(r2) &
  !is.na(names(r2)) &
  names(r2) != "" &
  !duplicated(names(r2))
]

r2 <- sort(r2, decreasing = TRUE)
gsea_p <- GSEA(r2, TERM2GENE = t2g, pvalueCutoff = 0.25, eps = 0, verbose = FALSE)
m <- merge(as.data.frame(gsea_h)[, c("ID", "NES")],
           as.data.frame(gsea_p)[, c("ID", "NES")], by = "ID", suffixes = c("_lfc", "_p"))
c(n = nrow(m), r = cor(m$NES_lfc, m$NES_p))

##          n          r 
## 14.0000000  0.9864682

# 3
r3 <- res_edger$logFC; names(r3) <- res_edger$symbol
r3 <- sort(r3[!is.na(r3) & names(r3) != "" & !duplicated(names(r3))], decreasing = TRUE)
gsea_e <- GSEA(r3, TERM2GENE = t2g, pvalueCutoff = 0.25, eps = 0, verbose = FALSE)
sets_d <- as.data.frame(gsea_h)$ID; sets_e <- as.data.frame(gsea_e)$ID
c(DESeq2 = length(sets_d), edgeR = length(sets_e),
  shared = length(intersect(sets_d, sets_e)))

## DESeq2  edgeR shared 
##     15     15     15

# 4
reactome <- get_sets("C2", "CP:REACTOME")
gsea_r <- GSEA(rank_stat, TERM2GENE = reactome[, c("gs_name", "gene_symbol")],
               pvalueCutoff = 0.05, eps = 0, verbose = FALSE)
c(hallmark = nrow(as.data.frame(gsea_h)), reactome = nrow(as.data.frame(gsea_r)))

## hallmark reactome 
##       15       61

Without a background, terms dominated by genes that were never expressed in lung tissue rise to the top. The two ranking statistics give highly correlated NES — GSEA is robust to that choice here. The DESeq2 and edgeR rankings give almost identical pathway results, which is the real lesson of Block 3: the tools disagree about borderline genes, and agree about biology. Hallmark’s 50 broad, deliberately non-redundant sets are far easier to summarise than Reactome’s ~1,600 hierarchical, overlapping ones.

5 Exporting and reporting

out <- data.frame(
  ensembl_gene_id = cmp$gene,
  symbol          = cmp$symbol,
  baseMean        = round(cmp$baseMean, 2),
  log2FC_DESeq2   = round(cmp$lfc_deseq, 3),
  padj_DESeq2     = signif(cmp$padj_deseq, 3),
  log2FC_edgeR    = round(cmp$lfc_edger, 3),
  FDR_edgeR       = signif(cmp$fdr_edger, 3),
  call            = cmp$call
)
out <- out[order(out$padj_DESeq2), ]
write.csv(out, file.path(res_dir, "DE_LUSC_vs_LUAD_DESeq2_edgeR.csv"), row.names = FALSE)
write.csv(as.data.frame(ego_up_s), file.path(res_dir, "GO_BP_up_in_LUSC.csv"), row.names = FALSE)
write.csv(as.data.frame(gsea_h),   file.path(res_dir, "GSEA_hallmark.csv"),    row.names = FALSE)
saveRDS(dds, file.path(data_dir, "dds_LUAD_LUSC.rds"))
head(out, 5)

5.1 A reporting checklist

A reader should be able to answer all of these from your methods alone:

install.packages("renv")
renv::init()      # at the start of a project
renv::snapshot()  # once everything works

Take-home messages

The design formula is the analysis. Check it is full rank and check covariate balance before running anything.
Give raw counts to the model — DESeq2 and edgeR normalise internally, and feeding them CPM or TPM breaks the variance model.
Read the diagnostics. Dispersion plot, MA plot, p-value histogram: each has a characteristic failure mode you can learn to recognise in seconds.
Shrunken fold changes for ranking and plotting; unshrunken p-values for significance.
DESeq2 and edgeR agree almost perfectly on effect sizes and disagree at the significance boundary, mostly for low-expression genes. Choose your primary method in advance and report the other as robustness.
Pathway results are far more stable than gene lists — and far more sensitive to your background set than to your DE tool.

Försöksdjursvetenskap (teoretisk) och avancerad databearbetning

DESeq2, edgeR, method comparison, and pathway analysis

Birong Zhang

07 September 2026

How today is organised

0.1 Learning objectives

0.2 Packages

0.3 Loading yesterday’s object

1 Block 1 — Design matrices and DESeq2

1.1 What the formula actually means

1.2 Filtering

1.3 Building the DESeq2 object

1.4 Step 1 — size factors (median of ratios)

1.5 Step 2 — dispersion

1.6 Step 3 — testing

1.7 Independent filtering and padj = NA

1.8 Log fold change shrinkage

1.9 Multiple testing

1.10 Testing against a fold-change threshold

1.11 Visualising the DESeq2 result

1.12 Exercise 1 — does the covariate matter?

1.12.1 Question

1.12.2 Solution

1.13 Exercise 2 — reading the diagnostics

1.13.1 Question

1.13.2 Solution

2 Block 2 — edgeR

2.1 The object and TMM normalisation

2.2 Dispersion

2.3 Quasi-likelihood test

2.4 Fold-change-thresholded testing

2.5 edgeR’s own MD plot

2.6 Exercise 3 — normalisation and test choice

2.6.1 Question

2.6.2 Solution

3 Block 3 — Comparing DESeq2 and edgeR

3.1 Correlation of effect sizes

3.2 Overlap of the significant sets

3.3 Who are the disagreements?

3.4 Concordance at the top

3.5 What to do with two lists

3.6 Exercise 4 — characterise the disagreement

3.6.1 Question

3.6.2 Solution

4 Block 4 — Pathway analysis

4.1 Two families of method

4.2 Preparing the inputs

4.3 Over-representation analysis: GO

4.4 Over-representation analysis: KEGG

4.5 GSEA

4.6 Pitfalls worth stating out loud

4.7 Exercise 5 — background and ranking

4.7.1 Question

4.7.2 Solution

5 Exporting and reporting

5.1 A reporting checklist

Take-home messages

Further reading

1.7 Independent filtering and `padj = NA`