bio-differential-expression-deseq2-basics

name	bio-differential-expression-deseq2-basics
description	Performs differential expression on bulk RNA-seq count data with DESeq2's negative-binomial GLM, Wald and LRT testing, apeglm/ashr/normal LFC shrinkage, independent filtering, Cook's outlier handling, VST/rlog transforms, and design formulas including paired, batch, and interaction terms. Use when running bulk DE, choosing DESeq2 over edgeR or limma-voom, building a paired or interaction design, applying LFC shrinkage for ranking or GSEA, choosing Wald vs LRT, troubleshooting padj=NA, picking VST vs rlog, importing salmon/kallisto via tximport, or analyzing prokaryotic RNA-seq.
tool_type	r
primary_tool	DESeq2

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

Test / estimator	What it tests	When mandatory	Failure mode
Wald	One coefficient = 0	Two-level factor or single contrast	Anti-conservative with many low-count outliers
LRT (test='LRT', reduced=)	Joint effect of dropped terms (>=1 df)	Multi-level factor, omnibus, interaction with >1 df	Reports LFC of the LAST coefficient, not omnibus -- read p-value but never report the LFC as "the effect"
lfcShrink(type='apeglm') (Zhu 2019)	Posterior LFC under heavy-tailed Cauchy prior	DEFAULT for ranking and visualization	Requires coef=; cannot use contrast= or numeric vectors
lfcShrink(type='ashr') (Stephens 2017)	Posterior LFC under unimodal prior; reports lfsr/svalue with svalue=TRUE	Arbitrary contrasts via contrast=	Slightly different inferential frame (sign-error rather than null-FDR)
lfcShrink(type='normal')	Posterior LFC under zero-centered normal; accepts coef= or contrast= (with res=)	Quasi-deprecated since v1.16; only path to get shrunken p-values	Cannot be used with formulas containing interaction terms
TREAT / lfcThreshold= (McCarthy & Smyth 2009)	LFC magnitude exceeds threshold tau	Want FDR control for "	LFC

Decision Tree by Scenario

Scenario	Recommended path	Why
Two-group bulk RNA-seq, n>=3/group	DESeq() + results(name=...) + lfcShrink(coef=..., type='apeglm')	Modern default; apeglm is the right prior
Factor with 3+ levels, "any change" question	DESeq(test='LRT', reduced=~1); read padj only	Wald + p-value combining is wrong
Interaction ~ genotype * treatment	Build combined factor group = paste(genotype, treatment) and design ~ 0 + group; contrast pairs of interest	Avoids the resultsNames trap; works with apeglm via relevel
Paired design (tumor/normal same patient)	~ patient + tissue; pairing variable FIRST	Absorbs subject variability; n_paired effective sample size
Salmon/kallisto input	tximport() -> DESeqDataSetFromTximport()	Carries length offsets automatically; DESeqDataSetFromMatrix(round(...)) loses length correction
n=2/group, no choice	Continue but report results as exploratory; consider edgeR QL F-test as sensitivity	Schurch 2016 RNA 22:839: all tools miss 20-40% of true positives at n=3
Single-cell pseudobulk (counts aggregated per donor)	DESeq2 standard pipeline on pseudobulk matrix	Crowell 2020 Nat Commun 11:6077: pseudobulk avoids the FDR inflation of cell-level DE
Many DE genes expected (>50% of genome)	estimateSizeFactors(controlGenes=stable) or spike-in normalization	Median-of-ratios assumes most genes unchanged
GSEA preranked input	Shrunken LFC OR stat (Wald Z) as the rank	Unshrunken LFC dominated by low-count noise
Cross-sample heatmap, PCA, ML feature	vst(dds, blind=FALSE) (or rlog if n<30 and library sizes vary >4x)	Raw counts make PC1 = library size

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')

summaryres
sig  subsetres padj  

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')

Method	Prior	Accepts	Use when
apeglm	Cauchy (heavy-tailed)	coef= only	Default; preserves large effects, suppresses low-count noise
ashr	Unimodal scale-mixture	coef= or contrast= (incl. numeric)	Need contrast= for interaction sums or pairwise from ~ 0 + group
normal	Zero-centered normal	coef= only; not interaction designs	Only when the old shrunken p-value is required (legacy)

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:

Situation	Use
Zero counts in some samples for many genes	type='poscounts' -- uses only positive entries per gene
Very small libraries, hard to converge	type='iterate'
Spike-ins (ERCC) or known stable housekeeping genes	controlGenes = indices
Majority-DE biology (prokaryotic stress, viral host shutoff, MYC amplification)	controlGenes with curated stable genes; or spike-in SBN (Jiang 2011 Genome Res 21:1543)

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

Per-Method Failure Modes

apeglm refuses arbitrary contrasts

Trigger: Question of the form "drug effect in KO genotype" from ~ genotype * treatment; user tries lfcShrink(dds, contrast=list(...), type='apeglm').

Mechanism: apeglm fits per-coefficient priors. A numeric or list contrast is a linear combination of coefficients, not a coefficient -- no prior to apply.

Symptom: Error: "type='apeglm' shrinkage only for use with 'coef'"

Fix: Rebuild design as ~ 0 + group with group = paste(genotype, treatment), relevel so the comparison is a coefficient, refit. Or use type='ashr' which accepts contrasts.

LRT reports the wrong LFC

Trigger: Multi-level factor analyzed with DESeq(dds, test='LRT', reduced=~1); user reports the log2FoldChange column as "the effect".

Mechanism: The LRT p-value is the omnibus test of "any difference among levels". The LFC reported by results() after LRT is for the LAST coefficient in resultsNames(dds), which is one specific level-vs-reference comparison.

Symptom: A 4-level factor produces a single LFC value per gene; reviewer asks "the effect of which condition?"

Fix: Treat LRT padj as a screen for "any change". For effect sizes, extract per-level Wald coefficients individually via results(dds, name='<specific coefficient>') for each non-reference level.

Cook's silences the gene of interest

Trigger: Rare-disease cohort or CNV-amplified patient where ONE sample drives the biology; that gene has padj=NA.

Mechanism: Cook's distance flagged the patient as an outlier and zeroed the gene's p-value (or replaced its counts if n>=7).

Symptom: A gene of clear biological interest comes back as NA in the results table even though the count matrix shows the expected pattern.

Fix: results(dds, cooksCutoff = FALSE). Optionally cross-validate the result by running a sensitivity analysis with and without the outlier sample.

Independent filtering kills the master regulator

Trigger: A transcription factor expressed at ~10 counts but consistently across all samples shows padj=NA despite obvious biology.

Mechanism: baseMean is below the data-driven independent filtering threshold; the gene was excluded from FDR adjustment.

Symptom: Low-count genes with clean signal end up NA.

Fix: results(dds, independentFiltering = FALSE); or filterFun = ihw from the IHW package (often less aggressive on low-count genes).

Parametric dispersion-mean trend doesn't fit the cloud

Trigger: plotDispEsts(dds) shows the red parametric trend curve nowhere near the cloud of gene-wise (blue) and final (black) dispersion estimates.

Mechanism: Default fitType='parametric' assumes dispersion ~ a/mean + b. Fails when the experiment has very few samples per group with highly heterogeneous biology, many very-low-count genes pulling the trend, or a continuous covariate driving large variability.

Symptom: Curved trend that doesn't match the cloud; mismatch shows up most clearly in low-baseMean genes.

Fix: Refit with DESeq(dds, fitType='local') (local regression) or fitType='mean' (flat trend); compare plotDispEsts between fits and pick the one tracking the cloud. Falling back to 'mean' is a sign the data is unusual; investigate before trusting results.

Median-of-ratios fails on prokaryotic stress

Trigger: Bacterial RNA-seq under stress where >50% of genes change in one direction; PCA shows huge global shift.

Mechanism: Median-of-ratios assumes most genes are not DE. Under massive global perturbation, the reference is dominated by DE genes; size factors absorb the biology.

Symptom: MA plot shows the bulk cloud shifted off zero; reported fold changes don't match qPCR.

Fix: Use controlGenes with curated stable housekeeping genes; or spike-in normalization (Jiang 2011); or RUVg with negative controls.

PyDESeq2 (Python alternative)

from pydeseq2.dds import DeseqDataSet
from pydeseq2.ds import DeseqStats

dds = DeseqDataSet(counts=count_df, metadata=metadata, design='~condition')
dds.deseq2()
stat_res = DeseqStats(dds, contrast=('condition', 'treated', 'control'))
stat_res.summary()
results_df = stat_res.results_df

PyDESeq2 0.5+ supports Wald, multi-factor designs, and apeglm shrinkage. No LRT yet. Results are numerically close to R DESeq2 but with small differences from MLE solver choice.

Prokaryotic RNA-seq

• Non-spliced aligners (BWA-MEM, Bowtie2) -- no introns.
• Polycistronic operons cause read-through between adjacent genes; confirm gene boundaries in the GFF.
• rRNA depletion essential (80-95% rRNA without poly-A selection, which prokaryotes lack anyway).
• Median-of-ratios fails under stress (see failure mode above) -- use controlGenes or spike-ins.
• KEGG organism codes are strain-specific (e.g., pae for P. aeruginosa PAO1): clusterProfiler::search_kegg_organism().
• Annotation comes from Prokka or Bakta GFF; Ensembl/biomaRt are eukaryote-only.

Common errors

Error / symptom	Cause	Fix
design matrix not full rank	Confounded covariates	alias(model.matrix(design, coldata))$Complete to find the redundant column
counts matrix should be integers	Salmon/kallisto/RSEM counts are fractional	Use DESeqDataSetFromTximport() -- it rounds AND carries length offsets
Wrong sign of LFC vs expected	Reference level set alphabetically	relevel() BEFORE DESeq()
padj = NA for biologically meaningful gene	Independent filtering or Cook's outlier	See padj=NA section
LRT LFC doesn't match Wald LFC for the same comparison	LRT reports last coefficient; Wald reports the named coefficient	Extract specific Wald per level for the effect size
summary(res) shows fewer DE genes than expected	summary() default alpha=0.1, NOT the alpha passed to results()	summary(res, alpha = 0.05)

References

• Anders S, Huber W. 2010. Differential expression analysis for sequence count data. Genome Biol 11(10):R106. doi:10.1186/gb-2010-11-10-r106
• Love MI, Huber W, Anders S. 2014. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15(12):550. doi:10.1186/s13059-014-0550-8
• Zhu A, Ibrahim JG, Love MI. 2019. Heavy-tailed prior distributions for sequence count data: removing the noise and preserving large differences. Bioinformatics 35(12):2084-2092. doi:10.1093/bioinformatics/bty895
• Stephens M. 2017. False discovery rates: a new deal. Biostatistics 18(2):275-294. doi:10.1093/biostatistics/kxw041
• Ignatiadis N, Klaus B, Zaugg JB, Huber W. 2016. Data-driven hypothesis weighting increases detection power in genome-scale multiple testing. Nat Methods 13(7):577-580. doi:10.1038/nmeth.3885
• Bourgon R, Gentleman R, Huber W. 2010. Independent filtering increases detection power for high-throughput experiments. PNAS 107(21):9546-9551. doi:10.1073/pnas.0914005107
• McCarthy DJ, Smyth GK. 2009. Testing significance relative to a fold-change threshold is a TREAT. Bioinformatics 25(6):765-771. doi:10.1093/bioinformatics/btp053
• Soneson C, Love MI, Robinson MD. 2015. Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences. F1000Res 4:1521. doi:10.12688/f1000research.7563.2
• Schurch NJ et al. 2016. How many biological replicates are needed in an RNA-seq experiment and which differential expression tool should you use? RNA 22(6):839-851. doi:10.1261/rna.053959.115
• Crowell HL et al. 2020. muscat detects subpopulation-specific state transitions from multi-sample multi-condition single-cell transcriptomics data. Nat Commun 11:6077. doi:10.1038/s41467-020-19894-4
• Jiang L et al. 2011. Synthetic spike-in standards for RNA-seq experiments. Genome Res 21(9):1543-1551. doi:10.1101/gr.121095.111

Related Skills

• edger-basics - Cross-check or use when n<5/group; QL F-test framework
• de-results - padj=NA handling, IHW, TREAT, GSEA input preparation, gene annotation
• de-visualization - MA, volcano (with shrunken LFC), PCA, heatmap, dispersion plot
• batch-correction - Include batch in design vs Nygaard 2016 cardinal sin
• timeseries-de - DESeq2 LRT with splines for time-course
• expression-matrix/counts-ingest - tximport, featureCounts, STAR output decisions
• expression-matrix/normalization - RLE/TMM/VST/rlog mechanics and failure modes
• expression-matrix/metadata-joins - Reference level, paired design, interaction parameterization
• expression-matrix/gene-id-mapping - Annotating DE results with symbols
• rna-quantification/tximport-workflow - Detailed tximport mechanics
• pathway-analysis/gsea - Ranked-list input from DE
• pathway-analysis/go-enrichment - ORA with proper background
• data-visualization/volcano-and-ma-plots - Custom volcano with apeglm-shrunken LFC