How do I check for batch effects in bulk RNA-seq?

Start by cross-tabulating your batch variables (prep date, flowcell, lane, kit lot) against condition to check the design isn't confounded, then run PCA and a sample-distance heatmap on the variance-stabilized matrix colored by both batch and condition. If PC1 tracks a batch variable instead of condition, or the table is block-diagonal, you have a problem to address before trusting any DE result.

Should I run ComBat before DESeq2?

No, not for differential expression testing. Pre-correcting the count matrix with ComBat or ComBat-seq and then feeding it into DESeq2 or edgeR causes incorrect degrees of freedom and artificially small p-values, because the model doesn't know the counts have already been adjusted. Put batch in the design formula instead (~ batch + condition) and let DESeq2 or edgeR subtract the effect during testing.

What's the difference between removeBatchEffect() and adding batch to the design formula?

removeBatchEffect() edits the expression matrix itself and is meant only for visualization: PCA plots, heatmaps, clustering figures. Adding batch to the design formula (~ batch + condition) leaves the raw counts untouched and instead has the negative binomial or linear model estimate a batch coefficient and subtract it during hypothesis testing, which is the statistically correct approach for anything that produces a p-value.

Can batch and condition be confounded in a way I can't fix?

Yes. If every sample in one batch belongs to only one condition and every sample in another batch belongs to the other condition, the design matrix is not full rank and no statistical method, including ComBat-seq or a design-formula batch term, can separate batch from biology. The only real fix is generating new samples that break the confound.

When should I not correct for a batch effect?

When what looks like a batch is actually a biological replicate structure (different patients, animals, or cohorts) or a real subtype difference rather than a technical artifact. Correcting it out in those cases removes genuine biological variability you may need to report, not noise.

Sanity check · Bulk RNA-seq

How to Detect Batch Effects in Bulk RNA-seq

A PCA plot that separates by prep date instead of treatment is telling you the truth about your experiment, not a bug to plot around.

By Ming "Tommy" Tang, Director of Bioinformatics in Big Pharma · Updated 2026-09-13 · 6 min read

You ran DESeq2, got a gene list, and then someone asked you to color the PCA plot by processing date. PC1 splits cleanly by date, not by control versus treatment. Now you're staring at the plot wondering whether the last three days of work are salvageable, and whether removeBatchEffect() or ComBat is the thing that fixes it.

What's at stake isn't cosmetic. If processing date, library prep batch, flowcell, or operator lines up with your biological groups, "condition differences" in your results table are really batch plus biology mixed together, and no downstream method can pull them apart after the fact. Leek and colleagues showed batch can explain more variance than the biology you're studying, and it can flip the sign of an apparent effect. A gene list built on a confounded design is not wrong in a way you can catch by filtering harder; it's wrong at the design level.

This page gives you the checks to run before you trust a PCA plot or a DE table, in order from a five-minute metadata audit to a proper variance-partition analysis, plus the design-formula fix that actually works versus the pre-correction move that quietly breaks your p-values.

What it looks like when it's happening

PC1 or PC2 in a PCA of the variance-stabilized matrix separates samples by processing date, sequencing lane, or operator rather than by condition
A sample-to-sample distance heatmap clusters replicates by which day they were prepped, not by treatment group
table(batch, condition) is block-diagonal: each batch contains only one condition, with empty off-diagonal cells
Two genes with no known relationship correlate tightly across samples, but only because both track RNA quality or a lane effect
The direction of a fold change flips, or a gene drops out of the results entirely, once batch is added to the design formula
Replicates from the same prep day sit closer together in PCA space than replicates from the same condition
PC1's loading correlates more strongly with RIN, %rRNA, or total library size than with any biological covariate you recorded

Why it happens

Every RNA-seq library goes through extraction, rRNA depletion or poly-A selection, fragmentation, adapter ligation, and amplification, then gets sequenced on a specific flowcell and lane. Each of those steps has lot-to-lot and day-to-day variability: a different kit lot changes rRNA depletion efficiency, a different operator changes fragmentation size distribution, a different flowcell changes base-call quality and GC bias. None of that is biology, but all of it shifts read counts across thousands of genes simultaneously, in a coordinated direction. PCA is a linear decomposition of total variance, so whichever source of variation moves the most genes together captures the top axes. With typical bulk designs of three to six replicates and modest true fold changes, a systematic technical shift across an entire batch is often larger in magnitude than the biological signal, so it wins PC1.

The design-level failure is worse than a noisy covariate. If every control sample was prepped in batch 1 and every treated sample in batch 2, the two terms are perfectly collinear: the design matrix used by DESeq2, edgeR, or limma cannot assign variance to "batch" and "condition" separately, because every batch-1 sample is also a control and vice versa. The GLM will fit a coefficient for "condition," but that coefficient is mathematically indistinguishable from "batch." This is not a statistical power problem you can fix with more genes or a different test; it's a rank problem in the design.

Batch effects also distort structure you'd trust for other reasons. A sample correlation heatmap or a co-expression network built on raw or naively normalized counts will show genes as "co-regulated" simply because both respond to the same technical artifact, like RNA degradation, across the same set of samples. And batch frequently correlates with sequencing depth, because samples prepped on the same day are often sequenced together: depth imbalance makes lowly-expressed genes look "up-regulated" in the deeper batch purely from better sampling, which stacks a second artifact on top of the first.

Correction methods can only work when batch and condition are not perfectly confounded, because they rely on seeing multiple conditions within each batch to estimate what "batch" looks like independent of biology. When that within-batch variation exists, DESeq2, edgeR, and limma all model batch as an additive term in the GLM and subtract its effect during hypothesis testing, without ever touching the raw counts. When it doesn't exist, the honest answer is that the experiment cannot answer the question it was designed to answer.

The checks

Run them in order. Each one tells you what healthy looks like and what the problem looks like.

0/8 checked · saved in this browser

01Audit recorded batch metadata
Before touching the count matrix, list every sample against every plausible batch variable: processing date, RNA extraction date, library prep kit lot, rRNA-depletion or poly-A kit lot, sequencing flowcell, lane, and operator. If any of these weren't recorded at the bench, flag it now; you cannot reconstruct it later from the counts.
Healthy
A complete colData table where every sample has a value for every batch variable, alongside condition.
Red flag
Batch metadata is missing, was reconstructed from memory after the fact, or only 'processing date' was kept while lane and kit lot were not.
02Cross-tabulate batch against condition
Run a contingency table of batch versus condition and look at whether every batch contains a mix of conditions.
r
```
table(coldata$batch, coldata$condition)
```
Healthy
Off-diagonal cells are populated: each batch has samples from more than one condition, so the design matrix is full rank and batch can be estimated separately from condition.
Red flag
The table is block-diagonal, i.e. each batch row has nonzero counts in only one condition column. This is complete confounding, not a modeling problem.
03Confirm batch is encoded as a factor, not numeric
Check the class of the batch column in colData. If it was read in as integers (1, 2, 3), the model will treat batch as a continuous ordinal variable, implying batch 3 is 'twice' batch 2 rather than three unrelated categories.
r
```
coldata$batch <- factor(coldata$batch)
class(coldata$batch)
```
Healthy
class(coldata$batch) returns 'factor', with levels corresponding to discrete batches.
Red flag
class(coldata$batch) returns 'numeric' or 'integer'.
04PCA colored by batch and by condition, side by side
On the variance-stabilized or rlog-transformed matrix, generate two PCA plots from the same principal components: one colored by each batch variable (prep date, flowcell, lane, kit lot, operator) and one colored by condition. Note the percent variance explained by PC1 and PC2 in each.
Healthy
PC1 and PC2 track condition; batch variables show no visible separation, or only appear on lower, minor-variance PCs.
Red flag
PC1 or PC2 cleanly separates by a batch variable instead of, or in addition to, condition.
05Sample-to-sample distance heatmap annotated with batch and condition
Compute pairwise sample distances on the same transformed matrix and cluster them, annotating rows and columns with both batch and condition.
Healthy
Samples cluster primarily by condition; any residual sub-clustering by batch is a minor effect nested within condition blocks.
Red flag
The dendrogram splits first by batch, with condition only forming sub-clusters inside each batch block.
06Correlate PC1/PC2 with library size and quality metrics
Numerically correlate PC1 and PC2 scores against total library size, number of detected genes, RIN, and %rRNA or %mitochondrial content, in addition to batch variables.
Healthy
Low correlation between PC1/PC2 and these technical metrics; the top components track condition, not depth or RNA quality.
Red flag
PC1 correlates strongly with library size or RIN rather than with any biological covariate you recorded.
07Quantify variance attributable to batch vs condition per gene
Use a quantitative tool such as variancePartition or BatchQC to fit a model with both batch and condition as terms, and extract the fraction of variance each explains per gene, rather than relying on PCA alone.
Healthy
Condition explains a meaningful share of variance for a plausible number of genes; batch's contribution is present but does not dominate across most of the genome.
Red flag
Batch explains a larger median fraction of variance than condition across most genes, or genes with high 'batch variance' overlap heavily with your top DE hits.
08Compare DE results with and without batch in the design formula
Fit the same DESeq2 model twice, once with only condition in the design and once with batch added, then compare the resulting gene lists and p-value distributions.
r
```
dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design = ~ batch + condition)
dds <- DESeq(dds)
res <- results(dds, name = "condition_trt_vs_untrt")
```
Healthy
Adding batch either leaves conclusions largely unchanged (if batch effects were small) or removes a chunk of previously 'significant' genes that were tracking batch, not biology, while the design remains full rank.
Red flag
resultsNames(dds) shows the batch coefficient could not be estimated (aliasing), or the gene list changes dramatically and unpredictably when batch is added, signaling near-confounding.

What to do about it

Add batch to the design formula (DESeq2)

When: Batch and condition are not confounded: your cross-tab from the checks shows multiple conditions represented within each batch.

Include batch as a factor term before condition in the design formula (design = ~ batch + condition), then run DESeq(dds) as usual and pull results for the condition contrast by name, e.g. results(dds, name="condition_trt_vs_untrt"). DESeq2 estimates a coefficient for batch and subtracts it during testing without ever modifying the raw counts.

Caveat: This consumes degrees of freedom, so you need enough samples per batch-by-condition cell to estimate both terms; with very few replicates per batch, the batch coefficient becomes unstable. Keep the variable of interest last in the formula to use DESeq2's default contrast behavior.

Additive batch term in edgeR

When: Same non-confounded scenario, but you're running edgeR instead of DESeq2, and you want quasi-likelihood testing.

Build a design matrix with model.matrix(~ group + batch), fit with glmQLFit(), and test with glmQLFTest() rather than glmLRT(), which gives fewer false positives when a batch term is present.

Caveat: An additive model assumes batch effects don't interact with treatment group; if you suspect batch changes the magnitude of the treatment effect itself, this assumption breaks and you'd need an interaction term, which requires more replicates than most bulk designs have.

removeBatchEffect for visualization only

When: You need a PCA plot or heatmap that shows biology cleanly for a figure or QC report, not for statistical testing.

Run limma::removeBatchEffect(expr, batch=batch, design=design) on the vst or rlog matrix, passing a design argument that specifies which biological terms to preserve so the function doesn't strip out real signal along with the batch effect.

Caveat: This output must never be fed into DESeq2, edgeR, or limma for hypothesis testing. Removing batch from the matrix collapses the degrees of freedom the model needs to correctly estimate uncertainty, producing artificially small p-values. Use it for plots only.

ComBat-seq for exploratory tools that need corrected integer counts

When: A downstream tool that isn't doing model-based DE testing (exploratory clustering, a co-expression network, an external classifier) requires an adjusted count matrix rather than a design-formula correction.

Run ComBat_seq(count_matrix, batch=batch_vector, group=condition_vector, full_mod=TRUE) from the sva package on the raw count matrix; it uses negative binomial regression to model both mean and dispersion batch effects and returns integer counts.

Caveat: Even though the output is integer-valued and technically loads into DESeq2 or edgeR, using it as DE testing input double-corrects: you've already removed batch from the counts and then the model estimates uncertainty as if it hadn't. Keep this path for exploratory or non-DE downstream uses, and use the design-formula approach for any results table you plan to report p-values from.

Redesign the experiment

When: Batch is completely confounded with condition: every sample in one batch is one condition and every sample in another batch is the other condition.

There is no statistical fix. Generate new samples that break the confound, for example by processing a subset of both conditions together in a new batch, or by re-running the comparison with a design that interleaves conditions across batches from the start.

Caveat: This costs time and possibly reagents, but any method applied to a confounded design (ComBat, removeBatchEffect, a batch term in the GLM) will produce a number, not a separable estimate of the biological effect. Reporting that number as if it were biology is worse than reporting the redesign delay.

When not to "fix" it

If your "batches" are actually your biological replicates, like different patient cohorts, different animals, or different tissue donors sampled and processed at different times because that's how the biology was collected, then the variance you're seeing between batches may be real inter-individual variability, not a technical artifact. Forcing that out with removeBatchEffect or ComBat erases genuine biological heterogeneity you may need to report or account for.

The same caution applies when what looks like a batch effect is actually a real subtype or cohort difference. When PCA on bulk RNA-seq mixes samples across expected group boundaries, for example some LUSC samples clustering with LUAD samples in TCGA data, that overlap can reflect genuine molecular subtype ambiguity rather than a processing artifact. Before you correct anything, check whether the "batch" variable actually correlates with a documented technical step, or whether it's standing in for a biological covariate (cohort, collection site, disease subtype) that you should be modeling explicitly instead of removing.

Five things experienced analysts do here

Record every batch variable at the bench when you prep the libraries, not after sequencing comes back; you cannot reconstruct kit lot or operator from the FASTQ files later.
Randomize which conditions get prepped and sequenced together. 'All controls this week, all treated next week' guarantees a confounded design no method can rescue.
Look at PCA colored by every batch variable you recorded before you run any differential expression test, not after you already have a gene list you like and are looking for reasons to trust it.
Never let removeBatchEffect() or ComBat output touch anything that computes a p-value; keep it strictly for plots, clustering, or non-statistical downstream tools.
Before trusting any batch coefficient, check that the design matrix is full rank with table(batch, condition) and that batch is stored as a factor, not a number.

Questions people ask

How do I check for batch effects in bulk RNA-seq?: Start by cross-tabulating your batch variables (prep date, flowcell, lane, kit lot) against condition to check the design isn't confounded, then run PCA and a sample-distance heatmap on the variance-stabilized matrix colored by both batch and condition. If PC1 tracks a batch variable instead of condition, or the table is block-diagonal, you have a problem to address before trusting any DE result.
Should I run ComBat before DESeq2?: No, not for differential expression testing. Pre-correcting the count matrix with ComBat or ComBat-seq and then feeding it into DESeq2 or edgeR causes incorrect degrees of freedom and artificially small p-values, because the model doesn't know the counts have already been adjusted. Put batch in the design formula instead (~ batch + condition) and let DESeq2 or edgeR subtract the effect during testing.
What's the difference between removeBatchEffect() and adding batch to the design formula?: removeBatchEffect() edits the expression matrix itself and is meant only for visualization: PCA plots, heatmaps, clustering figures. Adding batch to the design formula (~ batch + condition) leaves the raw counts untouched and instead has the negative binomial or linear model estimate a batch coefficient and subtract it during hypothesis testing, which is the statistically correct approach for anything that produces a p-value.
Can batch and condition be confounded in a way I can't fix?: Yes. If every sample in one batch belongs to only one condition and every sample in another batch belongs to the other condition, the design matrix is not full rank and no statistical method, including ComBat-seq or a design-formula batch term, can separate batch from biology. The only real fix is generating new samples that break the confound.
When should I not correct for a batch effect?: When what looks like a batch is actually a biological replicate structure (different patients, animals, or cohorts) or a real subtype difference rather than a technical artifact. Correcting it out in those cases removes genuine biological variability you may need to report, not noise.

Compare · DESeq2 vs edgeR: Which One Should You Use?

Sources

Batch Effect: To Correct or Not for Bulk RNA-seq Data — Distinction between removeBatchEffect for visualization vs. modeling batch in the design formula for statistical testing, and why pre-correction breaks degrees of freedom
Batch effects in DESeq2: determination and model building — Batch must be a factor, not numeric, in the DESeq2 design formula, and how the model subtracts batch without touching raw counts
Adjusting for confounding effects in edgeR — edgeR additive design matrix for batch, glmQLFit/glmQLFTest recommendation, and the confounded-design limit no method can fix
ComBat-seq: batch effect adjustment for RNA-seq count data — Why standard ComBat produces negative counts unsuitable for DESeq2/edgeR, and ComBat-seq's negative-binomial, integer-preserving alternative
What does evidence of batch effect in bulk rnaseq look like? — Recommendation to use variancePartition or BatchQC for quantitative batch assessment rather than visual inspection alone
Analyzing RNA-seq data with DESeq2 — Design formula syntax, variable-ordering convention (batch before condition), and lfcShrink usage

Part of the Batch effects series.