Batch Effect Correction in Omics Data: ComBat, ComBat-Seq and When Not to Correct

You prepare samples in two batches, six months apart, and the PCA separates perfectly by batch. Now every gene you care about correlates with the batch variable, and you cannot tell whether that is biology or the reagent lot. This is the most common genuinely hard problem in omics analysis, and the standard fixes range from correct to actively harmful.

The single most important idea in this tutorial is that batch is a property of your experimental design, not just your data. What you can legitimately do about it is largely determined before you ever open R.

What a batch effect actually is

A batch effect is variation that enters systematically with a technical grouping rather than with the biology you meant to study. Common sources: library preparation date, sequencing run, operator, reagent lot, culture incubator, sample storage time, and the machine a proteomics plate was run on.

It is not random noise, which is why it does not average out. It shifts whole groups of samples in a consistent direction, so it can dominate the first principal component and produce convincing-looking differential expression that is purely technical.

Step one: look, before you fix

library(SummarizedExperiment)
library(ggplot2)

# Variance-stabilise counts; raw counts are dominated by depth differences
vsd <- vst(dds, blind = FALSE)

pca <- plotPCA(vsd, intgroup = c("condition", "batch"), returnData = TRUE)
pct  <- round(100 * attr(pca, "percentVar"))

ggplot(pca, aes(PC1, PC2, colour = condition, shape = batch)) +
  geom_point(size = 3) +
  xlab(paste0("PC1: ", pct[1], "% variance")) +
  ylab(paste0("PC2: ", pct[2], "% variance")) +
  theme_bw()

Read the plot before deciding anything. If samples separate by condition within each batch, batch is a nuisance and small. If PC1 tracks batch, you have a real problem. If batch and condition are perfectly aligned — all controls in batch 1, all treated in batch 2 — you have no problem in the statistical sense, you have an unfixable design, and no amount of correction will recover it.

Quantify the confound

table(colData(dds)$condition, colData(dds)$batch)

# Is any condition present in only one batch?
xt <- table(colData(dds)$condition, colData(dds)$batch)
which(xt == rowSums(xt), arr.ind = TRUE)

Stop here if the table shows a condition that exists in exactly one batch. Correction algorithms assume they can separate the biological signal from the technical one; when those are perfectly confounded, the separation is mathematically impossible and any ‘successful’ correction is generating an answer from an assumption.

The preferred fix: put batch in the model

For differential expression, the correct answer is almost always to include batch as a covariate rather than modifying the data. You estimate the biological effect while accounting for batch, and you never invent expression values.

design(dds) <- ~ batch + condition
dds <- DESeq(dds)

res <- results(dds, contrast = c("condition", "treated", "control"))
summary(res)

This works when batch is a nuisance but not the main story, and it is statistically cleaner than transforming the data because it propagates uncertainty properly. Limma follows the same logic:

design <- model.matrix(~ batch + condition)
fit <- lmFit(expr, design)
fit <- eBayes(fit)
topTable(fit, coef = "conditiontreated", number = 20)

The order of terms matters in a sense: the biological condition should be the term you test, and batch should be a blocking variable that gets estimated but not tested. R’s model formula handles this, but you must extract the right coefficient.

When you do need to correct the data

Some analyses cannot take covariates: clustering, heatmaps, PCA, machine learning, sample-level distance matrices. All of those operate on the data matrix itself, so batch has to leave the matrix or it will dominate the structure you are looking at.

Note the asymmetry: correction for visualisation and unsupervised learning is legitimate. Correction before differential expression is usually a mistake, because you then fit a model that no longer needs the batch term and you have thrown away the uncertainty the correction introduced.

ComBat for transformed, approximately normal data

library(sva)

# mod must describe the biology you want to preserve
mod <- model.matrix(~ condition, data = colData(dds))

expr <- assay(vsd)
expr_corrected <- ComBat(
  dat   = expr,
  batch = colData(dds)$batch,
  mod   = mod,
  par.prior = TRUE,
  prior.plots = FALSE
)

The mod argument is not optional bookkeeping — it tells ComBat which variation to protect. Leave it out and ComBat will happily remove the difference between your conditions along with the batch, especially when the design is unbalanced. Always pass the biological model explicitly and inspect the PCA afterwards.

ComBat-Seq for raw counts

library(sva)

counts_corrected <- ComBat_seq(
  counts  = counts(dds),
  batch   = colData(dds)$batch,
  group   = colData(dds)$condition,
  full_mod = TRUE
)

Use ComBat_seq, not ComBat, on raw counts. Plain ComBat expects approximately Gaussian data and will produce negative values on counts — a clear sign you have applied the wrong method, and one people frequently work around with a log transform rather than switching tools.

removeBatchEffect: the one that must never touch a test

# For a PCA or heatmap only
expr_vis <- removeBatchEffect(expr, batch = colData(dds)$batch,
                              design = mod)

plotPCA(expr_vis, intgroup = "condition")

removeBatchEffect subtracts estimated batch coefficients and returns a matrix with no associated uncertainty. Running a t-test or a linear model on its output produces p-values that are too small and entirely fictional. Its documentation says this; the warning is widely ignored.

Overcorrection: the failure that looks like success

An overcorrected dataset has a beautiful PCA — samples mixed perfectly by batch — and no biology left. Signs to watch for:

• Differential expression results shrink dramatically after correction, especially for genes you know are real.
• The corrected data no longer cluster by any biological grouping you trust.
• Batch variance is reduced to near zero rather than to the level of other technical noise.
• A known positive control gene loses its expected effect.

The safeguard is a positive control set: genes or proteins whose behaviour you already know from the literature or from spike-ins. Check that they survive correction. If they do not, your correction is removing signal, not noise.

A sanity check worth running every time

# Compare known markers before and after
markers <- c("GAPDH", "ACTB", "IL6", "TNF")
rbind(
  before = rowMeans(expr[markers, ]),
  after  = rowMeans(expr_corrected[markers, ])
)

# And confirm the expected biological separation survived
plotPCA(expr_corrected, intgroup = "condition")

Design decisions that prevent the problem

If a study is already collected and badly confounded, the honest move is to report the confound and interpret conservatively — not to run a correction and present the output as clean. A reviewer who notices perfect confounding will reject the finding, and they will be right to.

What to write in the methods

State the batch structure, whether batch was included in the model or corrected in the data, the method and version, the biological covariates protected during correction, and what you used corrected values for. A sentence like batch was included as a covariate in the DESeq2 design; ComBat-Seq was applied only for the clustering shown in Figure 2 tells a reader exactly how much to trust each panel.

That clarity — separating what the statistics tested from what the figure displays — is the difference between a defensible analysis and one that merely looks tidy. If you want to build that habit across a whole study, from design through model to figure, Ampersand Academy works through the full chain.

Practice

1. Make a PCA coloured by batch and by condition for a public dataset and describe which dominates PC1.
2. Build the condition-by-batch table and identify whether any condition is confounded.
3. Run DESeq2 with and without a batch term and compare the number of significant genes.
4. Apply ComBat with and without mod and see how the results differ — then explain which is defensible.

Related tutorials

• Bulk RNA-Seq Data Analysis: From FASTQ to a Trustworthy Differential-Expression Result — the wider RNA-seq workflow batch effects appear in
• Volcano Plots and Heatmaps for RNA-Seq Results: A Practical Tutorial — where an uncorrected batch effect becomes visible
• R for Bioinformatics: Bioconductor, DESeq2, and What R Is Actually For — DESeq2 and Bioconductor, the tools these corrections run alongside