bio-differential-expression-batch-correction

What this skill does

## Version Compatibility

Reference examples tested with: sva 3.50+ (includes ComBat + ComBat_seq), DESeq2 1.42+, edgeR 4.0+, limma 3.58+, RUVSeq 1.36+, ggplot2 3.5+, harmony 1.2+ (single-cell context only)

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.

# Batch Effect Correction

**"Remove the batch effect before DE"** -> Almost always WRONG. Include batch as a covariate in the design formula (`~ batch + condition`) so DESeq2/edgeR/limma model it without subtracting. Subtraction is for visualization only.

## The Single Most Important Modern Insight -- The Nygaard 2016 cardinal sin

Nygaard, Rødland, Hovig 2016 *Biostatistics* 17(1):29-39, "Methods that remove batch effects while retaining group differences may lead to exaggerated confidence in downstream analyses." Translation: **never run ComBat (or ComBat-seq, or `removeBatchEffect`, or SVA-subtract-then-test) and then run DE on the corrected matrix.**

Mechanism: batch-correction methods fit a model `y_ij = alpha + X_ij beta + gamma_i + delta_i epsilon_ij` and subtract the batch terms. The downstream DE test then computes p-values as if those degrees of freedom had never been spent. Residual df is lower than what DESeq2/edgeR/limma assume. Type-I error inflates -- the gene list looks more significant than it should.

The right approach: **include batch in the design**. `~ batch + condition`. DESeq2/edgeR/limma will properly partial out the batch effect from the condition estimate and account for the spent degrees of freedom in the inference. The batch is corrected at the inference stage, not by mutating counts.

`removeBatchEffect` (limma) is for visualization only -- the function's help page says so. ComBat/ComBat-seq output is for visualization, clustering, or downstream tools that cannot take a design matrix (rare; mostly ML).

A second clarification: structural confounding (every "treated" sample in batch 1, every "control" in batch 2) is non-identifiable. No method fixes this -- the "treatment effect" and "batch effect" are mathematically the same vector. The fix is experimental design (randomize batches in advance). ComBat/SVA on a fully-confounded design silently removes the treatment effect along with the batch effect.

## Algorithmic Taxonomy

| Method | Input | Mechanism | Use for |
|--------|-------|-----------|---------|
| Design-matrix inclusion (`~ batch + condition`) | Raw counts; known batch | Partial out batch in the GLM, spent df accounted | DE testing -- the correct path |
| ComBat (Johnson, Li, Rabinovic 2007) | Log-transformed / continuous expression | Empirical-Bayes location and scale shifts per batch | Visualization of microarray / continuous data |
| ComBat-seq (Zhang, Parmigiani, Johnson 2020) | Raw RNA-seq counts | NB-GLM equivalent of ComBat; returns integer counts | Visualization of RNA-seq counts; cross-study harmonization for ML (with caveat) |
| `limma::removeBatchEffect` | Normalized expression | Linear regression subtraction with design protection | Visualization only -- explicit help-page warning |
| SVA (Leek, Storey 2007 + Leek 2012) | Normalized expression | Estimate latent surrogate variables explaining residual variance independent of variable-of-interest | Unknown batch / hidden technical structure -- add SVs to design |
| svaseq (Leek 2014) | Counts | Count-data variant of SVA | Counts version |
| RUVg (Risso, Ngai, Speed, Dudoit 2014) | Counts; negative control genes | Factor analysis of control-gene residual; W as covariate | Strong negative controls (ERCC, housekeeping) -- add W to design |
| RUVs | Counts; replicate samples | Factor analysis within replicate groups; assumes replicates differ only by unwanted variation | Multi-condition with biological replicates as "controls" |
| RUVr | Counts; design only | Factor analysis of residuals from initial fit | Most data-driven; most likely to absorb biology -- caution |
| Harmony, MNN, Scanorama, BBKNN | Single-cell embeddings | Iterative alignment of clusters across samples | Single-cell ONLY; not bulk |

## Decision Tree by Scenario

| Scenario | Recommended approach | Why |
|----------|---------------------|-----|
| Known batch, want DE | `~ batch + condition` in design; do NOT subtract | Cardinal sin avoidance |
| PCA shows batch separation | Include batch in design; for the figure, `removeBatchEffect` is OK (visualization only) | Two purposes, two tools |
| Unknown batch structure | `sva` / `svaseq`; add SVs as covariates in design | Captures latent technical factors |
| Have ERCC spike-ins or trusted housekeeping | `RUVSeq::RUVg` with control gene indices; add W to design | Most principled UV removal |
| Have replicate samples (technical reps within biological) | `RUVSeq::RUVs` | Replicate structure indicates "this differs only by UV" |
| Cross-study integration (TCGA + ICGC + own data) | ComBat-seq for counts, ComBat for log-expression, THEN meta-analysis -- do NOT pool then DE | Goh 2017 warning |
| Fully confounded batch and condition | Re-collect samples; no method fixes this | Non-identifiable |
| Visualizing batch removal for a figure | `removeBatchEffect(expr, batch = batch, design = model.matrix(~condition))` | Visualization is what it's for |
| Single-cell data | Harmony / MNN / Scanorama in the single-cell category, not here | Different problem |

## Standard Workflow -- Design-Matrix Inclusion

**Goal:** Test the condition effect while accounting for known batch structure with proper degrees-of-freedom accounting.

**Approach:** Include batch as a covariate in the design formula; DESeq2/edgeR/limma will partial it out and compute correct p-values.

```r
library(DESeq2)

dds <- DESeqDataSetFromMatrix(countData = counts, colData = coldata,
                               design = ~ batch + condition)
dds <- DESeq(dds)
res <- results(dds, name = 'condition_treated_vs_control')
```

```r
library(edgeR)
y <- DGEList(counts = counts, group = coldata$condition)
keep <- filterByExpr(y, design = model.matrix(~ batch + condition, coldata))
y <- y[keep, , keep.lib.sizes = FALSE]
y <- normLibSizes(y)

design <- model.matrix(~ batch + condition, coldata)
y <- estimateDisp(y, design, robust = TRUE)
fit <- glmQLFit(y, design, robust = TRUE)
qlf <- glmQLFTest(fit, coef = 'conditiontreated')
```

When known confounders are continuous (RIN, library prep date as days), include them as continuous covariates -- no batch correction needed for those:

```r
design = ~ RIN + library_prep_day + condition
```

## ComBat-seq (Visualization or Cross-Study Counts)

**Goal:** Adjust raw counts to remove known batch effects while preserving biology, for VISUALIZATION or for downstream tools that need a single corrected matrix.

**Approach:** `ComBat_seq(counts, batch, group)` returns batch-adjusted integer counts via NB-GLM. Use the output for PCA, clustering, ML -- NOT for DE testing.

```r
library(sva)

corrected_counts <- ComBat_seq(counts = as.matrix(counts),
                                batch = coldata$batch,
                                group = coldata$condition,
                                full_mod = TRUE)

vsd_corrected <- vst(DESeqDataSetFromMatrix(corrected_counts, coldata, ~1))
plotPCA(vsd_corrected, intgroup = 'condition')
```

`full_mod = TRUE` keeps biological covariates protected (the `group` argument). With `full_mod = FALSE`, ComBat-seq removes batch AND group differences -- a serious failure mode.

`ComBat` (the non-seq version) is for log-transformed or microarray data, NOT raw counts. Running `ComBat` on raw counts produces fractional values and assumes Gaussian residuals where counts are NB. Wrong tool, silent failure.

## SVA -- Unknown Batch

**Goal:** Discover latent technical factors when batch is not recorded, add them