How to Correct Batch Effects in Your Own RNA-Seq Data

By the end of this guide you will have three things: a gene-by-sample count matrix paired with a sample metadata table that records every technical variable you can recover, a diagnostic report that quantifies how much variance each of those variables explains, and either a statistical model that accounts for batch or a corrected expression matrix suitable for clustering and visualization. To follow along you need the raw counts rather than TPMs from someone else’s pipeline, R 4.3 or later with edgeR, limma, sva, and RUVSeq installed from Bioconductor, and whatever run-level metadata your sequencing provider will release: library preparation date, kit lot, flow cell and lane, RIN or DV200, input mass, and draw date. If you are profiling yourself across months or years, the draw date is usually the single most important covariate, and it is also the one most likely to be confounded with the biology you care about.

1. Write down the design before touching the counts

Before any code runs, it is worth being precise about what a batch effect is and why it can be unfixable. A batch effect is systematic, non-biological variation between groups of samples processed together: different library prep kit lots, different flow cells, different RNA extraction days, different technicians. The difficulty is not that the effect exists, because it always does, but that it can be correlated with the variable you want to measure. If every sample from January was prepped on one kit and every sample from July on another, no algorithm can separate season from kit. Correction methods redistribute variance according to a model you supply, and when that model is unidentifiable they will happily return a matrix that looks clean and means nothing.

For that reason, the first command you run should be a cross-tabulation rather than an analysis:

meta <- read.csv("metadata.csv", stringsAsFactors = FALSE)
table(meta$batch, meta$condition)
qr(model.matrix(~ batch + condition, data = meta))$rank
ncol(model.matrix(~ batch + condition, data = meta))

If the rank is less than the number of columns, the design is singular and batch and condition are aliased. In personal longitudinal profiling this is the default state of the world, because each draw tends to be its own batch. The fix is logistical rather than statistical: bank aliquots and prepare several timepoints in one library prep run, or carry a reference aliquot through every run so that batches share a common anchor. We return to that in step 8.

2. Assemble counts and metadata in a form you can re-run

The goal in this step is a reproducible starting point, which means raw counts at gene level plus a metadata table you can regenerate. If you have BAM files, featureCounts against a GENCODE GTF annotation is fast and predictable:

featureCounts -T 8 -p --countReadPairs -s 2 \
  -a gencode.v44.primary_assembly.annotation.gtf \
  -o counts.txt sample*.bam

The -s 2 flag sets reverse strandedness, which is correct for most Illumina stranded kits, TruSeq Stranded mRNA among them. Getting it wrong drops most of your reads, and the loss is often uneven across kits, so it manifests as a batch effect that is really a parameter error. If you quantified with salmon or kallisto instead, import with tximport(type = "salmon", tx2gene = t2g, countsFromAbundance = "lengthScaledTPM") so that the values behave like counts in the negative-binomial models used downstream.

Your metadata table should have one row per library, keyed on the same identifier used in the count matrix columns, with a column for every technical variable you can obtain. Record RIN, library concentration, total mapped reads, duplication rate, and intronic fraction alongside the batch labels. These continuous quality metrics frequently explain more variance than the nominal batch does, and they can be used as covariates directly.

3. Filter low-expression genes and normalize with TMM

Normalization and batch correction are separate problems that are easily conflated, so it helps to keep the distinction in view. Normalization removes differences in sequencing depth and in library composition, while batch correction removes gene-specific technical shifts. Do normalization first, always.

library(edgeR)
y <- DGEList(counts = cts, samples = meta)
keep <- filterByExpr(y, group = y$samples$condition)
y <- y[keep, , keep.lib.sizes = FALSE]
y <- calcNormFactors(y, method = "TMM")
logCPM <- cpm(y, log = TRUE, prior.count = 3)

TMM, the trimmed mean of M-values, starts by picking a reference sample. For every gene it computes the log-ratio M between the sample and the reference along with the mean log-expression A. It then trims the extremes of both distributions, where the edgeR defaults are logratioTrim = 0.3 and sumTrim = 0.05. The library’s scaling factor is a precision-weighted mean of the remaining M values. The assumption behind it is that most genes are not differentially expressed between any two libraries, so the bulk of the M distribution should center on zero. That assumption holds well for whole blood or a single tissue sampled across timepoints in one person.

It fails when a few transcripts dominate the library, which is exactly what happens with globin carryover in whole blood. HBB, HBA1, and HBA2 can take 50 to 70 percent of reads without globin depletion, and the fraction varies between draws. Check cpm(y)["ENSG00000244734", ] (HBB). If it swings widely, either move to a globin-depleted protocol going forward or exclude globin genes before calling calcNormFactors. The filtering step deserves more credit than it usually gets as well: genes near zero carry nearly all technical variance, and they inflate every batch diagnostic you are about to run, which is what filterByExpr is there to prevent.

4. Quantify the batch effect before deciding to correct it

Correction carries risk, so the decision to correct should follow a measurement rather than precede it. The minimum useful diagnostic is a principal component analysis of the log-CPM matrix restricted to the most variable genes, with the leading components regressed against each candidate covariate:

v <- apply(logCPM, 1, var)
top <- names(sort(v, decreasing = TRUE))[1:2000]
pc <- prcomp(t(logCPM[top, ]), scale. = FALSE)
pve <- pc$sdev^2 / sum(pc$sdev^2)

for (cv in c("batch", "condition", "RIN", "date", "lane")) {
  for (k in 1:5) {
    fit <- summary(lm(pc$x[, k] ~ meta[[cv]]))
    cat(sprintf("%s PC%d R2=%.2f p=%.2g\n", cv, k,
        fit$r.squared, pf(fit$fstatistic[1], fit$fstatistic[2],
                          fit$fstatistic[3], lower.tail = FALSE)))
  }
}

The shape of the result matters as much as its magnitude. A batch that explains 40 percent of PC1 and nothing of PC2 through PC5 is a different problem from one that smears across the first four components. For a more principled decomposition, variancePartition::fitExtractVarPartModel with a formula such as ~ (1|batch) + (1|subject) + RIN gives you the fraction of each gene’s variance attributable to each term, which is the number you want when deciding whether correction is worth the risk. Web tools such as BatchServer package the same evaluate-visualize-correct loop for people who would rather not write the regression themselves ¹. Work on public data has shown that even large, well-managed repositories like TCGA carry batch structure strong enough to change downstream conclusions if left unassessed ².

5. Model batch in the design rather than editing the counts

For differential expression, our recommendation is unambiguous: put batch in the design matrix and let the model estimate it alongside everything else. Doing so preserves the correct residual degrees of freedom, which a two-stage correct-then-test procedure does not, and it keeps the uncertainty in the batch estimate where it belongs.

design <- model.matrix(~ batch + condition, data = y$samples)
v <- voom(y, design, plot = TRUE)
fit <- lmFit(v, design)
fit <- eBayes(fit)
topTable(fit, coef = "conditionpost", number = 20)

When you have repeated measures of the same person, the equivalent structure is a block on subject:

design <- model.matrix(~ batch + condition, data = y$samples)
vw <- voom(y, design)
cor <- duplicateCorrelation(vw, design, block = y$samples$subject)
vw <- voom(y, design, block = y$samples$subject, correlation = cor$consensus)
cor <- duplicateCorrelation(vw, design, block = y$samples$subject)
fit <- eBayes(lmFit(vw, design, block = y$samples$subject,
                    correlation = cor$consensus))

The two-pass duplicateCorrelation is deliberate. The first pass estimates the within-subject correlation from a naive voom fit, and the second refits the precision weights under that correlation. A consensus correlation above roughly 0.4 in a personal time series is normal and means most genes are stable within you and variable between people.

6. Produce a corrected matrix only when you need one

There are legitimate reasons to want a corrected matrix: clustering, heatmaps, UMAP, correlation networks, or handing data to a tool that cannot accept a covariate. For log-scale work, limma::removeBatchEffect is the right tool, and the design argument is not optional:

design0 <- model.matrix(~ condition, data = y$samples)
logCPM_adj <- removeBatchEffect(logCPM, batch = y$samples$batch,
                                covariates = y$samples$RIN,
                                design = design0)

Passing design tells the function which variation to protect. Omit it and the batch estimate absorbs any biological signal that happens to align with batch, and you will have quietly deleted your result.

When a downstream tool requires integers, use ComBat-seq. It fits a negative binomial regression per gene with batch as a covariate and maps the observed counts back onto a batch-free distribution, returning integers rather than continuous adjusted values:

library(sva)
adj_counts <- ComBat_seq(as.matrix(y$counts),
                         batch = y$samples$batch,
                         group = y$samples$condition)

One caveat is worth understanding before you rely on the output. Empirical Bayes shrinkage of the dispersion is what makes ComBat-family methods work in small batches. It can misbehave when one batch has only a handful of samples, or when a gene is near-zero in one batch and expressed in another. ComBat-ref addresses part of this by pooling dispersion toward a selected reference batch rather than shrinking all batches toward a common center, which improved sensitivity both in simulation and in real datasets ³. More recent work such as AMDBNorm continues along the same line, adjusting the negative binomial parameterization used for count-level correction ⁴. We would use ComBat-seq or ComBat-ref when producing counts is a hard requirement, and the design-matrix approach otherwise.

Whichever route you take, run topTable on both the modeled and the corrected path and compare the results. If they disagree substantially, the correction is doing something you have not yet understood.

7. Handle batches you cannot label

Sometimes the batch variable is unknown or only partly recorded, because you pooled data from two providers or an instrument drifted mid-run. Surrogate variable analysis handles this case by estimating the hidden structure from the data itself:

mod  <- model.matrix(~ condition, data = y$samples)
mod0 <- model.matrix(~ 1, data = y$samples)
dat  <- cpm(y, log = FALSE)
svs  <- svaseq(dat, mod, mod0)   # or n.sv = 2 to fix the number
design <- cbind(mod, svs$sv)

The alternative is RUV, which uses genes assumed to be unaffected by your condition as an empirical control set:

library(RUVSeq)
set <- newSeqExpressionSet(round(dat))
controls <- rownames(dat)[!rownames(dat) %in% topTable(fit, n = 5000)$ID]
set <- RUVg(set, controls, k = 2)

Both approaches are sensitive to the choice of k. Increasing the number of factors monotonically reduces apparent batch structure and monotonically increases the chance that you have removed real signal. Choose k by the point at which batch-associated variance in the PCA stops dropping sharply, not by the point at which your gene of interest becomes significant. A systematic comparison of correction methods in human transcriptome data found that performance depends heavily on the structure of the study and that quality-control samples processed alongside the biological samples materially improve the outcome ⁵.

8. Anchor every batch with a shared reference sample

This is the single highest-return intervention in longitudinal personal profiling, and it is a wet-lab decision rather than a computational one. Extract a large quantity of RNA once or bank a pool of PBMCs. Split it into single-use aliquots and include one aliquot in every library preparation run. Because the reference is biologically identical across runs, any difference you observe in it is technical by construction, and it gives you a direct estimate of the batch shift that is independent of your condition. The same logic motivates including QC samples in epidemiological transcriptome studies ⁵, and the same principle underlies methods that use a shared modality or shared observations to bridge otherwise unlinkable batches across omics types ⁶.

In practice, that means drawing into PAXgene or Tempus tubes and keeping the freezer chain consistent. It also means avoiding re-freezing and randomizing the order in which timepoints are prepped rather than processing them chronologically. Randomization costs nothing and converts an unfixable confound into a nuisance variable.

9. Check that the correction preserved biology

Three checks are worth running every time, because a correction that looks successful in a PCA can still have damaged your data. First, confirm that variance attributed to batch has fallen in the PCA and that variance attributed to your variable of interest has not. Second, verify that known biology survives. Sex-linked genes such as XIST, RPS4Y1, and KDM5D should still separate cleanly if your data include both sexes. In a personal time series, genes with well-characterized circadian or postprandial behavior should still track draw time. Third, run a null test by permuting the condition labels within batch and refitting. A correctly specified pipeline yields a roughly uniform p-value histogram under permutation, whereas a pipeline that has laundered batch structure into apparent signal yields a spike near zero.

Interpretation stops here. An expression change in a pathway is a measurement, and mapping it to anything about your health requires a clinician with your full clinical picture. Nothing in a corrected count matrix is a diagnosis.

Common problems

The failure modes below account for most of the trouble people run into, and each has a specific response.

• Batch is perfectly confounded with condition. There is no statistical remedy. Report the effect as uninterpretable and redesign the collection with shared aliquots or randomized prep order.
• One batch has two samples. Empirical Bayes shrinkage in ComBat-family methods becomes unstable with very small batches, producing corrected values dominated by the prior. Merge small batches that share a kit lot and date, or drop batch as a factor and use a continuous surrogate such as RIN or duplication rate.
• The correction improved the PCA but the differential expression results got worse. This usually means design was omitted in removeBatchEffect, or that k in RUV or SVA was set too high. Refit with the protected design and step k down.
• Corrected values are negative. removeBatchEffect works on the log scale and returns continuous, sometimes negative, values. They are valid for clustering and correlation and invalid as input to any tool expecting counts. Use ComBat-seq for that path.
• Cell composition is masquerading as batch. In whole blood, neutrophil fraction routinely drives PC1 and varies with draw time, acute infection, and exercise. Estimate proportions by deconvolution and include them as covariates, since removing them as batch discards real biology.
• You are extending the method to single-cell or multi-omic data. The estimators differ enough that bulk tools transfer poorly. Single-cell correction methods now typically operate on learned embeddings with explicit structure-preservation objectives ⁷, and multi-omic integration adds the problem of aligning batches that share no common feature space ⁸. Use a single-cell-specific quality control workflow before correcting anything ⁹.

Woolf Software builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

1. Tiansheng Zhu, Rui Sun, Fangfei Zhang, et al. BatchServer: A Web Server for Batch Effect Evaluation, Visualization, and Correction. Journal of Proteome Research, 2020. https://doi.org/10.1021/acs.jproteome.0c00488 ↩

2. Sadra Salehi-Mazandarani. The Necessity of Batch Effect Assessment and Correction in TCGA RNA-Seq Data. 2026. https://doi.org/10.20944/preprints202609.0747.v1 ↩

3. Xiaoyu Zhang. Highly effective batch effect correction method for RNA-seq count data. Computational and Structural Biotechnology Journal, 2024. https://doi.org/10.1016/j.csbj.2024.12.010 ↩

4. Xu Zhang, Wenshi Qiu, Feng Qiao, et al. AMDBNorm powers batch effects correction for RNA-seq count data. International Journal of Biomathematics, 2026. https://doi.org/10.1142/s1793524526500816 ↩

5. Almudena Espín-Pérez, Chris Portier, Marc Chadeau-Hyam, et al. Comparison of statistical methods and the use of quality control samples for batch effect correction in human transcriptome data. PLOS ONE, 2018. https://doi.org/10.1371/journal.pone.0202947 ↩ ↩²

6. Manuel Ugidos, Sonia Tarazona, José M Prats-Montalbán, et al. MultiBaC: A strategy to remove batch effects between different omic data types. Statistical Methods in Medical Research, 2020. https://doi.org/10.1177/0962280220907365 ↩

7. Yongjie Xu, Zelin Zang, Jun Xia, et al. Structure-preserving visualization for single-cell RNA-Seq profiles using deep manifold transformation with batch-correction. Communications Biology, 2023. https://doi.org/10.1038/s42003-023-04662-z ↩

8. Tiezheng Qiao, Fei Teng, Huiyong Zhang, et al. scAQUA: A Quintuplet Constraint-Based Batch Effect Correction Method for Single-Cell Multi-omics Data. Lecture Notes in Computer Science, 2026. https://doi.org/10.1007/978-981-92-3716-6_30 ↩

9. Daihan Ji, Mei Han, Shuting Lu, et al. SingleCellMQC: A comprehensive quality control workflow for single-cell multi-omics. iScience, 2026. https://doi.org/10.1016/j.isci.2026.117398 ↩