Skip to content
Bench to Insights
Blog

edgeR and limma

DESeq2 is one of three standard routes from a count matrix to a differential expression table. The other two are edgeR and limma-voom, and they model the counts differently. edgeR keeps the negative binomial model that DESeq2 uses. limma-voom moves the counts onto a log scale and fits a linear model instead, which runs faster and handles complex designs more easily. The three agree on the genes that change and disagree at the margin, so a second fit is a check on the first rather than a waste.

edgeR

Section titled “edgeR”

The edgeR pipeline mirrors the DESeq2 one. Counts go into a DGEList, genes with low expression are filtered out, and the trimmed mean of M-values handles the normalization. edgeR then estimates the dispersion and fits a generalized linear model, and the comparison of interest is the treatment coefficient.

library(edgeR)
library(airway)


data(airway)
counts <- assay(airway)
group <- colData(airway)$dex


# Build the DGEList and filter low-expression genes.
y <- DGEList(counts = counts, group = group)
keep <- rowSums(cpm(y) > 1) >= 2
y <- y[keep, , keep.lib.sizes = FALSE]


# TMM normalization and dispersion estimation.
y <- calcNormFactors(y, method = "TMM")
design <- model.matrix(~ group)
y <- estimateDisp(y, design)


# Fit the model and test the treatment coefficient.
fit <- glmFit(y, design)
lrt <- glmLRT(fit, coef = 2)
head(topTags(lrt, n = 5)$table)
                    logFC   logCPM       LR       PValue          FDR
ENSG00000152583 -4.583391 5.539109 360.8828 1.808617e-80 2.721969e-76
ENSG00000179094 -3.166942 5.180005 207.8318 4.082315e-47 3.071942e-43
ENSG00000125148 -2.186302 7.415214 177.7747 1.483563e-40 7.442543e-37
ENSG00000120129 -2.930761 7.312678 173.7642 1.114491e-39 4.193273e-36
ENSG00000189221 -3.288546 6.770899 169.3914 1.004819e-38 3.024505e-35

The cpm(y) > 1 filter keeps a gene when at least two samples show at least one count per million, the edgeR analogue of the read-count floor in DESeq2. Shrinkage enters at the dispersion step, where empirical Bayes pulls the gene-wise estimates toward a trend, the same idea as in DESeq2, fitted by a different procedure.

limma-voom

Section titled “limma-voom”

limma was written for microarrays, and voom extends it to count data. The transform estimates the mean-variance relationship of the log-counts and turns it into a precision weight for each observation, so a weighted linear model replaces the count model. Because the model is linear, the fit is fast and it scales to complex designs and large studies.

library(limma)


# Voom transform the TMM-normalized counts.
voom_y <- voom(y, design, plot = FALSE)


# Fit the linear model and apply empirical Bayes.
fit_lm <- lmFit(voom_y, design)
fit_lm <- eBayes(fit_lm)
head(topTable(fit_lm, coef = 2, number = 5))
                    logFC  AveExpr         t      P.Value    adj.P.Val        B
ENSG00000134686 -1.372871 6.840212 -17.14849 4.826557e-08 0.0002233334 9.096686
ENSG00000152583 -4.566806 4.168205 -18.54747 2.469027e-08 0.0002233334 8.845068
ENSG00000125148 -2.189646 7.025373 -16.30438 7.419715e-08 0.0002233334 8.691968
ENSG00000179094 -3.177256 4.421723 -16.69263 6.072728e-08 0.0002233334 8.626749
ENSG00000148175 -1.431891 8.857331 -16.55524 6.515572e-08 0.0002233334 8.596265

The voom call reuses the DGEList from the edgeR step, so the normalization carries over. eBayes then shrinks the gene-wise variances toward a common estimate, which is what keeps limma stable on small sample sizes.

When to choose which

Section titled “When to choose which”
Package Model Use it when
DESeq2 negative binomial the standard choice, well-documented defaults
edgeR negative binomial you want control over the dispersion model
limma-voom linear model on log-counts the design is complex, or the study is large

In practice the three packages return overlapping gene lists. The two tables on this page share ENSG00000152583, ENSG00000179094 and ENSG00000125148 near the top, even where the ranking differs. The disagreement is near the significance threshold, and a gene that passes in two packages out of three is a safer call than one that passes in a single fit. The enrichment page takes the DESeq2 result forward, and the visualization page plots it.