bio-differential-expression-edger-basics

What this skill does

## 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

| Test | What it does | When mandatory | Failure mode |
|------|--------------|----------------|--------------|
| QL F-test (`glmQLFit` + `glmQLFTest`) | NB GLM with second-stage QL dispersion to model dispersion uncertainty | DEFAULT for any modern bulk DE | None within design assumptions |
| LRT (`glmFit` + `glmLRT`) | Likelihood-ratio of nested GLMs using point dispersion estimate | Only when n>=10/group AND simple design AND no QL F-test available | Anti-conservative with small n -- inflated false positives |
| Exact test (`exactTest`) | Conditional NB test for two groups, no covariates | Legacy; one-factor two-group ONLY | Cannot adjust for batch or covariates |
| TREAT (`glmTreat`) | Tests H0: |LFC| <= tau vs HA: |LFC| > tau (McCarthy & Smyth 2009) | Want FDR control for "biologically meaningful fold change" claim | Conservative; tau must be pre-specified |
| voom + lmFit + eBayes | Linear model with empirical-Bayes moderation on voom-weighted log-CPM | Heterogeneous library sizes (>3x); samples with widely varying quality | Assumes log-normal-after-weighting; less direct count model |
| voomWithQualityWeights | voom + per-sample quality weights from arrayWeights (Liu 2015) | Some samples markedly worse quality (low RIN, contamination) | With very small n, sample weights are noisy |

## Decision Tree by Scenario

| Scenario | Recommended path | Why |
|----------|------------------|-----|
| Standard bulk RNA-seq, n>=3/group, modern script | `filterByExpr -> normLibSizes -> estimateDisp(robust=TRUE) -> glmQLFit(robust=TRUE) -> glmQLFTest` | Modern default with proper FPR control |
| n=2-3/group | edgeR QL F-test (over DESeq2) | Schurch 2016 *RNA* 22:839 + edgeR v4 changes: QL is the tightest FPR control at small n |
| Library sizes vary >3x or RIN heterogeneous | `voom` or `voomWithQualityWeights` + lmFit + eBayes | Precision weights downweight noisy observations per sample |
| Pre-specified biologically meaningful fold-change threshold | `glmTreat(fit, coef=..., lfc=log2(1.5))` | Post-hoc filtering by |LFC| does NOT control FDR for the magnitude hypothesis |
| Multi-group, "any change" omnibus | `glmQLFTest(fit, coef=2:k)` | Joint F-test on multiple coefficients |
| Multi-group, all pairwise contrasts | `~ 0 + group` parameterization + `makeContrasts` | Clean, readable contrasts |
| Transcript-level DE (DTE) | `catchSalmon()` / `catchKallisto()` -> `DGEList` with overdispersion | Baldoni 2024 *NAR* 52:e13: properly handles inferential variance |
| Cross-tool sanity check | Run DESeq2 in parallel; expect >=70% overlap at top 500 | <60% overlap suggests modeling problem, not tool difference |
| Single-cell pseudobulk | Aggregate counts per donor; standard edgeR QL pipeline | Crowell 2020 *Nat Commun* 11:6077: pseudobulk avoids the FDR inflation of cell-level DE |
| Reproducing pre-2025 result | `glmQLFit(legacy=TRUE)` | edgeR v4 introduced bias-corrected APL with `legacy=FALSE` default; v3 was slightly biased |

## 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.

```r
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`.

```r
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.

```r
y <- normLibSizes