Before running a single line of R, understand why count data needs special statistics — and when to reach for edgeR, limma-voom, or DESeq2.
When you sequence a sorghum transcriptome, the sequencer produces millions of short reads. After alignment with STAR or HISAT2, you count how many reads overlap each gene. The result is a count matrix — genes as rows, samples as columns, integers as values.
This sounds simple, but the numbers have a peculiar statistical character that breaks classical methods:
In a drought-stress experiment comparing well-watered vs water-limited sorghum, you might have 30,000 genes but only 3 biological replicates per condition. Classical statistics assume large sample sizes to estimate variance accurately. With n=3, you need to borrow information across genes to make reliable estimates — this is exactly what edgeR and limma do differently from a simple t-test.
A raw count of 50 reads means something very different in a library of 1 million reads vs 50 million reads. This sequencing depth variation must be accounted for during normalisation, before any statistical test is applied.
The classic Student's t-test assumes that your measurements follow a normal (Gaussian) distribution. RNA-seq counts violate this assumption in three important ways:
| Assumption of t-test | Reality of RNA-seq counts | Consequence of ignoring it |
|---|---|---|
| Data is continuous | Counts are discrete integers (0, 1, 2 …) | Invalid probability calculations |
| Variance is constant (homoscedastic) | Variance increases with the mean | High-count genes dominate; low-count genes buried |
| Data is normally distributed | Heavily right-skewed; many zeros | Inflated false positives for lowly expressed genes |
| Large sample size | Typically 2–5 replicates per group | Poor variance estimation → unreliable p-values |
Historical note: Early microarray papers used t-tests and fold-change thresholds on log-transformed data. This worked reasonably well because microarray intensities are closer to normally distributed. The RNA-seq community learned the hard way — by comparing results — that importing microarray statistics into count data analysis produces poor results.
Some researchers attempt to log-transform counts and apply a t-test. While this fixes the skewness partially, it introduces a new problem: log(0) is undefined. Adding a pseudo-count (e.g., log(count + 1)) partially helps but distorts low-count genes. The dedicated methods described below handle this properly.
The distribution that best describes RNA-seq count data is the negative binomial (NB) distribution. It generalises the simpler Poisson distribution by adding an extra parameter — the dispersion — that captures biological variability between samples.
Imagine a gene that produces exactly 100 mRNA molecules per cell. If you count 1,000 cells, the Poisson distribution describes the variation you see (it equals the mean: variance = 100). But cells in a real sorghum leaf are not identical — some cells respond more strongly to drought stress than others. This extra biological variability is captured by the dispersion parameter φ: Var(Y) = μ + φμ². When φ=0, you get back to Poisson. edgeR and DESeq2 both use this NB model. The difference is in how they estimate φ.
The NB distribution has two parameters:
Estimating φ reliably is the central challenge. With only 3 replicates, you cannot estimate gene-specific dispersion accurately. The three major tools (edgeR, DESeq2, limma-voom) each take a different approach to this problem, and understanding those differences tells you when to use which tool.
# The NB variance formula — the key equation to understand # # Var(Y) = mu + phi * mu^2 # # mu = mean count for this gene in this condition # phi = dispersion (biological + technical extra-variation) # # When phi = 0 → Var = mu (Poisson — no overdispersion) # When phi > 0 → Var > mu (negative binomial — realistic) # # Example: a gene with mu = 100 reads mu <- 100 # Poisson variance var_poisson <- mu cat("Poisson variance:", var_poisson, "\n") # Poisson variance: 100 # Negative binomial variance with phi = 0.1 (10% dispersion) phi <- 0.1 var_nb <- mu + phi * mu^2 cat("NB variance (phi=0.1):", var_nb, "\n") # NB variance (phi=0.1): 1100 # The NB standard deviation is sqrt(1100) ≈ 33 # Much larger than the Poisson sd of sqrt(100) = 10 # This is why you need biological replicates — and proper statistics!
Key takeaway: The NB dispersion parameter φ is not a nuisance — it is biologically meaningful. Genes with high dispersion (e.g., stress-response genes in sorghum) are inherently noisier and need stronger evidence before you call them differentially expressed. edgeR and DESeq2 use this to penalise noisy genes automatically.
One of the most important diagnostic plots in RNA-seq is the mean–variance plot. It reveals how variance scales with the mean across all genes in your dataset.
In RNA-seq data you will always observe that highly expressed genes also have higher variance. This is expected from the NB model: Var = μ + φμ². The challenge is that this trend is not perfectly smooth — some genes deviate from the trend, especially at low counts.
edgeR fits the mean–variance trend to estimate dispersion using a weighted likelihood approach. It then shrinks gene-level dispersion estimates towards a common trend (empirical Bayes shrinkage), borrowing strength across genes.
limma-voom takes a different path: the voom() function computes precision weights for each observation (each count in each sample) based on where it falls on the mean–variance trend. These weights are then passed to limma's linear model framework, which was designed for microarray data. The weights convert the heteroscedastic count data into something that behaves like normally distributed data.
DESeq2 uses a slightly different parametrisation of the NB model and a different shrinkage estimator (the regularised log transform, rlog, or variance-stabilising transformation, VST). All three are valid — the choice depends on your experimental design and the size of your dataset.
The mean–variance trend also explains why normalisation must happen before DE testing. If Library A has 30M reads and Library B has 15M reads, every gene in Library A will have roughly double the counts — even if expression is identical. Normalisation methods (TMM in edgeR, median-of-ratios in DESeq2) correct for this before any statistical test.
edgeR, limma-voom, and DESeq2 are the three dominant tools for bulk RNA-seq differential expression. They are not competitors — they are complementary approaches developed by the same bioinformatics community (mostly from the Gordon Smyth group at WEHI, Melbourne, and the Love/Anders/Huber group). Understanding their differences makes you a better analyst.
| Feature | edgeR | limma-voom | DESeq2 |
|---|---|---|---|
| Statistical model | Negative binomial GLM | Linear model + precision weights (voom) | Negative binomial GLM |
| Dispersion estimation | Weighted conditional likelihood + empirical Bayes | Mean–variance trend weights (no explicit dispersion) | Local regression shrinkage (DESeq2 shrinkage) |
| Normalisation | TMM (Trimmed Mean of M-values) | TMM (via edgeR) or other | Median-of-ratios |
| Fold change shrinkage | Not by default (optional) | Not typically | Yes (apeglm, ashr, normal) |
| Best for small n | ✅ Yes (n=2 possible) | ✅ Yes | ✅ Yes, with shrinkage |
| Complex designs | ✅ Full GLM flexibility | ✅ Excellent (uses limma's lmFit) | ⚠️ Possible but less flexible |
| Gene set testing | ✅ fry, camera, goana, kegga built-in | ✅ Same (limma-based) | ⚠️ External packages needed |
| R package | edgeR (Bioconductor) |
limma (Bioconductor) |
DESeq2 (Bioconductor) |
Common misconception: Many beginners think you must pick one tool and stick with it. Professional bioinformaticians often run two or three tools and look at the intersection of their results. Genes that are significant in edgeR, limma-voom, AND DESeq2 are the most reliable hits. We cover this Venn diagram approach in Lesson 9.
edgeR is particularly well-suited to:
glmQLFTest), even n=2 gives reasonable results.glmQLFTest handles complex interaction terms, covariates (batch, sex, treatment × time), and contrasts between any combination of groups.camera, fry, goana, and kegga functions in edgeR/limma are among the best-implemented gene set tests in Bioconductor.In your BPV genomic selection work, you may have RNA-seq data from drought-stressed vs well-watered sorghum lines, possibly with multiple time points (seedling, tillering, anthesis). This is a multi-factor design — genotype × treatment × time. edgeR's GLM framework handles this naturally, letting you test: "Which genes respond to drought specifically at anthesis, after accounting for genotype effects?"
limma-voom is often the preferred choice for:
lmFit approach is extremely efficient for designs with many covariates (batch effects, continuous variables like flowering time, random effects modelled with duplicateCorrelation).voomWithQualityWeights detects outlier samples and automatically down-weights them — very useful when one replicate has a degraded library.The voom insight (Law et al., 2014): The genius of limma-voom is that it converts the heteroscedastic count problem into a problem limma already knows how to solve. By computing precision weights from the mean–variance trend and passing them to lmFit, it lets you use 20 years of linear model theory for RNA-seq. This makes it extremely flexible for unusual designs.
In this module, you will learn both edgeR and limma-voom thoroughly. By Lesson 9, you will be able to run both on the same dataset, compare their results, and explain to a collaborator which genes you trust the most and why.
These conceptual exercises test your understanding before writing any code. Think through each one carefully — the answers reinforce the statistical foundations that every subsequent lesson builds on.
You have a gene with mean count μ = 200 in the control group. You observe counts across 4 replicates: 185, 212, 193, 220. Calculate the observed variance. Is this consistent with a Poisson distribution (where variance = mean = 200), or does it suggest overdispersion?
Mean: (185 + 212 + 193 + 220) / 4 = 202.5 (close to 200, good)
Variance: [(185−202.5)² + (212−202.5)² + (193−202.5)² + (220−202.5)²] / 3 = [306.25 + 90.25 + 90.25 + 306.25] / 3 = 264.3
The Poisson distribution predicts variance = mean = 200. The observed variance is 264.3 — larger than the mean. This is overdispersion, consistent with the negative binomial model. With biological replicates, overdispersion is virtually always present in real RNA-seq data.
A colleague runs a t-test on raw counts to find DE genes in a sorghum drought experiment. They find 8,000 significant genes at FDR < 0.05. List three reasons why this result is likely unreliable.
1. Normality assumption violated: RNA-seq counts are discrete, right-skewed, and zero-inflated — far from normal. The t-test's p-values assume normality; violations inflate false positives especially for low-count genes.
2. Heteroscedasticity ignored: The t-test assumes equal (or at least constant) variance. RNA-seq variance grows with the mean. Highly expressed genes will have their variance underestimated relative to lowly expressed genes, producing biased results in both directions.
3. No library size normalisation: Raw counts are not normalised for sequencing depth. A gene with 100 counts in a 10M-read library and 100 counts in a 50M-read library has very different expression levels. Without TMM or median-of-ratios normalisation, the t-test cannot distinguish real expression differences from depth differences. 8,000 "significant" genes almost certainly includes thousands of false positives from library size effects alone.
For each scenario below, which tool would you use first — edgeR, limma-voom, or DESeq2? Justify your choice.
A) Sorghum experiment, n=2 per group (drought vs control), no budget for more sequencing.
B) Human clinical trial, n=40 per group, with continuous covariates (age, BMI, batch date).
C) A standard sorghum experiment, n=4 per group, and you want fold change shrinkage for a publication-quality MA plot.
A) edgeR — specifically the classic exact test or glmQLFTest. edgeR was designed for very small n and can even be used with n=2. DESeq2 and limma-voom can also work, but edgeR has the best-established track record with n=2.
B) limma-voom — with n=40, the normal approximation is excellent. limma's linear model framework handles many continuous covariates naturally via lmFit, and voomWithQualityWeights can flag outlier samples in a large clinical cohort.
C) DESeq2 — DESeq2's lfcShrink() with the apeglm method is specifically designed for shrinking log fold changes toward zero for lowly expressed or variable genes. This produces cleaner MA plots and more conservative, reproducible results for publication. (You might also run edgeR or limma-voom and compare the intersection.)
Gene A has dispersion φ = 0.01. Gene B has dispersion φ = 0.5. Both have mean expression μ = 100. Which gene requires stronger evidence (larger fold change, more replicates, or lower p-value threshold) to be called differentially expressed? Why?
Gene B (φ = 0.5) requires stronger evidence.
Var(Gene A) = 100 + 0.01 × 100² = 100 + 100 = 200. SD ≈ 14.1
Var(Gene B) = 100 + 0.5 × 100² = 100 + 5000 = 5100. SD ≈ 71.4
Gene B is five times noisier. A fold change of 2× in Gene B could easily be explained by biological noise between replicates. The same fold change in Gene A, which has low dispersion, is much more likely to represent a true biological difference. edgeR and DESeq2 automatically penalise high-dispersion genes, producing wider confidence intervals and higher p-values for noisy genes — even when the observed fold change is the same.