Bulk RNA-seq differential expression done right: DESeq2, edgeR, limma and batch effects
Bulk RNA-seq best practices: QC, Salmon vs STAR, design matrices, DESeq2 vs edgeR vs limma-voom, batch effects, shrinkage and enrichment pitfalls.
Key takeaways
- Check the design for confounding before analysis. If batch equals condition, no method can separate them.
- Model batch in the design matrix; do not feed batch-corrected counts into DESeq2 or edgeR.
- Use shrunken fold changes to rank genes, and the set of expressed genes as the enrichment background.
- Catch sample swaps with sex genes and known markers before you trust any result.
Bulk RNA-seq is mature. The methods are well validated and the software is excellent. Even so, bulk RNA-seq analyses still fail regularly, almost always because of design, QC or interpretation rather than the algorithm. This guide is the checklist our analysts follow.
1. Look at the design before the data
Tabulate every sample by condition, batch (extraction day, library prep, sequencing run), sex, and any pairing (patient, litter, donor). Then ask:
- Is any technical factor perfectly aligned with condition? If all treated samples were prepared on one day and all controls on another, the treatment effect cannot be separated from the day. See our batch-confounded case study.
- Are there enough replicates? Three per group is a minimum; four to six detects moderate effects reliably.
- Is there pairing to exploit? Paired designs (before and after in the same subject) gain a lot of power when modeled.
2. Processing and QC
We follow the nf-core/rnaseq standard: FastQC/fastp, alignment with STAR, quantification with Salmon, and a MultiQC report. Then sample-level checks:
- Library size and complexity, duplication and rRNA content.
- Strandedness confirmed from the data, not assumed.
- Gene body coverage: 3' bias suggests degraded RNA.
- Sex check with XIST and Y-chromosome genes. Mismatches with metadata reveal sample swaps.
- PCA and correlation heatmaps colored by every metadata variable, to see what drives variation.
3. Choosing the statistical model
DESeq2 and edgeR model counts with a negative binomial distribution and share information across genes to estimate dispersion. limma-voom transforms counts and models mean-variance trends with precision weights. All three handle multi-factor designs:
~ batch + condition # batch as covariate
~ patient + treatment # paired design
~ genotype + treatment + genotype:treatment # interaction
Use raw counts (or tximport-summarized estimates), never TPM or FPKM, as input.
What our experts avoid: running ComBat or another batch correction on counts and then passing the corrected values to DESeq2. Corrected values break the count model and understate uncertainty. Put batch in the design formula. Use corrected values only for visualization.
4. Unwanted variation you did not record
When PCA shows structure not explained by any metadata, estimate it with RUVSeq (using control genes or replicate samples) or SVA, and include the factors as covariates. Check that the estimated factors are not your condition in disguise.
5. Results that mean something
- Shrink fold changes with apeglm or ashr. Low-count genes otherwise produce huge, noisy fold changes that dominate rankings.
- Control the false discovery rate (Benjamini-Hochberg), and consider testing against a minimum fold-change threshold (
lfcThreshold) rather than zero when large gene lists are uninformative. - Report effect sizes with uncertainty, not just p-values.
6. Enrichment analysis
- Prefer rank-based GSEA (fgsea) on shrunken fold changes or test statistics over a hard-threshold gene list.
- For over-representation, use the set of expressed, tested genes as the background, not the whole genome. A whole-genome background inflates enrichment for any tissue-expressed pathway.
- Report pathway activity tools (PROGENy, decoupleR) as complementary evidence.
7. Beyond gene-level DE
- Differential transcript usage (DRIMSeq, DEXSeq on Salmon estimates) and splicing (rMATS, LeafCutter).
- Time courses with likelihood-ratio tests and clustering of temporal profiles.
- Deconvolution with single-cell references to separate composition changes from expression changes. This is often critical in tissue samples.
- Co-expression modules (WGCNA) linked to traits.
Common errors checklist
- Batch confounded with condition, discovered after sequencing.
- Batch-corrected values used for testing.
- TPM used as model input.
- Sample swaps not detected.
- Whole-genome background for enrichment.
- Ranking by p-value with unshrunken fold changes.
- Ignoring composition: "upregulated in disease" because more immune cells infiltrated.
How OmicsQuest can help
Our bulk RNA-seq analysis service covers QC through enrichment, isoforms and deconvolution, typically in 8-20 hours for standard projects. See typical costs.
Frequently asked questions
Should I use DESeq2, edgeR or limma-voom?
All three are well validated. DESeq2 and edgeR (quasi-likelihood) are strong defaults for typical experiments with a few replicates; limma-voom is fast, flexible for complex designs and performs well with larger sample sizes. The design and QC matter more than the choice among them.
Salmon or STAR for RNA-seq?
Salmon (pseudoalignment with bias correction) is fast and accurate for gene and transcript quantification, and pairs with tximport. STAR alignment is needed for splicing analysis, variant calling, fusion detection and genome-browser visualization. The nf-core/rnaseq pipeline runs STAR alignment with Salmon quantification, giving both.
Work with OmicsQuest. Related services:
References
- Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.
- Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139-140.
- Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43:e47.
- Zhu A, Ibrahim JG, Love MI. Heavy-tailed prior distributions for sequence count data: removing the noise and preserving large differences. Bioinformatics. 2019;35:2084-2092.
- Soneson C, Love MI, Robinson MD. Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences. F1000Research. 2015;4:1521.
- Patro R, Duggal G, Love MI, Irizarry RA, Kingsford C. Salmon provides fast and bias-aware quantification of transcript expression. Nat Methods. 2017;14:417-419.
- Dobin A, Davis CA, Schlesinger F, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29:15-21.
- Leek JT, Scharpf RB, Bravo HC, et al. Tackling the widespread and critical impact of batch effects in high-throughput data. Nat Rev Genet. 2010;11:733-739.
- Risso D, Ngai J, Speed TP, Dudoit S. Normalization of RNA-seq data using factor analysis of control genes or samples. Nat Biotechnol. 2014;32:896-902.
- Subramanian A, Tamayo P, Mootha VK, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. PNAS. 2005;102:15545-15550.
- Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc B. 1995;57:289-300.
- Ewels PA, Peltzer A, Fillinger S, et al. The nf-core framework for community-curated bioinformatics pipelines. Nat Biotechnol. 2020;38:276-278.