9. Differential Gene Expression Analysis with DESeq2 in R

RNA-seq read counts exhibit overdispersion—the variance across biological replicates exceeds the mean ($\sigma^2 > \mu$). Therefore, standard Poisson models underestimate false discovery rates. DESeq2 models count data using the Negative Binomial Generalized Linear Model (GLM) with Empirical Bayes shrinkage:

Production R Script: Complete DESeq2 Pipeline

library(DESeq2)
library(tidyverse)

# 1. Load Count Matrix and Sample Metadata Table
counts <- read.table("05_counts/counts_matrix_clean.tsv", header = TRUE, row.names = 1, sep = "\t")
colData <- read.csv("sample_metadata.csv", row.names = 1)

# Ensure sample names match exactly
all(colnames(counts) == rownames(colData))

# Set reference baseline condition (e.g. "Control")
colData$Condition <- factor(colData$Condition, levels = c("Control", "Treated"))

# 2. Construct DESeqDataSet
dds <- DESeqDataSetFromMatrix(
  countData = counts,
  colData = colData,
  design = ~ Batch + Condition   # Control for batch effects + main treatment
)

# 3. Filter ultra-low count genes (<10 reads across samples)
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep, ]

# 4. Execute Differential Expression Analysis
# Performs: Size factor estimation -> Dispersion estimation -> GLM Wald testing
dds <- DESeq(dds)

# 5. Extract Results with apeglm Log2 Fold Change Shrinkage
res <- results(dds, contrast = c("Condition", "Treated", "Control"), alpha = 0.05)
res_shrink <- lfcShrink(dds, coef = "Condition_Treated_vs_Control", type = "apeglm")

# Export Full Results to CSV
write.csv(as.data.frame(res_shrink), file = "06_deseq2/DEG_Treated_vs_Control_apeglm.csv")

# Filter Significant DEGs: FDR < 0.05 & |log2FC| >= 1
sig_degs <- subset(res_shrink, padj < 0.05 & abs(log2FoldChange) >= 1)
message("Total Significant DEGs: ", nrow(sig_degs))