Batch Effect

A batch effect is systematic variation in a measurement that tracks the processing run rather than the sample: reagent lot, flow cell, extraction day, mass spec column age, technician, sensor unit. It is not noise in the statistical sense. It is structured, reproducible within a run, and large enough that in most unadjusted omics datasets the first principal component separates batches, not conditions.

How it works

Every assay has a chain of physical steps between the biology and the number in your file, and each step has parameters that drift. In RNA-seq: RNA integrity at extraction, poly-A selection efficiency, PCR cycle count, cluster density on the flow cell. In proteomics: digestion completeness, LC column degradation over an injection series, mass spectrometer calibration, and in isobaric labeling, which TMT plex a sample landed on. In mass spectrometry generally, instrument drift across an acquisition sequence is well characterized and produces intensity shifts that can exceed the biological differences under study ¹. Even imaging MS, where there is no library prep at all, shows batch structure from matrix application and instrument state ¹.

The statistical model most correction tools assume is additive and multiplicative on the log scale: for feature $g$ in sample $i$ in batch $b$, $y_{gib} = \alpha_g + X_i\beta_g + \gamma_{gb} + \delta_{gb}\epsilon_{gi}$. ComBat estimates $\gamma$ and $\delta$ with an empirical Bayes shrinkage across features, which is why it behaves better than per-gene batch means when you have few samples per batch. ComBat-seq does the analogous thing on raw counts with a negative binomial, so downstream tools like DESeq2 still see integers. Tools like BatchServer wrap this with evaluation plots and a per-feature score so you can see which features carry the batch signal before and after correction ².

The other family removes latent structure rather than a declared label: svaseq estimates surrogate variables from the expression matrix itself, and RUVSeq uses control genes (spike-ins, or genes you assert are unaffected by your variable of interest) to estimate unwanted factors. These are the right choice when batch is unrecorded or partially recorded.

In your own data

Where it shows up:

• RNA-seq. In the salmon/kallisto output or the featureCounts matrix. Compute log-CPM with edgeR::cpm(y, log=TRUE, prior.count=2), run prcomp on the top 1000 variable genes, and color points by sequencing run and by library prep date. If PC1 separates runs, you have it. Also check the per-sample fraction of reads in mitochondrial genes and the 5’/3’ coverage bias from picard CollectRnaSeqMetrics (MEDIAN_5PRIME_TO_3PRIME_BIAS). Degradation differences between draws masquerade as expression changes in long transcripts.
• Proteomics. Missingness is the tell. Plot the number of quantified proteins per run against acquisition order. A downward slope is column or source fouling, not biology. If your panel is Olink or SomaScan, check the bridge/control samples: they should be flat across plates, and if they are not, apply the vendor bridging factors before anything else.
• WGS. Batch shows up as depth and GC-bias differences, and in variant calls as batch-specific low-quality calls. If you ever re-sequence, do not merge two separately called VCFs. Re-call jointly from the CRAMs with GATK GenotypeGVCFs over all gVCFs, or you will read the caller’s prior as a new mutation.
• CGM. Every sensor insertion is a batch. Sensor-to-sensor offsets of 10–15 mg/dL are ordinary. Fit a per-sensor intercept before comparing weekly means, and discard the first 12 hours after insertion.

The mistake we see most: running ComBat with batch = draw_date on a personal time series. Your intervention, your season, your age, and your batch are the same vector. The correction removes the biology. The fix is design, not software. Randomize samples across runs, which for an n-of-1 means freezing aliquots and processing several timepoints together in one batch. Bank plasma at −80 °C and run six months of draws in one LC-MS sequence. When that is impossible, include a pooled reference aliquot (an identical split from one early draw) in every run and use its drift as the correction target.

Limitations

Correction cannot separate batch from a variable it is confounded with, and no method advertises when it has failed. Adjusted values also have understated variance, so feeding a ComBat-corrected matrix into a t-test inflates significance. Prefer putting batch in the design matrix (~ batch + condition) and let the model account for it. Reserve removeBatchEffect output for plots and clustering. Integrative multi-omics methods that jointly correct and embed ³ and single-cell approaches with learned constraints ⁴ can align distributions convincingly and still discard a real shared signal, because alignment is the objective they optimize. Aging studies in particular are vulnerable, since collection year and age move together across cohorts ⁵. If a corrected result would change a clinical decision, take the raw values and the batch structure to a physician rather than the adjusted matrix.

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.