Bio Differential Expression Edger Basics

Version Compatibility

Reference examples tested with: edgeR 4.0+, limma 3.58+, statmod 1.5+ (for voom internals), tximport 1.30+ (for catchSalmon path)

Before using code patterns, verify installed versions match. If versions differ:

• R: packageVersion('<pkg>') then ?function_name to verify parameters

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

edgeR Basics

"Find differentially expressed genes between conditions" -> Fit a negative-binomial GLM per gene with empirical-Bayes-moderated dispersions, test coefficients with the quasi-likelihood F-test (proper finite-sample FPR control), and report ranked DE lists.

The Single Most Important Modern Insight -- edgeR v4 changed the QL framework and legacy=FALSE is now default

Chen, Chen, Lun, Baldoni, Smyth (2025) Nucleic Acids Res 53(2):gkaf018 introduced a bias-corrected adjusted profile likelihood (APL) for dispersion estimation that handles small/zero counts properly, and reworked the QL framework. glmQLFit() in v4 takes legacy=FALSE by default. Old (v3) pipelines that worked perfectly in 2023 now produce subtly different numbers in 2025 -- not wrong, just different. If a tutorial result doesn't match, check which version was used; set legacy=TRUE to reproduce v3 exactly OR (preferred) re-run with v4 defaults and accept the new -- better -- numbers.

Two other v4 changes worth knowing: (1) calcNormFactors() is deprecated in favor of normLibSizes() (same function, new name); (2) method='TMM' remains the documented default; method='TMMwsp' (TMM with singleton pairing) is an alternative introduced for samples with many zeros and is the preferred choice for sparse / low-count data. Pass method= explicitly for reproducibility. Old names still work, so old scripts run, but new code should use the new names.

The choice between Wald-equivalent (glmLRT) and QL-F (glmQLFTest) is not a stylistic preference: glmQLFTest is the modern default because it accounts for uncertainty in the dispersion estimate via an additional QL dispersion. glmLRT is anti-conservative with small n. The two p-values differ -- glmQLFTest p-values are typically >= glmLRT p-values for the same model (the second-stage QL dispersion correction is usually >= 1).

Algorithmic Taxonomy

Decision Tree by Scenario

Standard Workflow

Goal: Take a raw integer count matrix and group labels to a ranked DE table with proper finite-sample FPR control.

Approach: DGEList -> design-aware filtering -> TMM normalization (or TMMwsp for sparse data) -> robust dispersion estimation -> QL fit with robust dispersion shrinkage -> QL F-test on the coefficient of interest -> ranked table.

library(edgeR)
library(limma)

y <- DGEList(counts = counts, group = group)
design <- model.matrix(~ group)

keep <- filterByExpr(y, design)
y <- y[keep, , keep.lib.sizes = FALSE]

y <- normLibSizes(y)
y <- estimateDisp(y, design, robust = TRUE)
plotBCV(y)

fit <- glmQLFit(y, design, robust = TRUE)
qlf <- glmQLFTest(fit, coef = 2)

topTags(qlf, n = 20)
all_de <- topTags(qlf, n = Inf, sort.by = 'none')$table

The robust = TRUE flags propagate Phipson, Lee, Majewski, Alexander, Smyth (2016) Ann Appl Stat 10:946 robust hyperparameter estimation through both the NB dispersion shrinkage (estimateDisp) AND the QL dispersion shrinkage (glmQLFit). Set BOTH; setting only one is the single most common omission. The default robust=FALSE is a relic.

Filtering -- filterByExpr Internals

Goal: Remove genes with insufficient expression to reduce noise and multiple-testing burden, in a design-aware way.

Approach: filterByExpr(y, design) uses the smallest group size from the design matrix to set the minimum number of samples; threshold is CPM >= min.count / median.lib.size * 1e6 AND total count >= min.total.count.

Default parameters: min.count = 10, min.total.count = 15, large.n = 10, min.prop = 0.7.

keep <- filterByExpr(y, design)
y <- y[keep, , keep.lib.sizes = FALSE]

Filter ONCE, before normalization and dispersion estimation. Filtering after estimateDisp invalidates the trended dispersion (the trend was fit on the now-smaller gene set). edgeR has no automatic independent filtering like DESeq2 -- skipping filterByExpr is the single biggest reason an edgeR analysis underperforms a comparable DESeq2 analysis.

Normalization -- TMM/TMMwsp Is an Offset, Not a Division

Goal: Correct for library composition bias (the "few genes consume disproportionate reads" problem) via a per-sample scaling factor.

Approach: normLibSizes(y) defaults to method='TMM' in v4 (same as v3); method='TMMwsp' is the preferred alternative for sparse / single-cell-pseudobulk data with many zeros. Result is a vector of normalization factors stored in y$samples$norm.factors; the effective library size is lib.size * norm.factors. Counts are NOT divided.

y <- normLibSizes(y)
y$samples$norm.factors

cpm_vis <- cpm(y, normalized.lib.sizes = TRUE)
log_cpm <- cpm(y, log = TRUE, prior.count = 2)

The factor enters the GLM as part of the offset. Running TMM on pre-normalized values is wrong (and silent -- gives wrong numbers without error). prior.count = 2 is the modern edgeR default for cpm(log=TRUE); smaller priors (0.25) make low-count log values noisy; larger priors (5-10) shrink them toward zero.

When TMM/TMMwsp assumptions fail (see Failure Modes: prokaryotic stress, MYC amplification, viral host shutoff): use method='upperquartile' or supply known stable reference genes via offsets.

QL F-test vs LRT vs Exact Test

fit_ql <- glmQLFit(y, design, robust = TRUE)
qlf <- glmQLFTest(fit_ql, coef = 2)

fit_lrt <- glmFit(y, design)
lrt <- glmLRT(fit_lrt, coef = 2)

y <- estimateDisp(y)
et <- exactTest(y)

QL F-test is the modern default. LRT is anti-conservative with small n because it does not account for dispersion uncertainty. The exact test handles only two-group one-factor designs and exists for backward compatibility -- use the GLM pipeline for anything with covariates or blocking.

In v4, the QL F-test is bias-corrected via the new APL when legacy=FALSE (the v4 default). To exactly reproduce a v3 result: glmQLFit(y, design, robust = TRUE, legacy = TRUE).

TREAT -- Testing Against a Fold-Change Threshold

Goal: Control FDR for the hypothesis "biologically meaningful fold change", not "fold change non-zero".

Approach: glmTreat(fit, coef=..., lfc=log2(tau)) tests H0: |LFC| <= tau vs HA: |LFC| > tau. The threshold tau MUST be biologically pre-specified -- choosing tau after seeing the data is p-hacking.

tr <- glmTreat(fit, coef = 2, lfc = log2(1.5))
topTags(tr)

The cardinal sin TREAT avoids: filtering padj < 0.05 & abs(logFC) > 1 after a vanilla QL F-test does NOT control FDR for "the LFC exceeds 1". Post-hoc filtering is FDR-controlled for the |LFC|>0 hypothesis only. If reviewers ask "what is the FDR of the '2-fold up' gene list", the only honest answers are TREAT (FDR controlled at the threshold) or "FDR is for non-zero only; the magnitude filter has no FDR guarantee".

Equivalent in DESeq2: results(dds, lfcThreshold = log2(1.5), altHypothesis = 'greaterAbs'). Both implement the McCarthy & Smyth 2009 Bioinformatics 25:765 idea.

Contrasts via No-Intercept Designs

Goal: Define clean pairwise or multi-group contrasts.

Approach: Use ~ 0 + group parameterization so each coefficient is the group mean; build makeContrasts(...) for any pairwise or weighted-average comparison.

design <- model.matrix(~ 0 + group)
colnames(design) <- levels(group)
y <- estimateDisp(y, design, robust = TRUE)
fit <- glmQLFit(y, design, robust = TRUE)

con <- makeContrasts(
    TreatedVsControl = treated - control,
    DrugAVsDrugB     = drugA - drugB,
    ATvsBT           = (treated_A - control_A) - (treated_B - control_B),
    levels = design
)

qlf_t_vs_c <- glmQLFTest(fit, contrast = con[, 'TreatedVsControl'])
qlf_interaction <- glmQLFTest(fit, contrast = con[, 'ATvsBT'])

makeContrasts is more readable than numeric vectors for any non-trivial design.

voom and voomWithQualityWeights

Goal: Use limma's linear-model framework with empirical-Bayes moderation, applying voom precision weights for the mean-variance trend of log-CPM.

Approach: voom transforms counts to log-CPM with per-observation precision weights derived from the mean-variance trend; lmFit + eBayes proceeds as for microarrays.

v <- voom(y, design, plot = TRUE)
fit <- lmFit(v, design)
fit <- eBayes(fit, robust = TRUE)
tt <- topTable(fit, coef = 2, number = Inf)

voomWithQualityWeights (Liu, Holik, Su et al. 2015 NAR 43:e97) adds per-sample weights to downweight outlier samples. Use when:

• Library sizes vary >5x across samples
• RIN varies >2 units
• PCA shows one or more samples clearly off the main cluster

v <- voomWithQualityWeights(y, design, plot = TRUE)
fit <- lmFit(v, design)
fit <- eBayes(fit, robust = TRUE)

eBayes(robust = TRUE) uses Phipson 2016 robust hyperparameter estimation -- NOT the default but strongly recommended for RNA-seq.

Transcript-Level DE via catchSalmon

Goal: Test differential transcript expression with proper inferential variance from quantification bootstrap replicates.

Approach: catchSalmon() / catchKallisto() import per-transcript counts with overdispersion estimates from the bootstrap replicates; pass to DGEList; the rest of the pipeline is standard.

salmon <- catchSalmon(paths = file.path('salmon_out', samples$id))
y_tx <- DGEList(counts = salmon$counts / salmon$annotation$Overdispersion,
                genes  = salmon$annotation)

Baldoni, Chen, Hediyeh-zadeh et al. 2024 NAR 52:e13 ("Dividing out quantification uncertainty"): the per-transcript Overdispersion column captures the variance contribution from the Salmon EM. Dividing the counts by this scaling provides effective counts that limma/edgeR can model as if quantification was certain.