The Tool Choice That Decides Your Gene List
The Tool Choice That Decides Your Gene List
DESeq2, edgeR, and limma-voom will hand you three different answers from the same count matrix. Here is why they disagree, when each is right, and how to build a differential expression workflow you can defend.
Run the same RNA-seq count matrix through DESeq2, edgeR, and limma-voom and you will get three overlapping but non-identical lists of differentially expressed genes. The overlap is usually large — often 80 to 90 percent for a well-powered experiment with a clean design — but the disagreement at the margins is exactly where biology gets reported, follow-up experiments get funded, and clinical claims get made. A gene that clears significance in one tool and misses it in another is not a rounding error; it is a decision the analyst made by choosing a method, whether or not they realized they were making it. For a field that treats the differentially expressed gene list as a near-objective readout of an experiment, this dependence on tool choice is uncomfortably underexamined.
The discomfort is not a reason for paralysis. The three dominant methods are all statistically principled, all well-validated, and all converge on the same answer when the experiment is well-designed and adequately powered. They diverge in predictable ways under specific, identifiable conditions — small sample sizes, heavy-tailed dispersion, outlier-driven counts, complex multi-factor designs — and a practitioner who understands those conditions can choose the right tool deliberately rather than by habit or by whichever vignette they read first. This article unpacks the shared statistical foundation, profiles where the three methods part ways, and lays out the workflow decisions that separate a reproducible differential expression analysis from one that quietly depends on undocumented choices.
The Shared Foundation Everyone Forgets Is Shared
All three methods solve the same underlying problem, and they solve it with more in common than their competing reputations suggest. RNA-seq produces count data — the number of reads assigned to each gene in each sample — and count data has two properties that break ordinary linear regression. First, the variance grows with the mean: a gene averaging 10,000 counts varies far more in absolute terms than one averaging 10. Second, biological replicates show more variability than a simple Poisson count process would predict, a phenomenon called overdispersion. Both DESeq2 and edgeR address this by modeling counts with a negative binomial distribution, which adds a dispersion parameter on top of the Poisson mean to absorb that extra biological variability.1
The central statistical challenge is that a typical experiment has far too few replicates to estimate each gene's dispersion reliably on its own. With three samples per group, the per-gene variance estimate is so noisy that significance testing built on it would be hopeless. The shared insight that makes modern RNA-seq analysis work is information borrowing: the methods estimate dispersion for each gene partly from that gene's own data and partly from the trend across all genes with similar expression levels, shrinking the noisy individual estimates toward a stable, model-fitted expectation. This empirical Bayes shrinkage — pioneered in the microarray era and carried into sequencing — is the single most important reason these tools can find anything at all with three replicates.2
limma-voom takes a different route to the same destination. Rather than model raw counts, it transforms them to log-counts-per-million and estimates a precision weight for each individual observation that captures the mean-variance relationship, then runs the same well-tested linear modeling and empirical Bayes machinery that limma has used for microarrays for two decades. The counts never enter a negative binomial model at all; the weighting scheme makes the transformed data behave well enough for normal-theory linear models.3 Three roads, one mountain.
Where the Three Methods Actually Diverge
If the foundations are shared, the divergence comes from the modeling choices layered on top — and those choices have predictable consequences that an analyst can anticipate rather than discover after the fact.
DESeq2: conservative by construction
DESeq2's defining features are its outlier handling and its independent filtering. It flags genes whose counts in one sample are extreme relative to the fitted model using a Cook's distance cutoff, and by default it replaces or removes those outlier-driven calls rather than reporting them as significant. It also applies log-fold-change shrinkage, pulling the estimated effect sizes of low-information genes toward zero so that a gene with two counts in one group and zero in the other does not post an enormous, meaningless fold change. The net effect is a tool biased toward conservatism: DESeq2 tends to report fewer genes than edgeR at the same nominal threshold, with those it does report being more robust to single-sample artifacts. For small-sample clinical work where a false positive is expensive, this conservatism is a feature.1
edgeR: sensitive, with a quasi-likelihood upgrade
edgeR historically reported more genes than DESeq2 at the same threshold, trading a little specificity for sensitivity. The modern recommended workflow — the quasi-likelihood F-test, glmQLFit followed by glmQLFTest — substantially tightened edgeR's error-rate control by accounting for the uncertainty in dispersion estimation itself, rather than treating the estimated dispersion as known. The quasi-likelihood pipeline is now widely regarded as among the best-calibrated approaches for experiments with very few replicates, and is the version that belongs in any current edgeR analysis. The legacy exact test and likelihood-ratio test remain in the package but should generally be retired from new work.4
limma-voom: the scaling and flexibility champion
limma-voom's advantages show up at the two extremes of experimental complexity. For very large studies — hundreds or thousands of samples — the linear-model machinery is dramatically faster than refitting negative binomial models gene by gene, and the normal-theory framework makes large-cohort designs computationally tractable where the count-based methods strain. At the other end, for complex designs with multiple factors, interaction terms, blocking variables, or repeated measures, limma's mature design-matrix and contrast infrastructure — including duplicateCorrelation for handling correlated samples — is more flexible and better documented than its competitors. The cost is a mild loss of power at the smallest sample sizes, where the count-based models' explicit treatment of low counts gives them an edge.3
| Method | Model | Strength | Reach for it when |
|---|---|---|---|
| DESeq2 | Negative binomial + shrinkage | Outlier robustness; conservative calls | Small clinical studies; messy or single-outlier-prone data |
| edgeR (QL) | Negative binomial + quasi-likelihood F | Best-calibrated error rate at low n | Few replicates where FDR control is paramount |
| limma-voom | Log-CPM + precision weights + linear model | Speed; complex/blocked designs | Large cohorts; multi-factor, repeated-measures, or batch-blocked designs |
The Decisions That Matter More Than the Tool
Here is the uncomfortable truth that the tool-comparison literature often buries: the choice between DESeq2, edgeR, and limma-voom is usually a smaller source of variation in your results than several upstream decisions that analysts make without much deliberation. A defensible workflow controls these first.
Normalization is not optional and not interchangeable. Raw counts cannot be compared across samples because total sequencing depth differs and a handful of highly expressed genes can dominate a library. The composition-aware normalization methods — DESeq2's median-of-ratios and edgeR's TMM — both correct for this and tend to agree closely, but neither resembles a naive counts-per-million scaling, which assumes total RNA output is constant across conditions. That assumption fails exactly when it matters most, in experiments where a treatment globally shifts transcription. Choosing CPM scaling when a composition-aware method is warranted will distort the entire downstream analysis regardless of which DE tool runs afterward.5
Filtering low-count genes changes your false discovery rate. Genes with near-zero counts across all samples carry no information but inflate the multiple-testing burden, weakening the correction applied to the genes you care about. Filtering them out before testing — edgeR's filterByExpr is the standard implementation — measurably increases power. But the filtering threshold is a parameter, and an undocumented, ad hoc threshold is a reproducibility hole. State it explicitly.
The fold-change threshold is doing statistical work. Applying a log-fold-change cutoff after significance testing — keeping only genes that are both statistically significant and biologically substantial — is common and reasonable, but doing it as a post hoc filter on a p-value-ranked list does not control the error rate for the combined criterion. DESeq2 and edgeR both provide formal tests against a fold-change threshold (lfcThreshold in DESeq2; glmTreat in edgeR) that fold the effect-size requirement into the hypothesis test itself, which is the statistically correct way to ask for genes that are both real and large.1
Batch effects belong in the model, not in your imagination. If samples were processed in distinct batches, on different days, or by different operators, that structure must enter the design matrix as a covariate. Removing batch effects by eye, or worse, ignoring them, is one of the most common ways a published RNA-seq result fails to replicate. All three frameworks accommodate batch covariates natively; the failure is analytical discipline, not tooling.
Building a Workflow You Can Defend
A production-grade differential expression analysis — the kind that survives peer review, supports a regulatory submission, or anchors a clinical biomarker claim — is defined less by its choice of DE package than by the reproducibility scaffolding around it. Several practices separate a defensible workflow from a fragile one.
Pin every version. DESeq2, edgeR, and limma evolve, and so do their defaults; the median-of-ratios estimator, the quasi-likelihood machinery, and the voom weighting have all changed across releases. A result generated under one Bioconductor release may not reproduce numerically under another. Capturing the full software environment — every package version and the R version — in a container or lockfile is the only way to guarantee that the analysis you ran in March reproduces in December.
Run the sanity checks before the test. A principal component analysis of the normalized, variance-stabilized counts should show samples clustering by biological condition rather than by batch or library size before any differential testing happens. A sample that lands in the wrong cluster is a swap, a contamination, or a batch confound — and discovering it after you have a gene list is discovering it too late. Dispersion plots and mean-variance trend plots should look like the textbook examples; deviations signal a model that does not fit the data.
Pre-specify the comparison. The design matrix, the contrasts of interest, the significance threshold, and the fold-change criterion should be written down before the data are unblinded. RNA-seq offers an enormous garden of forking paths — which samples to include, which threshold to apply, which of three tools to trust — and the way to keep that flexibility from manufacturing false positives is to commit to the analysis plan in advance. When a result genuinely requires deviating from the plan, document the deviation and its rationale.
When two tools disagree on a gene that matters, look at the counts. The most useful response to a borderline gene that clears significance in edgeR but not DESeq2 is not to pick the answer you prefer — it is to plot the per-sample normalized counts and look. Disagreements at the margin are almost always explainable: a single high-leverage sample, a low-count gene at the edge of the filter, a fold change inflated by a near-zero denominator. The tools are not being capricious; they are encoding different, defensible judgments about how much to trust thin data. Reading the counts directly turns a tool disagreement into a biological decision you can actually justify.
The Honest Summary
The differentially expressed gene list is one of the most consequential outputs in modern molecular biology, and it depends on choices that are too often made implicitly. DESeq2, edgeR, and limma-voom are all excellent, all principled, and all correct in the regimes they were built for — the failure mode is not choosing the wrong one but failing to recognize that a choice was made at all. For a well-designed, adequately powered experiment, any of the three will land you in nearly the same place, and the upstream decisions about normalization, filtering, and batch structure will move your results more than the package name on the call. For the small-sample, outlier-prone, or complexly designed experiments where the methods genuinely diverge, the right move is to understand why they diverge, choose deliberately, and document the choice so the next analyst — or the reviewer, or the regulator — can follow the reasoning.
That is, in the end, what a well-scoped transcriptomics engagement is for: not running the tool, which anyone can do, but making the dozen upstream choices defensibly, documenting them so the result reproduces, and knowing which of the three answers to trust when they disagree.
Zetobit LLC Bioinformatics Insight Series · June 2026
Lexington, KY · Genomics · Transcriptomics · Proteomics
Publications
- Love M.I., Huber W. & Anders S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15, 550. doi:10.1186/s13059-014-0550-8
- Smyth G.K. (2004). Linear models and empirical Bayes methods for assessing differential expression in microarray experiments. Statistical Applications in Genetics and Molecular Biology, 3, Article 3. doi:10.2202/1544-6115.1027
- 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, R29. doi:10.1186/gb-2014-15-2-r29
- Chen Y., Lun A.T.L. & Smyth G.K. (2016). From reads to genes to pathways: differential expression analysis of RNA-seq experiments using Rsubread and the edgeR quasi-likelihood pipeline. F1000Research, 5, 1438. doi:10.12688/f1000research.8987.2
- Robinson M.D. & Oshlack A. (2010). A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biology, 11, R25. doi:10.1186/gb-2010-11-3-r25
- 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), 139–140. doi:10.1093/bioinformatics/btp616
- Soneson C. & Delorenzi M. (2013). A comparison of methods for differential expression analysis of RNA-seq data. BMC Bioinformatics, 14, 91. doi:10.1186/1471-2105-14-91