Skip to content
RNA-Seq

How to Use DESeq2 for Differential Expression Analysis

A rigorous walkthrough of the DESeq2 workflow: count matrix prep, dispersion estimation, LFC shrinkage, batch correction, and pitfalls that get past most tutorials.

PR Dr. Lin Wei 12 min read DESeq2 Differential Expression R

DESeq2 is the most cited differential expression method in RNA-Seq for a reason: its negative-binomial GLM with empirical Bayes shrinkage is both statistically principled and empirically robust across thousands of published studies. But most tutorials stop at results(dds). This guide covers the parts that actually matter in a real analysis: dispersion diagnostics, LFC shrinkage, batch correction, and interpreting adjusted p-values honestly.

The mental model

DESeq2 asks a single question, one gene at a time:

Given the count distribution I observe across all samples, is the change in this gene’s expression between condition A and condition B larger than I’d expect by chance under a negative-binomial model with shared dispersion trend?

Everything else — normalization, size factors, LFC shrinkage — supports that question.

1. From counts to a DESeqDataSet

Assume you already have a raw integer count matrix (from featureCounts, HTSeq-count, or tximport(txi, "gene")):

library(DESeq2)

counts <- as.matrix(read.table("counts/gene_counts.tsv",
                               header = TRUE, row.names = 1))
coldata <- read.csv("meta/samples.csv", row.names = 1)
stopifnot(all(colnames(counts) == rownames(coldata)))

dds <- DESeqDataSetFromMatrix(
  countData = counts,
  colData   = coldata,
  design    = ~ batch + condition
)

Note two design-formula rules that trip up beginners:

  • The variable of interest goes last. DESeq2 reports the LFC for the last variable by default.
  • Every level in every factor must have at least one replicate. If not, that coefficient is unestimable.

2. Pre-filtering low-count genes

Removing genes with almost no reads speeds up estimation and improves multiple-testing power:

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

This keeps genes with ≥10 counts in at least three samples — a good default for a 6-sample study. Do not filter more aggressively than this before DESeq(); DESeq2’s independent filtering step (built into results()) handles the rest.

3. Run the model

dds <- DESeq(dds)

Under the hood, DESeq() does three things:

  1. Estimate size factors (median-of-ratios normalization).
  2. Estimate dispersions per gene, then shrink toward a fitted mean-dispersion trend.
  3. Fit a negative-binomial GLM and run a Wald test on each coefficient.

Always look at the dispersion plot. If the fit is wrong, everything downstream is wrong:

plotDispEsts(dds)

Healthy plot: a scatter of black points (raw estimates) around a red curve (fit), with blue points (shrunken estimates) hugging the curve. Diagnostic red flags: a cloud that doesn’t decline with mean expression (batch effects), or the red curve pinned at a very high dispersion floor (extreme sample-to-sample noise).

4. Extract results, honestly

res <- results(dds,
               contrast = c("condition", "trt", "ctrl"),
               alpha    = 0.05)
summary(res)

summary(res) will report:

  • Total genes tested after independent filtering
  • Number with adjusted p < 0.05 up and down
  • Genes filtered out as low mean count
  • Outliers detected via Cook’s distance

5. LFC shrinkage — do it, always

Raw log2 fold-changes are noisy for low-count genes. lfcShrink() uses an empirical Bayes prior to shrink LFCs toward zero when the data don’t strongly support a large effect:

resLFC <- lfcShrink(dds, coef = "condition_trt_vs_ctrl", type = "apeglm")

Compare plotMA(res) and plotMA(resLFC) — the second is what you report and what you rank by for GSEA. apeglm is the modern default; use ashr for arbitrary contrast= arguments.

6. Handling batch effects

If your samples were sequenced in batches, hard-code the batch in the design:

design(dds) <- ~ batch + condition
dds         <- DESeq(dds)

DESeq2 will estimate a per-batch mean shift for every gene, and the reported LFC for condition is orthogonal to batch. For visualization (PCA, heatmaps) apply limma::removeBatchEffect() to the vst() transformed data.

vsd <- vst(dds, blind = FALSE)
assay(vsd) <- limma::removeBatchEffect(assay(vsd), vsd$batch)
plotPCA(vsd, intgroup = "condition")

7. Interpreting adjusted p-values

DESeq2 uses the Benjamini-Hochberg procedure over the genes that pass its independent filter. A gene with padj < 0.05 means: among the set of significant genes at this threshold, we expect ≤5 % to be false discoveries.

It does not mean the gene is real biology. Always inspect:

  • Sample-level counts (plotCounts(dds, gene = "ENSG00000141510", intgroup = "condition"))
  • Cook’s distance flags (res$maxCooks)
  • Whether the effect is driven by one outlier sample

8. Exporting and reporting

Save the results table with useful annotations:

library(AnnotationDbi); library(org.Hs.eg.db)
res_df <- as.data.frame(resLFC) |>
  tibble::rownames_to_column("ensembl_gene_id") |>
  dplyr::mutate(symbol = mapIds(org.Hs.eg.db,
                                keys = ensembl_gene_id,
                                column = "SYMBOL",
                                keytype = "ENSEMBL",
                                multiVals = "first"))
write.csv(res_df, "results/deseq2_trt_vs_ctrl.csv", row.names = FALSE)

Report the following in your methods:

  • DESeq2 version and R version
  • Full design formula
  • Filtering thresholds
  • Shrinkage method (apeglm, ashr, or normal)
  • Independent filtering: on (default) or off

Common mistakes

  • Passing normalized counts to DESeqDataSetFromMatrix(). Only raw integers are valid input.
  • Filtering after DESeq() to boost significance — this invalidates the FDR control.
  • Running DE separately per batch and hoping to combine results. Instead, fit one model with a batch term.
  • Reporting log2FoldChange from results() instead of from lfcShrink(). Ranked lists based on raw LFC are dominated by low-count noise.

DESeq2 rewards patience. Read the dispersion plot, use shrinkage, keep the whole design in one model, and your DE list will replicate.

FAQ

Q. Should I use lfcShrink() before or after filtering?

A. Always run DESeq() on the full dataset first, then filter results by adjusted p-value and shrunken LFC. Filtering the count matrix before dispersion estimation biases the mean-variance relationship that DESeq2 relies on.

Q. What's the difference between apeglm, ashr, and normal shrinkage?

A. apeglm is the modern default and works for any contrast that can be expressed as a coefficient (use coef=). ashr uses a mixture-of-normals prior and works with contrast=. normal is the original DESeq2 shrinkage — still fine, but overshrinks large effects. Prefer apeglm for two-level factors and ashr for arbitrary contrasts.

Q. Do I need to normalize my counts before feeding them to DESeq2?

A. No. DESeq2 performs its own median-of-ratios normalization internally. Pass raw integer counts. Passing TPMs or FPKMs invalidates every downstream test.

Related in RNA-Seq