DESeq2 & EdgeR RNA-Seq Differential Expression Workbench

📈 Robust Statistical Modeling for RNA-Seq Counts

RNA-seq count data exhibits biological variance exceeding Poisson expectations (overdispersion). DESeq2 models gene counts with a Negative Binomial generalized linear model (GLM), using empirical Bayes shrinkage for dispersions and fold-changes to reliably detect differential expression even with small replicate sizes ($N \ge 3$).

🌋 Interactive Volcano Plot & Significance Simulator

Interactive Visualizer

Log2 Fold Change Cutoff (|log2FC| ≥):

0.5 (1.4×) 1.0 (2.0×) 3.0 (8.0×)

Adjusted p-value / FDR Threshold (padj ≤):

🔴 Down-regulated: 128 genes

🟢 Up-regulated: 164 genes

⚪ Not Significant: 4,708 genes

📐 Experimental Design Formula & Contrast Generator

Interactive Tool

Experimental Structure

Experimental Design Type:

Reference Baseline Level (Control):

Generated R Formula & Contrast

design = ~ condition

results(dds, contrast = c("condition", "Treated", "Control"))

0. The Negative Binomial Model ($Y \sim \text{NB}(\mu, \alpha)$)

Because biological replicates exhibit variance exceeding Poisson sampling noise ($\text{Var} = \mu + \alpha \mu^2$), DESeq2 models raw unnormalized read counts with a Negative Binomial distribution where $\alpha$ represents dispersion. Information is pooled across genes to shrink dispersion estimates toward a common trend.

1. Count Matrix & Sample Metadata Ingestion

Input raw non-normalized integer counts and a matching sample annotation data frame.

library(DESeq2)
library(tidyverse)

# Read raw counts matrix (rows = genes, cols = sample IDs)
counts_matrix <- read.table("gene_counts.txt", header=TRUE, row.names=1, sep="\t")

# Read metadata (row names must exactly match column names of counts_matrix)
col_data <- read.csv("sample_metadata.csv", row.names=1)
stopifnot(all(colnames(counts_matrix) == rownames(col_data)))

2. Setting Factor Reference Levels (relevel)

Explicitly designate the control condition as the reference level so fold-changes are positive for up-regulated genes.

# Explicitly designate 'Control' as denominator
col_data$condition <- relevel(factor(col_data$condition), ref = "Control")

# Build the DESeqDataSet object
dds <- DESeqDataSetFromMatrix(countData = counts_matrix,
                              colData = col_data,
                              design = ~ condition)

3. Pre-filtering Low-Count Noise Genes

Remove genes with zero or negligible read counts across samples to reduce multiple testing burden.

# Keep only genes with at least 10 reads in at least 3 samples (min group size)
keep <- rowSums(counts(dds) >= 10) >= 3
dds <- dds[keep, ]

4. Running Core Differential Expression (DESeq)

Executes size factor estimation, dispersion estimation, and negative binomial Wald tests in one function.

# Run DESeq2 statistical workflow
dds <- DESeq(dds)

# Inspect estimated dispersion curve
plotDispEsts(dds)

5. Variance-Stabilizing Transformation (VST / rlog)

Apply VST to normalize count heteroskedasticity for distance clustering and exploratory PCA plots.

# Variance-stabilizing transformation for N > 30 samples (or rlog for small cohorts)
vsd <- vst(dds, blind = FALSE)

6. Exploratory PCA & Sample Distance Heatmaps

Verify that biological replicates cluster tightly together and identify potential outlier samples or batch effects.

library(pheatmap)
library(RColorBrewer)

# Principal Component Analysis
plotPCA(vsd, intgroup = c("condition"))

# Sample Euclidean distance matrix
sampleDists <- dist(t(assay(vsd)))
sampleDistMatrix <- as.matrix(sampleDists)
pheatmap(sampleDistMatrix,
         clustering_distance_rows = sampleDists,
         clustering_distance_cols = sampleDists,
         col = colorRampPalette(rev(brewer.pal(9, "Blues")))(255))

7. Extracting Differential Results & FDR Correction

Extract log2 fold changes, Wald test statistics, and Benjamini-Hochberg adjusted p-values.

# Extract results table at 5% FDR
res <- results(dds, contrast=c("condition", "Treated", "Control"), alpha=0.05)
summary(res)

# Export significant genes table
res_sig <- subset(res, padj < 0.05 & abs(log2FoldChange) >= 1)
write.csv(as.data.frame(res_sig), "DESeq2_significant_genes.csv")

8. Log2 Fold Change Shrinkage (apeglm / ashr)

Shrink noisy fold changes from low-count genes without modifying high-confidence signal.

# Shrink LFC with apeglm algorithm
res_shrunk <- lfcShrink(dds, coef = "condition_Treated_vs_Control", type = "apeglm")
plotMA(res_shrunk, ylim = c(-5, 5))

9. Enhanced Volcano Plotting

Generate publication-quality volcano plots highlighting top candidates.

library(EnhancedVolcano)

EnhancedVolcano(res_shrunk,
  lab = rownames(res_shrunk),
  x = 'log2FoldChange',
  y = 'padj',
  pCutoff = 0.05,
  FCcutoff = 1.0,
  pointSize = 2.0,
  labSize = 4.0,
  title = 'Treated vs Control Volcano Plot')

10. Pathway & Gene Ontology Enrichment (clusterProfiler)

Identify overrepresented biological processes (GO-BP) and KEGG metabolic pathways in up/downregulated gene sets.

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

# Extract significant upregulated Entrez IDs
up_genes <- rownames(subset(res_sig, log2FoldChange > 0))

# Over-representation analysis (GO Biological Process)
ego <- enrichGO(gene          = up_genes,
                OrgDb         = org.Hs.eg.db,
                keyType       = 'SYMBOL',
                ont           = "BP",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.05)

dotplot(ego, showCategory=15, title="Top Enriched Biological Processes")

🔀 Decision Architecture: DESeq2 vs EdgeR vs Limma-voom

Decision Tree

Select Replicate Cohort Size:

Recommended: DESeq2 or EdgeR (Quasi-Likelihood glmQLFit)
For typical biological experiments with 3-10 replicates, DESeq2 provides the most reliable error control through empirical Bayes dispersion shrinkage and apeglm fold-change moderation.