Differential Gene Expression Analysis Done Right

Q: Can I run DESeq2 or edgeR on TPM or FPKM values?

No. Both require un-normalised counts or estimated counts; DESeq2’s documentation states that “transformed or normalized values such as counts scaled by library size should not be used as input.” The model needs the raw count scale to judge measurement precision, and it handles library size internally through offsets. If you only have TPMs, go back to the quantifier output and import the estimated counts with tximport.

Q: Why do some of my genes have an NA adjusted p-value?

Two causes. Independent filtering removes low-expression genes from the multiple-testing adjustment to increase power for the rest, setting their padj to NA. Separately, genes with an extreme count in one sample are flagged by Cook’s distance and returned as NA, unless you have seven or more replicates per group, in which case the outlying value is imputed instead. NAs of the second kind clustering in one sample usually mean that sample is an outlier.

Q: Should I filter my results by fold change as well as adjusted p-value?

Not after the fact. Use the built-in threshold test instead — lfcThreshold in DESeq2’s results() or treat() in limma — which tests the composite null that the effect is no larger than your threshold and therefore controls error rates against the hypothesis you actually mean. Post-hoc magnitude filtering of a p-value list does not.

Q: How many biological replicates do I need?

The 48-replicate benchmark recommends at least six per condition, rising to at least 12 when you need to catch differentially expressed genes across all fold-change sizes; three replicates recovered only 20–40% of the genes found with the full set for most tools. Three remains publishable for large effects, but describe the resulting list as a high-confidence subset rather than a complete one.

Q: Does it matter whether I use DESeq2, edgeR or limma-voom?

Much less than your design and replicate count. In the 48-replicate benchmark, edgeR and DESeq2 gave the best true-positive/false-positive trade-off below 12 replicates; limma-voom is documented as giving accurate type I error control at small sample sizes and is comfortable with complex designs. Pick one, and do not run all three and report whichever gives the longest list.

Q: Can I compare the DE gene lists from two different comparisons?

Not as a test. A gene reaching significance in one comparison and not the other is not evidence that the effects differ — that requires an interaction term in a single model, tested directly, as in a design of the form ~ group + treatment + group:treatment . Overlapping gene lists are a description, not an inference.

Q: Do I need a different method for single-cell data?

For comparisons between biological conditions, aggregate cells within each replicate into pseudobulk profiles and use the bulk tools. The six best-performing methods in a benchmark against eighteen gold-standard datasets all did this; cell-level tests treat cells from one donor as independent replicates and generated hundreds of differentially expressed genes in comparisons where no perturbation existed.

A differential gene expression (DGE) test is not a comparison of two expression values. It is a statement about a population of biological units — animals, donors, independently grown cultures — inferred from a handful of them. Almost every way a DGE analysis goes wrong is a way of losing sight of that: feeding the model numbers that no longer carry counting precision, writing a design formula that answers a different question than the one asked, ranking genes by a fold change the data cannot support, or treating cells as if they were replicates. The software will return a tidy table of significant genes in all four cases.

This guide works through the analysis in the order the decisions actually bind: what goes into the model, what the model is, how the design formula encodes your question, how effect sizes are estimated, and how to read every column of the output — including the blank ones.

Give the model raw counts, and only raw counts

The single most common fatal input error is handing DESeq2 or edgeR a matrix of TPM, FPKM, RPKM, or counts already divided by library size. The DESeq2 documentation is explicit: “The values in the matrix should be un-normalized counts or estimated counts of sequencing reads…transformed or normalized values such as counts scaled by library size should not be used as input.” The Bioconductor gene-level workflow puts the reason plainly — “the statistical model is most powerful when applied to un-normalized counts, and is designed to account for library size differences internally…only counts allow assessing the measurement precision correctly.”

That last clause is the whole point. A count of 8 and a count of 8,000 carry very different Poisson-level uncertainty. Once you divide by library size, that information is destroyed, and the model can no longer tell a confidently measured small difference from a noisily measured large one.

Normalisation is internal, and it is not just depth

Both tool families estimate a per-sample scaling factor and use it as an offset inside the model rather than rescaling the counts. DESeq2 uses the median-of-ratios method; edgeR uses the trimmed mean of M-values (TMM), applied through calcNormFactors. Neither is simply “divide by total reads,” and the reason is composition bias. Robinson and Oshlack’s original TMM paper states it directly: “if a large number of genes are unique to, or highly expressed in, one experimental condition, the sequencing ‘real estate’ available for the remaining genes in that sample is decreased.” Total-count scaling then makes every other gene in that sample look coordinately down-regulated even when nothing changed.

If you quantified with Salmon, kallisto or RSEM

Transcript-level quantifiers emit estimated counts and effective transcript lengths. Import them with tximport/tximeta rather than summing a TPM column, because the import carries the length information into the model as an offset and, in the Bioconductor workflow’s words, “corrects for any potential changes in gene length across samples (e.g. from differential isoform usage).” A gene whose dominant isoform switches between conditions changes its effective length; without the offset, that shows up as a spurious expression change.

The model: a negative binomial GLM that borrows strength across genes

RNA-seq counts are over-dispersed — their variance exceeds their mean, because biological variation between samples adds to counting noise. DESeq2 describes counts for gene i in sample j with “a GLM of the negative binomial family with a logarithmic link,” giving a variance of µ + αµ², where α is the gene’s dispersion.

Estimating α well is the hard part, and it is why you cannot replace this with a t-test on normalised values. As the DESeq2 authors put it: “A core challenge is the small number of samples in typical HTS experiments – often as few as two or three replicates per condition. Inferential methods that treat each gene separately suffer here from lack of power.” The fix in both major packages is empirical Bayes shrinkage: DESeq2 assumes “genes of similar average expression strength have similar dispersion,” fits a smooth trend across all genes, and pulls each gene-wise estimate toward that trend, with the strength of pull depending on how tightly genes cluster around the fit and on the available degrees of freedom. edgeR’s quasi-likelihood pipeline does the analogous thing: “The raw QL dispersion estimates are squeezed towards a global trend.”

Which test, and does the package choice matter

Wald test (DESeq2 default) — “compares the beta estimate…divided by its estimated standard error…to a standard normal distribution.” Right for a single named contrast.
Likelihood ratio test — compares a full model against a reduced one. Use it when you want a single p-value for “does this factor matter at all” across more than two levels, or across a time course.
edgeR quasi-likelihood F-test — the edgeR authors recommend QL F-tests “instead of the more usual likelihood ratio tests (LRT) as they give stricter error rate control by accounting for the uncertainty in dispersion estimation.”
limma-voom — converts counts to log-CPM, fits a LOWESS mean–variance trend, and turns it into per-observation precision weights so ordinary linear models apply. Its main documented advantage is “accurate type I error rate control even when the number of RNA-seq samples is small.”

For tool choice, the largest empirical benchmark available is Schurch and colleagues’ 48-replicate-per-condition yeast experiment in RNA (2016). Their conclusion: “If fewer than 12 replicates are used, a superior combination of true positive and false positive performances makes edgeR and DESeq2 the leading tools.” Below that replicate count, the design and the input matter far more than which of the three you pick.

The design formula is the experiment, written down

More analyses answer the wrong question here than fail anywhere else. Three rules cover most of it.

1. Know which contrast you are being handed by default

DESeq2’s documentation is unambiguous and routinely missed: “With no additional arguments to results(), the log2 fold change and Wald test p value will be for the last variable in the design formula, and if this is a factor, the comparison will be the last level of this variable over the reference level.” Two consequences follow. Order your formula so the variable of interest comes last, and set the reference level deliberately — “it is prefered in R that the first level of a factor be the reference level (e.g. control, or untreated samples).” R orders factor levels alphabetically by default, which is how analyses end up silently reporting treated-over-treated or, worse, a fold change with the sign inverted throughout the manuscript.

2. Put batch in the design; do not pre-correct the data

If samples were processed in identifiable batches, add the batch factor to the formula (~ batch + condition) so the model estimates and accounts for it while testing condition. Do not run a batch-correction function first and then test the corrected values. The limma documentation says so in its own note on removeBatchEffect: “This function is not intended to be used prior to linear modelling. For linear modelling, it is better to include the batch factors in the linear model.” Pre-correcting hides the degrees of freedom the correction consumed, so the downstream test believes it has more evidence than it does. Batch-corrected values are for plots and clustering; the model gets the covariate.

None of this rescues a batch that is fully confounded with condition — all controls extracted Monday, all treated samples Tuesday. That is a design failure, not an analysis problem, and it is settled before sequencing. See our guide to RNA-seq experimental design through analysis for the blocking and randomisation decisions that prevent it.

3. To ask whether an effect differs between groups, test the interaction

The frequent error is running two separate analyses — treated vs. control in males, treated vs. control in females — then declaring the genes significant in one list and not the other to be sex-specific responders. That comparison is not a test of anything; a gene can sit at p = 0.04 in one group and p = 0.06 in the other with no meaningful difference between the two effects. The Bioconductor workflow gives the correct construction: “if the research aim is to determine for which genes the effect of treatment is different across groups, then interaction terms can be included and tested using a design such as ~ group + treatment + group:treatment.” Test the interaction term and let it produce its own p-value.

Shrink the fold changes before you rank, plot or filter on them

The maximum-likelihood log2 fold change (LFC) for a gene with a handful of reads is close to noise. DESeq2’s authors describe “the strong variance of LFC estimates for genes with low read count,” where weakly expressed genes show exaggerated apparent differences purely from count heteroskedasticity. Shrinkage pulls estimates “toward zero in a manner such that shrinkage is stronger when the available information for a gene is low, which may be because counts are low, dispersion is high or there are few degrees of freedom.”

The current recommended estimator is apeglm, which uses a heavy-tailed Cauchy prior. Its authors’ framing of the problem is worth quoting because it names the two workarounds people reach for instead: “One approach is to introduce filtering thresholds and pseudocounts to exclude or moderate estimated LFCs. Filtering may result in a loss of genes from the analysis with true differences in expression, while pseudocounts provide a limited solution that must be adapted per dataset.” The Cauchy tail is what lets apeglm shrink uninformative genes hard while leaving genuinely large effects essentially untouched, giving “lower bias than previously proposed shrinkage estimators.”

Practically: run the test, then call lfcShrink with type="apeglm", and use the shrunken LFCs for your MA plot, your volcano plot, your ranked list for GSEA, and any figure that shows an effect size. Use the unshrunken p-values for significance. A volcano plot built from unshrunken LFCs has a characteristic flared base of low-count genes with enormous fold changes and no significance — that flare is an artefact, not biology.

Fold-change thresholds belong inside the test, not after it

Filtering a results table to padj < 0.05 & abs(log2FoldChange) > 1 is the near-universal habit, and it does not do what people believe. The p-value was computed against a null of “LFC = 0.” Screening the survivors by magnitude afterwards does not give you error-rate control against the null you actually care about, which is “LFC is no bigger than 1.”

DESeq2 provides the correct construction and explains why it exists: “With sufficient sample size, even genes with a very small but non-zero LFC will eventually be detected as differentially expressed. A change should therefore be of sufficient magnitude to be considered biologically significant.” The package therefore tests composite null hypotheses of the form |β| ≤ θ, letting you “evaluate statistically directly whether there is sufficient evidence that the LFC is above the chosen threshold.” Pass the threshold to results(lfcThreshold = ...); limma users have the equivalent in treat(). The resulting gene list is shorter and defensible.

Reading the results table, including the blanks

baseMean — mean of normalised counts across all samples. Use it as a sanity check: a “significant” gene with a baseMean of 3 deserves a look at the raw counts before it goes in a figure.
log2FoldChange — the effect for the contrast you actually requested. Confirm the direction; do not assume.
lfcSE — standard error of that estimate. Two genes with the same LFC and very different lfcSE are not equivalent findings.
pvalue — unadjusted, per gene. Never the column you filter on.
padj — the adjusted p-value, reported as a Benjamini–Hochberg false discovery rate. Selecting padj < 0.05 means roughly 5% of the genes in that selected list are expected to be false positives. It does not mean each individual gene has a 5% chance of being wrong. Our guides to q-values and the positive false discovery rate and to Bonferroni correction versus FDR cover why FDR rather than family-wise control is the right frame for a 20,000-gene screen.

Why padj is NA — two different causes

Blank adjusted p-values are not a bug, and the two causes mean different things.

Independent filtering. “The results() function automatically performs independent filtering based on the mean of normalized counts for each gene, optimizing the number of genes which will have an adjusted p value below a given FDR cutoff.” Genes below the chosen expression threshold are dropped from the multiple-testing adjustment entirely, which lifts power for everything else, and their padj is set to NA. That is a gene the experiment had no realistic chance of calling, not a gene that was tested and failed.

Cook’s distance outliers. DESeq2 “flags, for each gene, those samples that have a Cook’s distance greater than the 0.99 quantile of the F(p, m−p) distribution.” The behaviour then depends on your replicate count: with two or fewer replicates per condition outliers cannot be detected at all; with seven or more, flagged counts are replaced with imputed values; between those, the flagged gene is removed and returns NA. A cluster of NAs here is a symptom — usually one sample that does not belong, which a PCA on variance-stabilised counts would have shown before you ran the test.

Run the sample-level QC before, not after

The Bioconductor workflow deliberately puts “transformations of the counts in order to visually explore sample relationships” ahead of testing. Use a variance-stabilising transform (VST or rlog) — not raw counts, not log-CPM without moderation — and plot sample distances and PCA. If the first principal component separates extraction day rather than condition, you have found the analysis you actually need to run before any gene list is worth reading. These transformed values are for visualisation only; the test still gets the counts.

“Not significant” is not “not changed”

A gene absent from your DE list may be unchanged, or may simply be one your replicate count could not resolve. This distinction matters most when a reviewer asks you to confirm that a pathway was not affected. Absence from a DE list is not evidence of absence of an effect; that claim needs an equivalence test or an explicit power argument.

How many replicates the evidence actually supports

Schurch and colleagues sequenced 48 biological replicates in each of two conditions in yeast and then subsampled, which makes their numbers unusually concrete. With three biological replicates, “nine of the 11 tools evaluated found only 20%–40% of the significantly differentially expressed (SDE) genes identified with the full set of 42 clean replicates.” That rises to over 85% if you restrict attention to genes changing more than fourfold — but recovering more than 85% of SDE genes regardless of fold change required more than 20 replicates. Their recommendation: “at least six biological replicates should be used, rising to at least 12 when it is important to identify SDE genes for all fold changes.”

Two caveats before you quote that at a budget meeting. It is a yeast study, where biological variability is far lower than in human cohorts or outbred animals, so it is a floor rather than a ceiling. And the practical implication for a three-replicate experiment is not that the results are invalid — the same nine tools controlled FDR at approximately 5% at every replicate count — but that the gene list is badly incomplete, and strongly biased toward large fold changes. Report it as such.

The corollary for pre-filtering: keep genes with a floor of counts in at least as many samples as your smallest group. DESeq2 suggests keeping “only rows that have a count of at least 10 for a minimal number of samples,” noting that “the count of 10 is a reasonable choice for bulk RNA-seq” and that the sample count should be the smallest group size — so a gene expressed only in one condition survives.

Single-cell data: cells are not replicates

If your comparison is between biological conditions, not between cell types within a sample, the bulk logic still governs — and this is where the field’s largest correction has landed. Squair and colleagues surveyed 500 recent scRNA-seq studies and found that 14 DE methods accounted for almost 90% of them, with “the most widely used analysis packages in the field currently employ[ing] DE methods prone to false discoveries by default.” Splitting unperturbed control samples into random groups, cell-level methods “identified hundreds of DE genes in the absence of any perturbation.”

Against eighteen curated gold-standard datasets with matched bulk and single-cell measurements, the six top-performing methods shared one property: they aggregated cells within each biological replicate into pseudobulk profiles before testing. Treating each cell as an independent observation understates gene-expression variance, because cells from the same donor are not independent, and biases the result toward highly expressed genes. This is textbook pseudoreplication — mistaking the measurement unit for the experimental unit. Sum counts per gene per sample per cell type, then run DESeq2 or edgeR on the resulting matrix with donors as replicates.

Before you publish the list

State the software and versions, the design formula verbatim, the contrast, the pre-filtering rule, and whether LFCs are shrunken. Any of these missing makes the analysis unreproducible.
Report the significance rule you used, including any lfcThreshold, not just “padj < 0.05.”
Deposit raw counts and metadata alongside the FASTQ files — the Gene Expression Omnibus is the usual destination and most journals now require the accession before acceptance.
Validate a subset independently where the claim carries weight. RT-qPCR on the same RNA confirms the quantification, not the biology; a fresh cohort tests the finding.

Frequently asked questions

Can I run DESeq2 or edgeR on TPM or FPKM values?

No. Both require un-normalised counts or estimated counts; DESeq2’s documentation states that “transformed or normalized values such as counts scaled by library size should not be used as input.” The model needs the raw count scale to judge measurement precision, and it handles library size internally through offsets. If you only have TPMs, go back to the quantifier output and import the estimated counts with tximport.

Why do some of my genes have an NA adjusted p-value?

Two causes. Independent filtering removes low-expression genes from the multiple-testing adjustment to increase power for the rest, setting their padj to NA. Separately, genes with an extreme count in one sample are flagged by Cook’s distance and returned as NA, unless you have seven or more replicates per group, in which case the outlying value is imputed instead. NAs of the second kind clustering in one sample usually mean that sample is an outlier.

Should I filter my results by fold change as well as adjusted p-value?

Not after the fact. Use the built-in threshold test instead — lfcThreshold in DESeq2’s results() or treat() in limma — which tests the composite null that the effect is no larger than your threshold and therefore controls error rates against the hypothesis you actually mean. Post-hoc magnitude filtering of a p-value list does not.

How many biological replicates do I need?

The 48-replicate benchmark recommends at least six per condition, rising to at least 12 when you need to catch differentially expressed genes across all fold-change sizes; three replicates recovered only 20–40% of the genes found with the full set for most tools. Three remains publishable for large effects, but describe the resulting list as a high-confidence subset rather than a complete one.

Does it matter whether I use DESeq2, edgeR or limma-voom?

Much less than your design and replicate count. In the 48-replicate benchmark, edgeR and DESeq2 gave the best true-positive/false-positive trade-off below 12 replicates; limma-voom is documented as giving accurate type I error control at small sample sizes and is comfortable with complex designs. Pick one, and do not run all three and report whichever gives the longest list.

Can I compare the DE gene lists from two different comparisons?

Not as a test. A gene reaching significance in one comparison and not the other is not evidence that the effects differ — that requires an interaction term in a single model, tested directly, as in a design of the form ~ group + treatment + group:treatment. Overlapping gene lists are a description, not an inference.

Do I need a different method for single-cell data?

For comparisons between biological conditions, aggregate cells within each replicate into pseudobulk profiles and use the bulk tools. The six best-performing methods in a benchmark against eighteen gold-standard datasets all did this; cell-level tests treat cells from one donor as independent replicates and generated hundreds of differentially expressed genes in comparisons where no perturbation existed.