Chapter 8: Differential Gene Expression Analysis with DESeq2: Negative Binomial GLMs

1. The Theoretical Architecture of DESeq2

Simple statistical tests like the Student's \(t\)-test or standard Poisson regression fail when applied to raw RNA-Seq counts because sequencing data exhibits extreme overdispersion: the variance among biological replicates is significantly larger than the mean. DESeq2 models raw counts \(K_{ij}\) using a Negative Binomial (Gamma-Poisson) distribution with mean \(\mu_{ij}\) and gene-specific dispersion parameter \(\alpha_i\):

$$K_{ij} \sim \text{NB}(\mu_{ij}, \alpha_i)$$

$$\text{Var}(K_{ij}) = \mu_{ij} + \alpha_i \mu_{ij}^2$$

Empirical Bayes Dispersion Shrinkage

When an experiment has only 3–5 replicates per condition, individual gene-wise estimates of variance are highly unreliable. DESeq2 uses an Empirical Bayes shrinkage procedure: it fits a global mean-dispersion trend across all 20,000 genes, and then shrinks each gene's individual dispersion estimate toward the trend line. Genes with unusually high variance (true outliers) are preserved, while noisy low-count genes have their dispersions stabilized.

2. Multi-Factor Experimental Design Formulas

DESeq2 utilizes standard R design formulas (~ terms) to model complex biological and technical variables simultaneously:

3. Complete End-to-End DESeq2 Production Workflow

library(DESeq2)
library(apeglm) # Recommended modern empirical Bayes LFC shrinkage

# 1. Build DESeqDataSet from count matrix and metadata
dds <- DESeqDataSetFromMatrix(
    countData = mock_counts,
    colData = clinical_clean,
    design = ~ Extraction_Batch + Condition
)

# 2. Pre-filter low count rows
keep <- rowSums(counts(dds) >= 10) >= 3
dds <- dds[keep, ]

# 3. Execute DESeq Pipeline
# (Estimates size factors, estimates dispersions, fits GLM models, Wald tests)
dds <- DESeq(dds)

# 4. Check Dispersion Shrinkage Diagnostic Plot
# Points should shrink toward the red trend line
plotDispEsts(dds, main = "DESeq2 Dispersion Estimates & Empirical Bayes Shrinkage")

# 5. Extract Results with Empirical Bayes LFC Shrinkage (apeglm)
# Lists available coefficient names:
resultsNames(dds)

# Apply apeglm shrinkage to eliminate noisy large fold changes in low-count genes
res_shrunk <- lfcShrink(
    dds,
    coef = "Condition_Treated_vs_Control",
    type = "apeglm"
)

# Convert to tidy dataframe
res_df <- as.data.frame(res_shrunk)
res_df$Gene_Symbol <- rownames(res_df)

# Annotate significance call
res_df$Call <- case_when(
    res_df$padj < 0.05 & res_df$log2FoldChange >= 1.0  ~ "UP",
    res_df$padj < 0.05 & res_df$log2FoldChange <= -1.0 ~ "DOWN",
    TRUE ~ "NS"
)

# Top 10 Differentially Expressed Genes
top_degs <- res_df %>% filter(Call != "NS") %>% arrange(padj) %>% head(10)
cat("=== Top 10 High-Confidence Differentially Expressed Genes ===\n")
print(top_degs[, c("Gene_Symbol", "baseMean", "log2FoldChange", "pvalue", "padj", "Call")])