Bio Differential Expression Batch Correction

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

Decision Tree by Scenario

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.

library(DESeq2)

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

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:

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.

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 as covariates.

Approach: Estimate the number of surrogate variables; estimate the SVs themselves; add them to the design and re-run DE.

library(sva)
library(DESeq2)

dds <- DESeqDataSetFromMatrix(counts, coldata, design = ~ condition)
dds <- estimateSizeFactors(dds)
norm_counts <- counts(dds, normalized = TRUE)

mod  <- model.matrix(~ condition, coldata)
mod0 <- model.matrix(~ 1, coldata)

n_sv <- num.sv(norm_counts, mod, method = 'leek')
svobj <- svaseq(norm_counts, mod, mod0, n.sv = n_sv)

for (i in seq_len(ncol(svobj$sv))) {
    colData(dds)[[paste0('SV', i)]] <- svobj$sv[, i]
}
sv_formula <- as.formula(paste('~', paste(paste0('SV', seq_len(ncol(svobj$sv))),
                                          collapse = ' + '), '+ condition'))
design(dds) <- sv_formula
dds <- DESeq(dds)

CRITICAL CHECK: compute the correlation between each SV and the variable of interest. If any SV correlates with treatment at r > 0.3, do NOT include it -- doing so would partial out biology, deflating the effect of interest.

sapply(seq_len(ncol(svobj$sv)),
       function(i) cor(svobj$sv[, i], as.numeric(coldata$condition)))

svaseq is the count-data version (Leek 2014); sva is for log-transformed / microarray data.

RUVSeq -- Negative-Control-Anchored UV

Goal: Estimate unwanted variation using ERCC spike-ins, housekeeping genes, replicate samples, or post-fit residuals; add as covariates in design.

Approach: Pick the RUV variant matching the available controls; add the resulting W matrix to the design.

library(RUVSeq)
library(edgeR)

control_idx <- which(rownames(counts) %in% ercc_genes)

set <- newSeqExpressionSet(as.matrix(counts), phenoData = coldata)
set <- betweenLaneNormalization(set, which = 'upper')

ruv <- RUVg(set, control_idx, k = 2)

design <- model.matrix(~ pData(ruv)$W_1 + pData(ruv)$W_2 + condition,
                        data = coldata)
y <- DGEList(counts = counts, group = coldata$condition)
y <- normLibSizes(y)
y <- estimateDisp(y, design, robust = TRUE)
fit <- glmQLFit(y, design, robust = TRUE)
qlf <- glmQLFTest(fit, coef = 'conditiontreated')

Choice of k (number of unwanted factors): no automatic procedure. Try k = 1, 2, 3 and inspect PCA after RUV correction. Treatment effect should remain visible; if it disappears, k is too high. 5-10 is typical; >10 is suspicious.

removeBatchEffect -- Visualization Only

library(limma)

design <- model.matrix(~ condition, coldata)
log_expr_corrected <- removeBatchEffect(log_expr,
                                         batch = coldata$batch,
                                         design = design)

library(ggplot2)
pca <- prcomp(t(log_expr_corrected), scale. = TRUE)
ggplot(data.frame(PC1 = pca$x[,1], PC2 = pca$x[,2],
                  condition = coldata$condition, batch = coldata$batch),
       aes(PC1, PC2, color = condition, shape = batch)) +
    geom_point(size = 3) +
    ggtitle('After removeBatchEffect (visualization only)')

The design argument protects condition while removing batch. Common error: passing batch ALSO in the design ("Coefficients not estimable" warning) -- pass batch via batch= only, NOT also in design=.

DO NOT feed log_expr_corrected back into limma lmFit for DE testing. The function's own help page says so. The right approach: lmFit(log_expr, model.matrix(~ batch + condition, coldata)) -- include batch as a covariate.

Detecting Confounding Early

ct <- table(coldata$condition, coldata$batch)
ct

alias() reveals collinear columns in a design matrix:

alias(model.matrix(~ batch + condition, coldata))$Complete

The Single-Cell Boundary

Harmony (Korsunsky et al. 2019 Nat Methods 16:1289), MNN (Haghverdi et al. 2018 Nat Biotechnol 36:421), Scanorama, BBKNN are designed for single-cell integration where the goal is aligning cell-type clusters across samples/batches. They modify the embedding (PCA coordinates) for downstream UMAP/clustering -- they assume the structural alignment problem of single-cell (many similar cells; align clusters).

DO NOT apply to bulk RNA-seq. Bulk lacks the cluster structure these methods assume.

Conversely, DO NOT use ComBat / ComBat-seq for single-cell. ComBat assumes much more homogeneity than scRNA exhibits; it does not handle zeros well.

See single-cell/batch-integration for the single-cell methods.

The Goh 2017 Critique on "Harmonization"

Goh, Wang, Wong 2017 Trends Biotechnol 35(6):498-507 warn that batch-correction methods can:

• Introduce NEW batch-like artifacts (false positives) when the batch structure is unclear
• Fail when the "most genes not DE" assumption is violated (heat shock, immune activation, viral host shutoff)
• Inflate downstream confidence (echoing Nygaard 2016)

Their recommendation for cross-study work: per-cohort DE first, then meta-analyze the effect sizes (metafor, limma treat across cohorts). DO NOT pool data, harmonize, then re-test.

"Harmonization" in clinical genomics often means ComBat across cohorts then DE on the harmonized matrix -- the exact failure mode Nygaard 2016 describes. Same problem, different vocabulary.