Bio Differential Expression Deseq2 Basics

Version Compatibility

Reference examples tested with: DESeq2 1.42+, apeglm 1.28+, ashr 2.2+, IHW 1.34+, tximport 1.30+, edgeR 4.0+ (for cross-comparison), PyDESeq2 0.5+

Before using code patterns, verify installed versions match. If versions differ:

• R: packageVersion('<pkg>') then ?function_name to verify parameters
• Python: pip show <package> then help(module.function) to check signatures

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

DESeq2 Basics

"Find differentially expressed genes between conditions" -> Fit a negative-binomial GLM per gene with shared dispersion shrinkage, test the coefficient of interest (Wald) or the joint effect of a factor (LRT), and report a shrunken effect-size estimate for ranking.

The Single Most Important Modern Insight -- Shrunken LFC and the Wald p-value come from different models

lfcShrink() returns LFCs from a Bayesian posterior with apeglm/ashr/normal priors, BUT the p-value column it carries forward is still the unshrunken Wald p-value from results(). This is a deliberate design choice (Zhu, Ibrahim, Love 2019 Bioinformatics 35:2084) -- the shrunken estimate is for ranking and visualization; the p-value is for inference. Reporting "shrunken LFC = 0.4, padj = 1e-8" mixes two models, which is fine because both are correct for their stated purpose. What is NOT fine: using the shrunken LFC in a downstream filter and then claiming FDR control on that filter (it has none). For threshold-based FDR claims, use lfcThreshold= or TREAT (glmTreat in edgeR).

A second consequence: results(dds) with no name= or contrast= argument silently returns the last coefficient in resultsNames(dds) -- which depends on factor level order and design formula order. Always specify the contrast explicitly. Tutorials that hard-code results(dds) are setting an example that breaks the moment another factor is added.

Algorithmic Taxonomy

Decision Tree by Scenario

Standard Workflow

Goal: Take a raw integer count matrix and a sample table to a ranked, shrunken DE result table.

Approach: Construct DESeqDataSet with the design formula, set reference levels explicitly, run the pipeline, extract by explicit contrast, shrink for downstream use.

library(DESeq2)
library(apeglm)

dds <- DESeqDataSetFromMatrix(countData = counts, colData = coldata, design = ~ condition)
dds$condition <- relevel(dds$condition, ref = 'control')

keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep, ]

dds <- DESeq(dds)
resultsNames(dds)

res <- results(dds, name = 'condition_treated_vs_control', alpha = 0.05)
res_shrunk <- lfcShrink(dds, coef = 'condition_treated_vs_control', type = 'apeglm')

summary(res)
sig <- subset(res, padj < 0.05)

The reference level fix is non-cosmetic: DESeq2 picks alphabetically if not told otherwise, so c('Treated','Untreated') makes 'Treated' the reference and the LFC reads inverted. Set it BEFORE DESeq().

Tximport (Salmon / kallisto / RSEM)

For salmon/kallisto/RSEM input, use DESeqDataSetFromTximport() which carries the per-sample length matrix as an offset automatically:

library(tximport)
txi <- tximport(files, type = 'salmon', tx2gene = tx2gene)
dds <- DESeqDataSetFromTximport(txi, colData = samples, design = ~ condition)
dds <- DESeq(dds)

The tximport(..., countsFromAbundance='lengthScaledTPM') form is for limma-voom (no offset mechanism). Full mechanics, the four countsFromAbundance options, RSEM zero-length traps, and tximeta provenance: see expression-matrix/counts-ingest.

Design Formulas and the resultsNames Trap

Goal: Encode batch, paired, and interaction structure correctly and extract the intended contrast.

Approach: Put the variable of interest LAST for readability, but never trust the default results(dds) -- inspect resultsNames(dds) and pass name= or contrast= explicitly.

design(dds) <- ~ batch + condition
dds <- DESeq(dds)
resultsNames(dds)
# "Intercept" "batch_B_vs_A" "condition_treated_vs_control"

res <- results(dds, name = 'condition_treated_vs_control')

Interaction design with the canonical trap:

design(dds) <- ~ genotype + treatment + genotype:treatment
dds <- DESeq(dds)
resultsNames(dds)
# "Intercept" "genotype_KO_vs_WT" "treatment_drug_vs_vehicle" "genotypeKO.treatmentdrug"

• results(dds, name='treatment_drug_vs_vehicle') returns the drug effect IN THE WT REFERENCE only, NOT a marginal average. This is the single most common misinterpretation.
• results(dds, name='genotypeKO.treatmentdrug') returns the DIFFERENCE in drug effect between KO and WT.
• Drug effect in KO requires summing: results(dds, contrast=list(c('treatment_drug_vs_vehicle','genotypeKO.treatmentdrug'))).

Cleaner alternative for interactions when many pairwise contrasts are needed: combined factor + ~ 0 + group.

dds$group <- factor(paste(dds$genotype, dds$treatment, sep = '_'))
design(dds) <- ~ 0 + group
dds <- DESeq(dds)
res_drug_in_ko <- results(dds, contrast = c('group', 'KO_drug', 'KO_vehicle'))

Wald vs LRT

Goal: Choose Wald for single-coefficient hypotheses, LRT for joint hypotheses involving more than 1 df.

Approach: Wald is default. LRT is mandatory for multi-level factors tested as "any change", interactions with >1 df, and ANOVA-style omnibus tests. With LRT, the reported LFC is for the LAST coefficient in resultsNames(dds) -- not the omnibus effect.

dds <- DESeq(dds, test = 'LRT', reduced = ~ batch)
res_lrt <- results(dds)
# res_lrt$padj is the LRT joint p-value (correct).
# res_lrt$log2FoldChange is for the last coefficient in resultsNames -- NOT an omnibus summary.

When reporting LRT results, name what the LFC actually represents or extract specific Wald coefficients per level for the effect-size table.

LFC Shrinkage -- Three Flavors, Three Failure Modes

Goal: Get a stable effect-size estimate appropriate for ranking, GSEA input, volcano-plot x-axis, and reporting.

Approach: Default to apeglm. Switch to ashr when arbitrary contrasts are needed. Never use normal for new analyses.

res_apeglm <- lfcShrink(dds, coef = 'condition_treated_vs_control', type = 'apeglm')
res_ashr   <- lfcShrink(dds, contrast = c('condition','treated','control'), type = 'ashr')

The apeglm-cannot-use-contrast footgun: if the question is "drug effect in KO" from ~ genotype * treatment, apeglm cannot directly shrink that contrast. Workarounds: (a) rebuild as combined factor ~ 0 + group and relevel so the desired comparison is a coefficient; (b) use ashr; (c) accept the unshrunken LFC for that one comparison.

p-values do NOT change when shrinking. lfcShrink() preserves the Wald p-value from results(). See the Single Most Important Insight at the top.

Independent Filtering, Cook's Outliers, padj=NA

padj = NA has three distinct causes (independent filtering, Cook's outlier, all-zero in a group), each with a different remediation -- see de-results for the full diagnostic table, IHW alternative, and recovery code.

Two DESeq2-specific points worth knowing at this layer:

• Cook's distance filtering is NOT computed for continuous covariates -- a continuous-covariate analysis has effectively no automatic outlier filtering. Disable Cook's only when the outlier IS the signal (results(dds, cooksCutoff = FALSE)). At n>=7 per group, DESeq() REPLACES outliers via replaceOutliers() and refits (minReplicatesForReplace = 7).
• Pre-filtering (rowSums(counts(dds)) >= 10) is for memory and speed ONLY. It does NOT replace independent filtering, which operates downstream at results() time. Independent filtering can be swapped for IHW via results(dds, filterFun = ihw).

VST vs rlog vs normTransform (Visualization Only)

Goal: Produce homoskedastic log2-scale counts for PCA, heatmaps, clustering, ML features.

Approach: Use vst() by default. Switch to rlog() only for n<30 with size factors varying >4x. Never use normTransform() for distance-based plots. Never use VST/rlog values as input to DE.

vsd <- vst(dds, blind = FALSE)
rld <- rlog(dds, blind = FALSE)

The vst() function default is blind=TRUE, but the current DESeq2 vignette recommends blind=FALSE for any downstream visualization AFTER the model is fit (it uses the design when fitting dispersions; appropriate when the design is already settled). Reserve blind=TRUE for unsupervised QC where the design should not influence the transformation (e.g., "is this sample consistent with its group?"). The vignette has flip-flopped over the years on which to recommend by default -- pass blind= explicitly.

vst() uses 1000 most-variable genes to fit the dispersion trend by default. With <1000 genes after filtering, set nsub lower.

betaPrior Deprecation Timeline

DESeq2 v1.0-v1.15: betaPrior = TRUE was the default. Shrinkage was baked in and p-values were computed on the shrunken estimates. v1.16 (2017) flipped the default to FALSE and introduced lfcShrink(). Today: betaPrior = TRUE is quasi-deprecated and is the ONLY path to a p-value of the shrunken estimate. Most users do not want this. lfcShrink() with apeglm and the unshrunken Wald p-value is the modern compromise.

The contrast= vs name= numerical difference that older tutorials describe was specific to betaPrior = TRUE. With the current default, name= and contrast= for the same comparison return identical LFCs.

Size Factor Alternatives

estimateSizeFactors() defaults to type='ratio' (median-of-ratios). Edge cases:

Single-cell pseudobulk: most genes have zeros across donors, so type='poscounts' is often required for pseudobulk DESeq2.

Per-Method Failure Modes