How do I check for batch effects in spatial transcriptomics data?

Start by confirming batch metadata (slide, capture area, processing date) was actually recorded, then color a PCA or UMAP plot by sample_id before any correction. If PC1 or PC2 separates by slide instead of cell type or condition, you have a batch effect. Back that up with a VlnPlot of nCount_Spatial per sample and a sample-to-sample correlation heatmap.

Should I correct batch effects before or after clustering spots?

Run your batch diagnostics (PCA colored by sample, depth violin plots) before any clustering, then decide. If correction is warranted, apply it to the embedding (Harmony) or the neighbor graph (bbknn) before clustering so the clusters themselves reflect corrected structure, not raw batch-driven groupings.

Does ComBat-seq work for spatial transcriptomics?

It works on count matrices in general, using a negative-binomial model to adjust batch effects while preserving integer counts for edgeR or DESeq2. It was developed for bulk-style RNA-seq counts, so applying it to sparse spot-level or single-cell-level spatial counts needs care, and the field hasn't settled on the best stage (raw counts, normalized, or pseudobulk) to apply it at for spatial data specifically.

What's the difference between Harmony and ComBat-seq for spatial data?

Harmony corrects the PCA embedding only, iteratively adjusting cell/spot coordinates across batches; it never touches the expression matrix, so it's suited for clustering and visualization. ComBat-seq adjusts the actual count matrix using a negative-binomial model, producing corrected integer counts you can feed directly into a differential expression tool.

When should I not correct for batch effect in spatial transcriptomics?

Don't correct when batch is perfectly confounded with your biological condition, for example all controls on one slide and all treated on another; no method can separate the two, and the honest fix is redesigning the experiment. Also don't correct when the pattern reflects real regional biology, like a necrotic tumor core, rather than a technical artifact.

Sanity check · Spatial Transcriptomics

How to Detect Batch Effects in Spatial Transcriptomics

A slide is not a replicate, learn to tell tissue-section artifact from real biology before you trust a single cluster.

By Ming "Tommy" Tang, Director of Bioinformatics in Big Pharma · Reviewed September 2026 · 5 min read

You ran Visium or Xenium on several tissue sections, some control and some treated, and now you're staring at a PCA or UMAP plot expecting to see cell types separate out. Instead PC1 splits cleanly by slide. Or a cluster that looked like an exciting new cell state turns out to contain every spot from one capture area and nothing else. This is the single most common way a spatial transcriptomics analysis goes wrong before it even gets to biology.

The stakes are not abstract. Each tissue section is its own batch and its own biology at the same time, which means the confound is baked in from the moment you loaded the slide. If processing date, permeabilization time, or capture area lines up with your experimental groups, no downstream tool can cleanly separate "this gene changed because of treatment" from "this gene changed because this spot came from a lower-quality section." The genotyping longevity study that got retracted from Science is the cautionary tale: platform artifacts got read as biology, made the news, and didn't replicate.

This page gives you an ordered set of checks you can run in under an hour on a multi-sample Visium or Xenium object, from checking whether batch was even recorded, to looking at PCA colored by sample, to running Harmony and checking whether it actually helped. It ends with the harder call: when correction is the right move, and when the honest answer is that your comparison is confounded and no amount of RunHarmony fixes that.

What it looks like when it's happening

PC1 or PC2 in a PCA/UMAP plot separates by sample_id or slide, not by cell type or condition
SpatialDimPlot colored by sample_id shows each tissue section as its own tight island instead of mixed neighborhoods of shared cell types
VlnPlot of nCount_Spatial or nFeature_Spatial shows one slide's median UMI count several-fold higher or lower than the others
A cluster that looks biologically distinct turns out to be composed almost entirely of spots from a single capture area
Top PC1 loadings are dominated by mitochondrial (MT-) or ribosomal (RPS/RPL) genes rather than known marker genes
Differential expression between conditions returns hundreds of genes, and most of them correlate with per-spot depth rather than the biology you're testing
A sample-to-sample correlation heatmap clusters tightly by processing date or slide instead of by experimental condition

Why it happens

Every tissue section on a Visium slide or Xenium run goes through its own dissociation, permeabilization, and capture chemistry, and those steps are never perfectly identical between sections even when the same protocol is followed on the same day. Permeabilization time changes RNA capture efficiency spot by spot; a slightly under- or over-permeabilized section yields systematically lower or higher UMI counts across the whole capture area, not in a pattern that tracks biology. That single technical variable, sequencing depth per spot, is one of the strongest drivers of variance in single-cell and spot-based data, and when it differs across sections it looks exactly like an "on/off" biological pattern: cells or spots with more UMIs simply detect more genes and cluster together in PCA/UMAP purely by library size.

The mechanism compounds downstream. Once depth or capture efficiency differs systematically between sections, genes that have nothing to do with each other biologically start appearing co-regulated, because they're all responding to the same technical driver, RNA quality or capture efficiency, at the same time. This is the same failure mode documented for bulk RNA-seq: batch can explain more variance than the biological factor of interest, and it can even flip the direction of an apparent group difference. In spatial data it's worse in one specific way, because "batch" and "sample" are usually the same axis, there's no way to average over replicates within a batch the way you can pool multiple bulk libraries from one prep.

The genotyping-platform story is the same problem in a different assay: combining data from multiple platforms without accounting for their technical differences let platform-specific artifacts get read as genetic signal, and it took independent replication with consistent platforms to make the associations vanish. Spatial transcriptomics inherits this risk directly, because slide, capture area, kit lot, and operator are all platform-adjacent variables that vary section to section, and unless they're recorded in metadata at collection time, there is no way to test for their effect after the fact.

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

01Confirm batch metadata was actually recorded
Check that your object's metadata includes the variables that can create a batch effect: slide, capture area, permeabilization time, sequencing lane, kit lot, processing date, and operator. If these aren't columns in colData or obs, you cannot diagnose a batch effect after the fact, only guess at one.
r
```
colnames(colData(spe))
# expect to see: sample_id, slide, capture_area,
# permeabilization_time, seq_lane, kit_lot, processing_date, operator
```
Healthy
All relevant technical variables are present as explicit metadata columns, one row per spot/cell.
Red flag
Only sample_id exists, or slide/date/operator were never captured, meaning you're diagnosing blind.
02Cross-tabulate batch against biological condition
Build a table of condition versus batch variable (sample_id, slide, or processing date). This is the single most important check, because it tells you whether correction is even possible before you run anything.
r
```
table(colData(spe)$condition, colData(spe)$sample_id)
```
Healthy
Every condition appears across multiple batches, so batch and biology are not perfectly aligned.
Red flag
Each condition maps to exactly one batch or one processing date (all controls on slide 1, all treated on slide 2), a perfect confound that no algorithm can undo.
03Check per-spot depth distributions across samples
Plot nCount_Spatial and nFeature_Spatial grouped by sample_id. This is the cheapest diagnostic for the depth-driven batch effects that dominate spatial data.
r
```
VlnPlot(brain, features = "nCount_Spatial", group.by = "sample_id")
VlnPlot(brain, features = "nFeature_Spatial", group.by = "sample_id")
```
Healthy
Distributions overlap substantially across samples, with variation you'd expect from normal tissue heterogeneity.
Red flag
One or more slides show a median UMI count several-fold higher or lower than the rest, a depth batch effect, not biology.
04Overlay depth on the tissue image
Plot nCount_Spatial spatially to see whether the depth gradient follows tissue anatomy or follows the geometry of the capture area.
r
```
SpatialFeaturePlot(brain, features = "nCount_Spatial")
```
Healthy
A gradient that tracks known tissue structure, for example higher counts in gray matter than white matter.
Red flag
A uniform gradient across the entire capture area, or a sharp boundary that matches the edge of the slide rather than any tissue feature, pointing to permeabilization or capture artifact.
05Color PCA/UMAP by batch before any correction
Run PCA or UMAP on the merged, normalized object and color by sample_id/batch, then separately by cell type or condition. Do this before touching Harmony or any integration method, you need a documented "before" to know if correction did anything.
r
```
DimPlot(brain_merged, reduction = "pca", group.by = "sample_id")
DimPlot(brain_merged, reduction = "pca", group.by = "condition")
```
Healthy
Spots from different samples overlap and mix within shared cell-type clusters; condition, not sample, drives any separation you care about.
Red flag
PC1 or PC2 cleanly separates spots by sample_id/slide rather than by cell type or condition.
06Build a sample-to-sample correlation heatmap
Pseudobulk each sample (sum or average counts per gene), compute pairwise correlation, and hierarchically cluster with rows/columns labeled by both batch and condition.
Healthy
Samples cluster primarily by biological condition.
Red flag
Samples cluster tightly by processing date or slide regardless of condition, batch is the dominant axis in the correlation structure.
07Inspect what's driving PC1
Pull the gene loadings for the top PCs and check whether the highest-magnitude genes are known cell-type markers or technical-quality genes.
r
```
sort(Loadings(brain[["pca"]])[, 1]) |> tail(20)
```
Healthy
Top loadings are biologically interpretable marker genes for the tissue you're profiling.
Red flag
Top loadings are dominated by mitochondrial genes (MT-), ribosomal genes (RPS/RPL), or other RNA-quality/stress markers, a sign that PC1 is tracking tissue quality, not cell identity.
08Run Harmony and check before/after agreement with ground truth
Apply Harmony to the PCA embedding specifying your batch variables, then compare clustering before and after against any manual histology annotation you have, using Adjusted Rand Index (ARI) as the agreement metric.
r
```
harmony_embed <- harmony::RunHarmony(
  data_mat = reducedDim(spe, "PCA_M1_lam0.2"),
  meta_data = colData(spe),
  vars_use = c("sample_id", "subject_id"),
  max_iter = 20,
  verbose = FALSE
)
SpatialDimPlot(brain, group.by = "sample_id")  # before vs after
```
Healthy
After Harmony, samples mix within shared clusters and ARI against manual annotation goes up relative to the uncorrected embedding.
Red flag
Samples still separate along PCs after Harmony, or ARI doesn't improve, a sign batch is confounded with biology rather than merely offset from it.

What to do about it

Harmony on the PCA embedding

When: Batch is not confounded with condition (your cross-tab check showed every condition spread across multiple batches), and you mainly need clean clustering and visualization.

Merge samples, normalize with SCTransform, run PCA, then run harmony::RunHarmony() with vars_use set to sample_id (and subject_id if applicable), and re-cluster on the corrected embedding.

Caveat: Harmony adjusts the embedding only, never the expression values themselves, it will clean up your UMAP/clustering but you still can't run differential expression directly on Harmony-corrected coordinates.

ComBat-seq on the count matrix

When: You need batch-adjusted counts for differential expression testing, not just a cleaner embedding, typically after pseudobulking spots or cells by sample.

Run ComBat_seq(count_matrix, batch = batch_vector, group = design_vector); it models batch effects with a negative-binomial approach and quantile mapping, and returns integer counts you can feed straight into edgeR or DESeq2.

Caveat: It was built for bulk-style count matrices; applying it to highly sparse spot- or cell-level counts needs care, and there's no settled guidance yet on whether to correct at the raw count, normalized, or pseudobulk stage for spatial data specifically.

Batch-balanced kNN (bbknn) for multi-section integration

When: You're integrating many tissue sections in Scanpy and want a batch-aware neighbor graph for joint clustering rather than an embedding correction.

sc.pp.bbknn(adata, batch_key='sample_id', neighbors_within_batch=15) finds top neighbors within each section separately, then merges results into one graph before clustering.

Caveat: It changes the neighbor graph, not the expression values, and cluster granularity is sensitive to neighbors_within_batch, retune it if clusters look artificially fragmented or merged.

SCTransform normalization before any batch correction

When: Symptoms point to depth-driven separation (VlnPlot shows sharply different UMI/gene counts per slide) rather than a full-blown slide-level confound.

Normalize each section with SCTransform (Seurat v5's default for spatial data), which fits a regularized negative-binomial model to remove technical variance in depth before you touch PCA or integration.

Caveat: Normalization reduces depth-driven separation but doesn't remove slide-specific chemistry differences like permeabilization time, you'll likely still need Harmony or bbknn on top of it.

Spatially-aware QC to catch local, not global, artifacts

When: The artifact is confined to a region of one section (a torn corner, an over-permeabilized patch) rather than affecting an entire slide uniformly.

Use a spatially-aware QC method (e.g. SpotSweeper) that computes neighborhood-aware z-scores for total counts, detected genes, and mitochondrial fraction, comparing each spot to its local neighbors instead of a tissue-wide threshold.

Caveat: This is newer, less broadly validated tooling; tune the spatial neighborhood size deliberately, since too large a neighborhood will wash out the local artifact you're trying to catch.

Randomize the design for your next experiment

When: You're planning the next round of sample collection, not patching an already-confounded dataset.

Spread conditions across slides, sequencing lanes, and operators instead of processing all controls on one slide and all treated on another; if full randomization isn't feasible, at minimum balance groups across batches, and record every technical variable (slide, capture area, permeabilization time, lane, kit lot, date, operator) in metadata at collection time.

Caveat: This only prevents future confounding, it does nothing for data you've already collected with condition and batch perfectly aligned.

When not to "fix" it

Don't correct when your cross-tab shows batch and condition perfectly aligned, for example every control section processed on one date and every treated section on another. At that point "correcting for batch" and "correcting for treatment" are the same operation, and any method that removes the batch axis removes your biology with it. The honest move is to flag the confound, report it, and redesign the next round of sample collection with randomization, not to apply Harmony or ComBat-seq and report the result as if it were clean.

Also don't correct when the pattern you're seeing is real regional biology and not a technical artifact, a necrotic tumor core versus a viable edge will show a genuine gradient in RNA quality and depth that happens to look like the depth-driven batch signatures described above, but it reflects actual tissue state, not permeabilization variance. Check whether the "artifact" boundary follows the capture-area geometry (technical) or known tissue anatomy (biological) before deciding it needs correcting.

Five things experienced analysts do here

Record slide, capture area, permeabilization time, sequencing lane, kit lot, processing date, and operator in metadata the moment you collect the data, you cannot diagnose a batch effect after the fact if the variable was never captured.
Cross-tabulate batch against condition before you touch any correction tool; a diagonal table means you have a redesign problem, not an analysis problem, and no algorithm changes that.
Always generate the "before" PCA/UMAP colored by batch before running Harmony or bbknn, and look at the "after" plot side by side, a correction you can't show working before-and-after is a claim, not a result.
Keep embedding correction (Harmony, bbknn) and count correction (ComBat-seq) mentally separate: use the former for clustering and visualization, the latter when you actually need adjusted counts for a DE test.
When you have ground truth, like manual histology annotation, validate correction against it using something like ARI; a correction that doesn't improve agreement with known anatomy may just be smoothing away real biology.

Questions people ask

How do I check for batch effects in spatial transcriptomics data?: Start by confirming batch metadata (slide, capture area, processing date) was actually recorded, then color a PCA or UMAP plot by sample_id before any correction. If PC1 or PC2 separates by slide instead of cell type or condition, you have a batch effect. Back that up with a VlnPlot of nCount_Spatial per sample and a sample-to-sample correlation heatmap.
Should I correct batch effects before or after clustering spots?: Run your batch diagnostics (PCA colored by sample, depth violin plots) before any clustering, then decide. If correction is warranted, apply it to the embedding (Harmony) or the neighbor graph (bbknn) before clustering so the clusters themselves reflect corrected structure, not raw batch-driven groupings.
Does ComBat-seq work for spatial transcriptomics?: It works on count matrices in general, using a negative-binomial model to adjust batch effects while preserving integer counts for edgeR or DESeq2. It was developed for bulk-style RNA-seq counts, so applying it to sparse spot-level or single-cell-level spatial counts needs care, and the field hasn't settled on the best stage (raw counts, normalized, or pseudobulk) to apply it at for spatial data specifically.
What's the difference between Harmony and ComBat-seq for spatial data?: Harmony corrects the PCA embedding only, iteratively adjusting cell/spot coordinates across batches; it never touches the expression matrix, so it's suited for clustering and visualization. ComBat-seq adjusts the actual count matrix using a negative-binomial model, producing corrected integer counts you can feed directly into a differential expression tool.
When should I not correct for batch effect in spatial transcriptomics?: Don't correct when batch is perfectly confounded with your biological condition, for example all controls on one slide and all treated on another; no method can separate the two, and the honest fix is redesigning the experiment. Also don't correct when the pattern reflects real regional biology, like a necrotic tumor core, rather than a technical artifact.

Guide · How to Detect Integration Over-Correction in Spatial Transcriptomics
Guide · How to Choose Cell QC Thresholds in Spatial Transcriptomics
Guide · How to Sanity-Check Marker Genes and Cell Type Labels in Spatial Transcriptomics
Guide · How to Tell If You Overclustered in Spatial Transcriptomics
Guide · How to Avoid Pseudoreplication in Spatial Transcriptomics
Glossary · Pseudobulk
Glossary · Spatial transcriptomics
Glossary · Visium

Sources

Analysis, visualization, and integration of spatial datasets with Seurat — SCTransform normalization and VlnPlot/SpatialFeaturePlot depth diagnostics
Spatial data integration with Harmony (10x Visium Human DLPFC) — RunHarmony code and before/after ARI validation approach
ComBat-seq: batch effect adjustment for RNA-seq count data — ComBat-seq negative-binomial modeling for count-level batch correction
Scanpy bbknn documentation — Batch-balanced kNN parameters for multi-section integration
Spotsweeper-py: spatially-aware quality control metrics for spatial omics data in the Python ecosystem — Neighborhood-aware QC for local batch artifacts
Towards a Better Understanding of Batch Effects in Spatial Transcriptomics: Definition and Method Evaluation — Alignment vs. integration framing for spatial batch correction

Part of the Batch effects series.