DESeq2 vs edgeR vs limma-voom: Which to Use for RNA-seq (by Sample Size, with R Code)

If you've started learning RNA-seq analysis, you've probably hit the same wall every researcher does: there are three "standard" tools for differential expression — DESeq2, edgeR, and limma-voom — and the documentation for each insists it's the best. Which should you actually use? This guide cuts through the noise with a practical, statistics-grounded comparison: how each tool models your data, when one outperforms the others, and a clear decision framework for choosing the right tool for your specific experiment.

TL;DR

• I ran all three on the public airway dataset (4 vs 4). They agreed on most genes but not all: 65-79% of the combined DEG list was called by all three. Getting the design right (pairing by cell line) changed the results far more than the choice of tool.
• DESeq2: Best general-purpose default; conservative, well-documented, easiest for beginners.
• edgeR: Slightly more sensitive for small sample sizes (3-5/group); excellent for complex experimental designs.
• limma-voom: Fastest, most flexible model formulas, best for large studies (>50 samples) and continuous covariates.
• Never use FPKM/TPM as input. DESeq2 and edgeR need counts; limma-voom also accepts non-integer estimated counts (e.g. from tximport).
• Consider RUVSeq or SVA if you suspect hidden technical confounds.

Why Differential Expression Analysis Is Statistically Hard

RNA-seq differential expression seems simple — compare gene expression between two groups, run a statistical test, pick the significant genes. But three properties of count data make it surprisingly difficult:

1. Counts are non-negative integers, not normally distributed: Standard t-tests assume normality, which fails badly for low-count genes.
2. Mean-variance dependency: Genes with higher mean expression have higher variance — and the relationship is non-linear. Standard parametric tests assume constant variance.
3. Small sample sizes: Most RNA-seq experiments have 3-10 replicates per condition, far too few to estimate gene-specific variance accurately.

DESeq2, edgeR, and limma-voom are all responses to these challenges. Each takes a slightly different statistical approach, leading to slightly different results — but the fundamental problem they solve is the same: how to detect real expression changes when individual genes have unreliable variance estimates from few replicates.

Statistical Foundations

DESeq2: Negative Binomial GLM with Wald or LRT Tests

DESeq2 (Love, Huber & Anders, 2014) models read counts using a Negative Binomial distribution:

counts ~ NegBin(μ_ij, α_i)

Where μ_ij is the expected count for gene i in sample j (related to library size and the gene's true expression), and α_i is the gene-specific dispersion parameter capturing biological variability.

The key insight: DESeq2 shares dispersion information across genes. Instead of estimating each gene's variance independently from few replicates, it pools information across genes to estimate a "trend" (dispersion as a function of mean expression), then shrinks individual gene dispersions toward this trend. This dramatically improves variance estimation accuracy.

For testing, DESeq2 fits a Generalized Linear Model and uses either:

• Wald test (default): Tests a single coefficient (e.g., treatment vs control)
• Likelihood Ratio Test (LRT): Compares two model fits, useful for testing multiple coefficients simultaneously

edgeR: Negative Binomial with Empirical Bayes Dispersion

edgeR (Robinson, McCarthy & Smyth, 2010) also uses a Negative Binomial model but with a different dispersion estimation strategy. Three variants:

• Common dispersion: Single dispersion shared across all genes (rarely used, too simple)
• Trended dispersion: Like DESeq2's mean-dispersion trend
• Tagwise dispersion: Gene-specific, with empirical Bayes shrinkage toward the trend (recommended default)

edgeR offers two main testing approaches:

• Exact test (exactTest): For simple two-group comparisons, computationally fast
• Quasi-likelihood F test (glmQLFTest): For complex designs, uses quasi-likelihood to better account for dispersion uncertainty (recommended for most analyses)

limma-voom: Linear Modeling with Mean-Variance Trend

limma-voom (Law, Chen, Shi & Smyth, 2014) takes a fundamentally different approach. Instead of modeling counts directly, it:

1. Transforms counts to log2 counts-per-million (CPM)
2. Computes the mean-variance trend across all genes
3. Estimates gene/sample-specific weights from this trend
4. Fits linear models with weights (like classical limma for microarrays)
5. Uses empirical Bayes moderation of standard errors

The advantage: linear models are extremely flexible — limma can handle continuous covariates, complex interactions, and time-course designs more naturally than negative binomial GLMs. The disadvantage: the count-to-CPM transformation discards some information, especially for low-count genes.

When Each Tool Has an Advantage

Small Sample Sizes (3 per group)

For very small experiments, the most-cited benchmark is Schurch et al. (2016, RNA), which sequenced 48 replicates per condition in yeast and subsampled them. With three replicates, nine of the eleven tools found only 20–40% of the genes the full data called significant. For fewer than 12 replicates they found edgeR and DESeq2 the leading tools, and they recommended at least six replicates. Law et al. (2014), who introduced voom, reported that limma-voom controls the false discovery rate well at small sample sizes too.

Recommendation: at 3 vs 3, any of the three is defensible. More replicates help more than switching tools. Validate top hits independently.

Standard Sample Sizes (6-20 per group)

For typical experiments, all three agree on the core genes and differ at the edges; see the measured comparison below. The choice often comes down to:

• Familiarity / lab tradition
• Specific design needs (interaction effects, continuous covariates)
• Output format preferences

Recommendation: DESeq2 is a solid default. The shrunken log2FoldChange values it provides via lfcShrink() are particularly useful for visualization and downstream ranking.

Large Studies (>50 samples)

For larger cohorts, especially clinical studies with continuous covariates (age, BMI, tumor stage), limma-voom has clear advantages:

• Much faster computation
• More flexible model formulas
• Better handling of mixed continuous/categorical covariates
• Can correct for known confounders elegantly via design matrix

Recommendation: Use limma-voom. For traditional case/control with categorical groups, DESeq2 or edgeR still work but limma is cleaner.

Single-Cell RNA-seq (Pseudobulk Analysis)

For pseudobulk-aggregated single-cell data (recommended approach for cell-type-specific DE in scRNA-seq):

• DESeq2 or edgeR: Standard choices, performance comparable
• limma-voom: Works but the variance-weight assumption may break with high zero-inflation
• muscat: built specifically for multi-sample, multi-condition single-cell DE

Practical Code Comparison

Same input data, three tools — see how the workflows compare.

Setup (Common to All)

library(SummarizedExperiment)

# Your count matrix: genes (rows) × samples (columns), integer counts
counts <- as.matrix(read.csv("counts.csv", row.names = 1))

# Sample metadata
coldata <- read.csv("samples.csv", row.names = 1)
coldata$condition <- factor(coldata$condition, levels = c("control", "treated"))

DESeq2 Workflow

library(DESeq2)

# Construct DESeq2 dataset
dds <- DESeqDataSetFromMatrix(
  countData = counts,
  colData   = coldata,
  design    = ~ condition
)

# Pre-filter very low-count genes (optional but recommended)
keep <- rowSums(counts(dds) >= 10) >= 3
dds <- dds[keep, ]

# Run differential expression analysis
dds <- DESeq(dds)

# Extract results
res <- results(dds, contrast = c("condition", "treated", "control"))

# Shrink log2 fold changes for better visualization/ranking
# (type = "apeglm" needs the apeglm package: BiocManager::install("apeglm"))
res_shrunk <- lfcShrink(
  dds,
  coef = "condition_treated_vs_control",
  type = "apeglm"
)

# Order by adjusted p-value
res_ordered <- res_shrunk[order(res_shrunk$padj), ]
summary(res_ordered)

# Visualization
plotMA(res_shrunk, ylim = c(-5, 5))

edgeR Workflow

library(edgeR)

# Construct DGEList
dge <- DGEList(counts = counts, group = coldata$condition)

# Filter low-count genes (edgeR's recommended approach)
keep <- filterByExpr(dge, group = coldata$condition)
dge <- dge[keep, , keep.lib.sizes = FALSE]

# TMM normalization (edgeR standard)
dge <- calcNormFactors(dge, method = "TMM")

# Design matrix
design <- model.matrix(~ condition, data = coldata)

# Estimate dispersions (tagwise with trend)
dge <- estimateDisp(dge, design)

# Fit GLM with quasi-likelihood
fit <- glmQLFit(dge, design)

# Test for differential expression
qlf <- glmQLFTest(fit, coef = 2)  # coef 2 = condition effect

# Top DEGs
top_degs <- topTags(qlf, n = Inf, sort.by = "PValue")
summary(decideTests(qlf))

# Visualization
plotMD(qlf)

limma-voom Workflow

library(limma)
library(edgeR)  # for DGEList and filterByExpr

# Same filtering approach as edgeR
dge <- DGEList(counts = counts, group = coldata$condition)
keep <- filterByExpr(dge, group = coldata$condition)
dge <- dge[keep, , keep.lib.sizes = FALSE]
dge <- calcNormFactors(dge, method = "TMM")

# Design matrix
design <- model.matrix(~ condition, data = coldata)

# voom transformation: counts → log2(CPM) with weights
v <- voom(dge, design, plot = TRUE)

# Linear model fitting
fit <- lmFit(v, design)

# Empirical Bayes moderation
fit <- eBayes(fit)

# Top DEGs
top_degs <- topTable(fit, coef = 2, number = Inf, sort.by = "P")
summary(decideTests(fit, p.value = 0.05))

# Visualization
plotSA(fit, main = "Mean-variance trend")

Three Outputs Compared: What Actually Happened

[Rewritten October 2026. An earlier version said to "expect ~80-85% of the same top DEGs" without any measurement. These are measured.]

I ran the three workflows above on the Bioconductor airway dataset (Himes et al. 2014): four airway smooth-muscle cell lines, each untreated and treated with dexamethasone, so 4 vs 4 and naturally paired by cell line. I ran the code twice: once as written above (~ condition), and once with the pairing in the design (~ cell + condition).

FDR < 0.05	DESeq2	edgeR (QL)	limma-voom	In all three	All three / any
~ condition	2,782	1,952	1,878	1,811	65%
~ condition, |log2FC| ≥ 1	828	661	635	594	69%
~ cell + condition	4,099	5,130	5,170	4,055	77%
~ cell + condition, |log2FC| ≥ 1	1,000	1,169	1,155	960	79%

What the numbers say:

• The design mattered more than the tool. Adding the cell line to the model gave DESeq2 1,414 extra DEGs and cost it 97. Ignoring the pairing cost DESeq2 about a third of its calls and edgeR and limma-voom more than half.
• Which tool calls "more" depends on the design. Without pairing DESeq2 called 40% more genes than the other two; with pairing it called the fewest. Any rule like "DESeq2 is conservative" did not hold up here.
• The core is shared. Almost every gene DESeq2 called in the paired design was also called by the other two (4,055 of 4,099).
• Speed is not a reason to choose. On this machine each tool took between 0.4 and 3.1 seconds.

If your samples are paired (same donor, cell line or batch in both conditions), put that in the design. Here, getting the design right helped far more than switching tools did.

Critical Pre-Analysis Steps (Apply to All Three)

The biggest factor in DE analysis quality isn't tool choice — it's whether you've done these prerequisites correctly.

1. Use Raw Counts, Never FPKM/TPM

DESeq2 and edgeR need count data as input; limma-voom also works with non-integer estimated counts. Common mistake: feeding FPKM, TPM, or already-normalized values produces nonsense results.

If you only have FPKM/TPM, you should re-quantify from raw FASTQ files using salmon, kallisto, or STAR-htseq-count. There is no statistically sound way to "convert" FPKM to counts.

2. Quality Control Plots Before DE

Always run before testing:

• PCA plot on vst() (DESeq2) or cpm(log = TRUE) (limma) to detect outliers and batch effects
• Sample distance heatmap to identify mislabeled samples
• Library size visualization to spot failed sequencing runs
• Count distribution plots (boxplots, density curves) per sample

If PCA shows obvious batch effects or outliers, address them before running DE — otherwise your DEGs may be technical artifacts.

3. Batch Effect Correction

If samples were processed in batches (different days, kits, sequencing runs), include batch as a covariate:

# DESeq2
design = ~ batch + condition  # condition tested while controlling for batch

# edgeR / limma
design <- model.matrix(~ batch + condition, data = coldata)

For unknown technical confounders, use SVA (svaseq) or RUVSeq to estimate hidden surrogate variables, then include them in the design.

A batch term only helps when every batch contains both conditions; for what happens when run order lines up with the groups, see injection order in LC-MS.

4. Sensible Filtering

Pre-filter genes with very low counts before DE:

• Light filtering: e.g. genes with ≥10 counts in at least as many samples as your smallest group. This is a common pre-filter, not a DESeq2 default; DESeq2's own default is independent filtering inside results()
• Stricter filtering (edgeR filterByExpr): Considers minimum group size and library composition
• Don't over-filter: Removing too many genes inflates the multiple testing penalty for the rest

5. Check Multiple Testing Correction

All three tools use Benjamini-Hochberg (BH) FDR by default. For most analyses this is appropriate. For very strict applications (clinical biomarkers), consider Bonferroni (more stringent) or the qvalue package, which gains power over BH when a sizeable fraction of genes are truly changed (when almost everything is null, it and BH agree).

Common Pitfalls Across All Three Tools

Treating log2FoldChange Magnitudes as Reliable for Low-Expression Genes

A gene with 1 count in control vs 10 in treated has a calculated log2FC = 3.3, but this is meaningless statistical noise. DESeq2's lfcShrink() partially solves this by shrinking small effects toward zero. Always use shrunken estimates for ranking and visualization.

Confusing Statistical Significance with Biological Importance

A gene with FDR = 1e-50 but log2FC = 0.05 is statistically significant but probably biologically irrelevant. Use a fold-change threshold as well as an FDR threshold. Filtering on |log2FC| after testing is common; testing against the threshold directly (results(dds, lfcThreshold = 1) in DESeq2, glmTreat() in edgeR, treat() in limma) is the statistically cleaner version.

Multiple Testing Across Multiple Comparisons

If you run 5 separate condition comparisons, the FDR within each is controlled — but the family-wise error rate across all 5 is not. For studies with many comparisons, consider hierarchical FDR control or a study-wise correction.

Forgetting to Set Reference Level

The reference level of your factor determines which group is the "control" in fold change calculations:

# Without setting levels: alphabetical order — "control" first by chance, but...
coldata$condition <- factor(coldata$condition)  # alphabetical: A, B, C, ...

# Always set explicitly to avoid surprises:
coldata$condition <- factor(coldata$condition, levels = c("control", "treated"))

If your fold changes are in the wrong direction, this is usually why.

Decision Framework Summary

Here's a simple decision tree for choosing your DE tool:

Start with: How many samples per group?

• 3 per group: edgeR (glmQLFTest) or DESeq2 — both work, edgeR may be slightly more sensitive
• 5-20 per group: DESeq2 (sensible default) or edgeR
• 20+ per group: limma-voom (faster, more flexible)
• 50+ per group with continuous covariates (age, BMI, etc.): limma-voom (clearly best)

Then check: Is your design simple or complex?

• Two-group comparison: Any of the three works
• Multi-factor (interaction terms, batch + condition + sex): All three handle this; limma's syntax is cleanest
• Time course or continuous predictor: limma-voom strongly preferred
• Pseudobulk single-cell: DESeq2 or edgeR

Finally: Are you new to RNA-seq?

• Yes: Use DESeq2 — best documentation, largest user community, most tutorials, friendly error messages

Beyond Differential Expression

DEG analysis is the start, not the end. Once you have your DEG list:

1. Quality check: Do top DEGs make biological sense? Use volcano plot, MA plot, heatmap of top 50.
2. Functional analysis: GO + Pathway enrichment to understand biological themes (covered in Part 2 of our functional analysis series)
3. Network analysis: PPI Hub identification to find regulatory drivers (Part 1 of our series)
4. TF activity: DoRothEA + viper to identify upstream regulators
5. Validation: qPCR, Western blot, IHC for top candidates

Conclusion

Choosing between DESeq2, edgeR, and limma-voom is less consequential than people think. All three are mature, well-validated tools that produce comparable results for typical experiments. What matters far more:

• Correct experimental design and sample size
• Proper raw count input (not FPKM/TPM)
• Thoughtful batch effect handling
• Meaningful filtering thresholds
• Appropriate fold-change cutoffs alongside FDR

Pick a tool, learn its idioms well, and focus your energy on the biology and the experimental design. Then move to the downstream functional analysis that turns a DEG list into biological insight.

Code and outputs: the scripts and result tables for this post are on GitHub in sbmlab-benchmarks/deseq2-edger-limma-airway.

Further Reading

• Love, M. I., Huber, W. & Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15(12).
• Robinson, M. D., McCarthy, D. J. & Smyth, G. K. (2010). edgeR: A Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics, 26(1).
• Law, C. W., Chen, Y., Shi, W. & Smyth, G. K. (2014). voom: Precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biology, 15(2).
• Soneson, C. & Delorenzi, M. (2013). A comparison of methods for differential expression analysis of RNA-seq data. BMC Bioinformatics, 14(1).
• Schurch, N. J. 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).