/deseq2-differential-expression

name: "deseq2-differential-expression"
description: "Bulk RNA-seq DE with R/Bioconductor DESeq2. Negative binomial GLM, empirical Bayes shrinkage, Wald/LRT tests, multi-factor designs, Salmon tximeta import, apeglm LFC shrinkage, MA/volcano/heatmap viz. R gold standard. Use pydeseq2-differential-expression for Python; use edgeR for TMM normalization."
license: "LGPL-3.0"

DESeq2 Differential Expression Analysis (R/Bioconductor)

Overview

DESeq2 is the Bioconductor R package for differential gene expression analysis from bulk RNA-seq count data. It fits a negative binomial generalized linear model per gene, estimates dispersion parameters using empirical Bayes shrinkage across genes, and tests differential expression using Wald tests (two-group) or likelihood ratio tests (complex designs). DESeq2 is the R gold standard for RNA-seq DE analysis, with native Bioconductor integration for seamless import from Salmon (tximeta/tximport), featureCounts, or HTSeq.

When to Use

• Identifying differentially expressed genes between two experimental conditions (treated vs. control, disease vs. healthy) from bulk RNA-seq count data
• Analyzing multi-factor designs that account for batch effects or covariates (e.g., `~ batch + condition`)
• Testing complex hypotheses with interaction terms (e.g., time × treatment) or reduced models using likelihood ratio tests (LRT)
• Importing Salmon pseudoalignment output via tximeta or tximport for transcript-level uncertainty propagation
• Performing LFC shrinkage with apeglm for ranked gene lists, volcano plots, and downstream pathway analysis
• Conducting time-series experiments or any design with more than two levels requiring model comparison
• Working in an R/Bioconductor ecosystem where integration with SummarizedExperiment, clusterProfiler, or EnhancedVolcano is needed
• Use **pydeseq2-differential-expression** instead for Python-based pipelines with the same statistical model
• Use **edgeR** for negative binomial DE with TMM normalization, quasi-likelihood F-tests, or TREAT testing
• Use **omics-plotting** SKILL after DE and for publication-quality plots of DESeq2 results
• Use **gseapy-gene-enrichment** SKILL after DE to interpret results at the pathway level

Prerequisites

• **R packages**: `DESeq2` (Bioconductor), `tximeta` or `tximport` (Salmon import), `apeglm` (LFC shrinkage), `pheatmap`, `ggplot2`, `EnhancedVolcano`
• **Data requirements**: Raw (unnormalized) integer count matrix — gene rows × sample columns — plus a sample metadata data frame with matching column names. If using Salmon: per-sample `quant.sf` files and a transcript-to-gene mapping
• **Environment**: R ≥ 4.2, Bioconductor ≥ 3.16

# Install Bioconductor packages
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install(c("DESeq2", "tximeta", "tximport", "apeglm",
                        "EnhancedVolcano"))
install.packages(c("pheatmap", "ggplot2", "dplyr"))

Pre-flight Interview

Settle these with the user before writing any analysis code.

decisions:
  - id: D1
    param: design
    kind: required
    source: data
    ask: "Which column of the sample sheet separates the groups you want to compare?"
    default: null

  - id: D2
    param: contrast
    kind: required
    source: user
    depends_on: [D1]
    ask: "Which two groups should be compared, and which group should have positive log2 fold changes?"
    default: null

  - id: D3
    param: design
    kind: required
    source: data
    depends_on: [D1]
    ask: "Which available variables should be included to adjust for nuisance variation (batch, donor, sex, sequencing run)?"
    default: "none"
    skip_if: "sample sheet carries no column besides the grouping variable"

  - id: D4
    param: test
    kind: required
    source: user
    depends_on: [D1, D2]
    ask: "Is the primary hypothesis a named comparison, or an overall effect of this factor across its levels?"
    default: "named comparison (Wald test)"

  - id: D5
    param: alpha
    kind: optional
    source: user
    ask: "How strong must the evidence be before a gene counts as changed (false-discovery rate)?"
    default: 0.05

  - id: D6
    param: lfcThreshold
    kind: optional
    source: user
    ask: "Test against a zero change, or against a minimum fold change worth caring about?"
    default: 0

  - id: D7
    param: lfcShrink.type
    kind: derived
    source: upstream
    ask: "Shrink fold-change estimates so low-count genes stop dominating the rankings and plots?"
    default: "apeglm, applied for ranking and visualization"

  - id: D8
    param: blind (vst/rlog)
    kind: derived
    source: upstream
    ask: "Should the variance-stabilizing transform ignore the design?"
    default: "TRUE for QC/clustering, FALSE for design-aware downstream uses"

  - id: D9
    param: independentFilter
    kind: optional
    source: user
    ask: "Should DESeq2 use its adaptive mean-count filter to maximize discoveries at the chosen FDR threshold?"
    default: true

`design` appears twice because the formula is assembled from both D1 and D3; D2 defines the direction of the reported effect. The pipeline explicitly passes `alpha: 0.05` to `results()` so independent filtering is optimized for the reported FDR threshold rather than DESeq2's default of `0.1`.

Quick Start

Complete two-group comparison from a count matrix in under 20 lines.

library(DESeq2)

# Load count matrix (genes x samples) and metadata
counts <- as.matrix(read.csv("counts.csv", row.names = 1))
coldata <- read.csv("metadata.csv", row.names = 1)
coldata$condition <- factor(coldata$condition)

# Build DESeqDataSet, run full pipeline
dds <- DESeqDataSetFromMatrix(countData = counts,
                               colData   = coldata,
                               design    = ~ condition)
dds <- dds[rowSums(counts(dds)) >= 10, ]   # pre-filter
dds <- DESeq(dds)

# Extract results (treated vs. control)
res <-