How do I check for batch effects in scATAC-seq?

Color your LSI or UMAP embedding by known technical variables, lane, prep date, operator, not just by condition, and look for a dominant axis that tracks those instead of cell type. Back that up with DepthCor() to see which LSI components correlate with per-cell fragment counts, and cross-tabulate cluster membership by batch to catch skewed proportions numerically, not just visually.

Should I use Harmony or ComBat for scATAC-seq batch correction?

Harmony is the standard choice because it operates directly on the LSI embedding that Signac or ArchR already produce, iteratively adjusting coordinates per batch without needing a peak-level statistical model. ComBat-style methods built for continuous, roughly normal bulk data don't map cleanly onto sparse, near-binary peak counts. Note Harmony corrects only the embedding, not the peak matrix; for cases where linear methods like Harmony or MNN fail outright, epiConv is a documented alternative that skips the embedding step entirely.

When should I not correct for batch effects in scATAC-seq?

When batch is confounded with your biological condition, every control sample in one batch, every treatment sample in another, no algorithm, Harmony included, can separate the two. Forcing a correction hides the confound rather than resolving it. The honest fix is redesigning the experiment or explicitly reporting the confound as a limitation, not running a batch-correction tool and trusting the output.

What does a batch effect look like on a scATAC-seq PCA or LSI plot?

A dominant component, frequently LSI2, after you've already dropped the depth-driven LSI1, separates cells by sequencing date, lane, or operator rather than by cell type. DepthCor often shows that same axis still correlating with fragment counts, meaning what looks like 'batch' is frequently sequencing depth wearing a batch label.

Sanity check · Single-Cell ATAC-seq

How to Detect Batch Effects in Single-Cell ATAC-seq

Q: Is a batch effect in scATAC-seq the same problem as in bulk or single-cell RNA-seq?

The underlying idea is the same, systematic technical variation masquerading as biology, but the mechanism differs. scATAC-seq batch effects are driven heavily by nuclei:Tn5 stoichiometry at tagmentation, and the peak-by-cell matrix is sparser and closer to binary than a gene-by-cell matrix, so depth swamps LSI component 1 more severely than it swamps PC1 in RNA-seq. Peak sets themselves can also disagree between batches, a failure mode fixed-annotation RNA-seq doesn't have.

Your LSI plot separates by processing date instead of cell type, here's the order of checks that tells you whether that's a fixable artifact or a confound no tool can undo.

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

You extended your scRNA-seq pipeline to the 10x ATAC or Multiome kit, ran iterative LSI the way the Signac or ArchR vignette walks you through, and colored the UMAP by sample. Instead of the expected islands of T cells, B cells, and monocytes, you get islands labeled "batch 1" and "batch 2." Your first instinct, reach for Harmony, is reasonable, but premature.

If you correct before checking whether batch lines up with your biological groups, you can either erase a real biological difference or paper over a confound that no algorithm can fix. Batch effects in scATAC-seq are also a harder problem than the RNA-seq version you may already know: the peak-by-cell matrix is sparser and closer to binary, so depth differences between batches dominate the first LSI component almost every run, and peak sets called separately per batch can disagree badly enough that the same cell type looks like two different things.

This page gives you an ordered set of checks, cheap metadata audits first, embedding diagnostics second, peak-set comparisons last, so that within the hour you know whether you're looking at a correctable technical artifact, a confound you need to redesign around, or biology that a batch-correction tool would wrongly smooth away.

What it looks like when it's happening

UMAP or LSI plot where cells cluster by sample, sequencing date, or operator instead of by expected cell type
DepthCor() shows a component beyond the usual LSI1 (often LSI2) still strongly correlated with per-cell fragment counts
Fragments-per-cell distributions differ several-fold between batches (e.g. median 3,000 vs 15,000 fragments/nucleus)
A cell type abundant in one batch is nearly absent from another when you cross-tabulate cluster membership by batch
Peak sets called separately per batch overlap poorly (low Jaccard) for what should be the same cell type
After running Harmony, batches mix in the embedding but known marker gene activity scores also flatten out across clusters
QC pass rates (cells retained after fragment/TSS enrichment thresholds) differ noticeably by processing batch

Why it happens

Tn5 is a single-turnover enzyme: once it cuts, it's done. The number of fragments a nucleus yields is set largely by the nuclei:Tn5 stoichiometry at tagmentation. Any batch-to-batch drift in that ratio, a different operator pipetting slightly differently, a new reagent lot, a different nuclei prep day, produces a systematic difference in fragments per cell between batches, not a random one. That's the raw material every downstream artifact is built from.

LSI is TF-IDF-normalized peak counts run through SVD, and its first component tracks total signal per cell, i.e. depth, almost by construction. This is true in scRNA-seq too, but it's more severe in scATAC-seq: the peak-by-cell matrix is far sparser and closer to binary than a gene-by-cell matrix, so depth heterogeneity has more room to dominate the decomposition. Practitioners drop LSI component 1 for exactly this reason. The trap is assuming that's the only place depth hides, when depth differences are large, correlation with depth can bleed into component 2 as well, and that's the axis people mistake for "batch."

There's a second mechanism specific to ATAC: the feature space itself, the peak set, depends on the data you call it from. Peaks called on pooled data under-call anything specific to a rare cell type present in only one batch. Peaks called separately per batch, meant to fix that, instead produce peak sets that disagree with each other for the same abundant cell types, because differing depth and signal-to-noise between batches change which regions clear the calling threshold. So before you've touched normalization or integration, the columns of your matrix already encode batch.

Harmony, the standard fix, only adjusts the LSI/PCA embedding coordinates, it does not recompute the peak matrix or produce a corrected peak set you can project new cells onto. If the underlying peak-by-cell matrix is already batch-confounded from the calling step, "corrected" clusters still rest on a shaky feature space, and running differential accessibility on the raw counts within those clusters can reintroduce the exact artifact the embedding correction appeared to remove.

The checks

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

0/7 checked · saved in this browser

01Audit whether batch is even annotated, and whether it's confounded with condition
Before plotting anything, confirm your per-cell or per-sample metadata records sequencing lane, library prep date, operator, and reagent/kit lot, not just sample ID and condition. Then cross-tabulate each batch variable against your biological condition.
Healthy
Every batch variable has more than one biological condition represented in it, and vice versa, the design is not confounded.
Red flag
Any batch variable perfectly predicts condition (all controls in batch 1, all treated in batch 2), or you discover lane/date was never recorded at all, you cannot adjust for a batch you never annotated.
02Compare fragments-per-cell (depth) distributions across batches
Compute median fragments per cell per batch and plot the distributions side by side.
r
```
atac@meta.data %>%
  group_by(batch) %>%
  summarise(median_frags = median(nCount_peaks), n_cells = n())

VlnPlot(atac, features = "nCount_peaks", group.by = "batch", pt.size = 0)
```
Healthy
Median fragments/cell are within roughly 2-fold of each other across batches, with substantial overlap in the violins.
Red flag
One batch's median is several-fold higher or lower than another's (e.g. 3,000 vs 15,000 fragments/nucleus), a sign the nuclei:Tn5 ratio differed at tagmentation.
03Confirm exactly which LSI component tracks depth
Run DepthCor to see how each LSI component correlates with total fragment count per cell, and inspect the correlation plot component by component rather than assuming only component 1 is affected.
r
```
depth_cor <- DepthCor(atac, k = 30)
depth_cor
```
Healthy
Component 1 shows the well-known strong correlation with depth (expected, normal, and why it's dropped), and no later component in your working range exceeds roughly |r| > 0.5 with depth.
Red flag
Component 2 or 3 is also strongly depth-correlated, depth is leaking further into the embedding than the standard 'drop component 1, use 2:30' fix accounts for.
04Color LSI/UMAP embeddings by batch and by cell type side by side
Plot the LSI embedding (excluding component 1) and the UMAP colored by batch, then again colored by cell type, and compare which axis separates what.
r
```
DimPlot(atac, reduction = "lsi", dims = c(2, 3), group.by = "batch")
DimPlot(atac, reduction = "umap", group.by = "batch")
DimPlot(atac, reduction = "umap", group.by = "celltype")
```
Healthy
Batch labels are mixed within each cluster; the axis separating cell types is not the same axis separating batch.
Red flag
A dominant axis, often LSI2 or UMAP1, cleanly separates sequencing date, lane, or operator rather than cell type.
05Cross-tabulate cluster composition by batch
Build a contingency table of cluster assignment against batch, convert to proportions, and test for independence.
r
```
tab <- table(atac$seurat_clusters, atac$batch)
round(prop.table(tab, margin = 2), 2)
chisq.test(tab)
```
Healthy
Each cluster draws cells from all batches in roughly the proportions expected from cell-type biology and batch size, not skewed toward one batch.
Red flag
A cluster is >90% from a single batch, or a cell type common in one batch is nearly absent from another, exactly the pattern that produces fake differential accessibility once you compare 'batch A vs batch B' as if it were biology.
06Check whether per-batch peak calls actually agree
Call peaks separately for each batch (or for the same cluster within each batch) and compare the resulting peak sets directly.
r
```
peaks_A <- CallPeaks(atac, group.by = "batch", idents = "batchA", macs2.path = "/path/to/macs2")
peaks_B <- CallPeaks(atac, group.by = "batch", idents = "batchB", macs2.path = "/path/to/macs2")
```
Healthy
Jaccard overlap between batch-specific peak sets for the same presumed cell type is high (well above 0.5 as a rough guide).
Red flag
Peak sets barely overlap, or a rare cell type's peaks show up in only one batch, your consensus peak-by-cell matrix already encodes batch before normalization even runs.
07Run Harmony and verify it mixed batch without erasing cell type
Correct the LSI embedding with Harmony, recluster on the corrected dimensions, and recheck both the batch cross-tab (check 5) and known marker gene activity scores on the new clusters.
r
```
atac <- RunHarmony(atac, group.by.vars = "batch", reduction = "lsi",
                    dims.use = 2:30, reduction.save = "harmony")
atac <- FindNeighbors(atac, reduction = "harmony", dims = 2:30) %>% FindClusters()
atac <- RunUMAP(atac, reduction = "harmony", dims = 2:30)
```
Healthy
Post-Harmony clusters mix batches evenly while marker gene activity/motif scores still cleanly separate known cell types.
Red flag
Batches mix but marker signal flattens too (over-correction blending distinct cell types), or clusters still separate mostly by batch (under-correction).

What to do about it

Correct with Harmony on the right LSI dimensions

When: Batch is not confounded with condition, every biological group has cells in every batch, and check 4 shows a genuine batch axis distinct from cell-type structure.

Run Harmony on LSI dims 2:30 (never include the depth-dominated component 1), grouped by your batch variable, then cluster and build UMAP on the harmony reduction instead of raw LSI.

Caveat: Harmony adjusts only the embedding, never the underlying peak counts. Never run differential accessibility on Harmony coordinates or 'corrected' values, always test on raw counts per cluster after using Harmony only to define the clusters.

Call peaks per cell type, not per batch or pooled

When: Check 6 shows batch-specific peak sets disagreeing, or you suspect rare cell types are missing from a pooled peak set.

Do a first coarse clustering pass (even on uncorrected data), call peaks per cluster/cell type with CallPeaks(group.by = 'celltype'), then merge into one consensus peak set and requantify all cells against it before reclustering.

Caveat: This is circular, you need clusters to call good peaks and good peaks to get reliable clusters. Iterate once: coarse clusters, then peaks, then a final clustering pass on the consensus matrix; don't chase more rounds than that.

Exclude the depth-dominated component(s) explicitly, every time

When: DepthCor (check 3) shows more than just component 1 correlated with depth.

Use dims 2:N (or whatever range clears the depth-correlation threshold) for FindNeighbors, RunUMAP, and Harmony, and rerun DepthCor whenever you change the feature set or add samples.

Caveat: Don't hardcode 'drop component 1' as a rule, verify it on your own data each time, since extreme depth differences between batches can push the artifact into component 2 or beyond.

Redesign, don't correct, when batch and condition are confounded

When: Check 1 shows a batch variable that perfectly predicts biological condition (e.g. all controls sequenced in batch 1, all treated in batch 2).

No embedding correction can separate two variables that vary together in every sample. Rebalance remaining samples across batches if you can still generate data, or explicitly report the confound as a limitation and avoid drawing condition-specific conclusions from that axis.

Caveat: Costs real time and sequencing budget, but reporting 'corrected' results from a confounded design is worse, it produces confident, wrong biology.

Standardize nuclei:Tn5 stoichiometry, and downsample only for specific depth-sensitive comparisons

When: Check 2 shows large, systematic fragments-per-cell differences traced back to tagmentation batch.

Fix the protocol going forward (consistent nuclei count and Tn5 volume per reaction). For the dataset you already have, downsample fragments per cell to the lowest-depth batch's median only when making a specific depth-sensitive comparison, not as your default preprocessing step.

Caveat: Downsampling throws away real signal and can introduce new dropout artifacts of its own; use it narrowly, not as a blanket fix applied to the whole dataset.

When not to "fix" it

If batch corresponds to donor identity and your question is about donor-specific chromatin states rather than a shared cell-type map, don't integrate it away, that variation is the signal, not noise, the same logic that applies to deciding whether to integrate scRNA-seq datasets across donors. And if check 1 shows batch fully confounded with your biological condition, "fixing" it computationally is not a fix: Harmony, ComBat, or any other correction will blend the two indistinguishable sources of variation and hand you a plot that looks clean while hiding an unanswerable comparison. In that case the correct move is redesign or an honest limitation, not a smoother PCA plot.

Five things experienced analysts do here

Record batch metadata (lane, prep date, operator, kit lot) at the bench, not after the fact, you cannot correct for or even detect a batch you never annotated.
Run DepthCor before trusting any LSI embedding, on every new dataset, never assume component 1 is the only depth-correlated axis without checking.
Cross-tabulate batch against biological condition before you touch Harmony; if any cell of that table is empty or near-zero, stop and flag confounding instead of correcting.
Never run differential accessibility on Harmony-corrected embeddings or values, Harmony only adjusts coordinates, so DA testing always goes back to raw peak counts within clusters.
Treat per-batch peak-calling disagreement (Jaccard overlap) as a first-class QC number, not a downstream footnote, it tells you whether your consensus peak set is even valid across samples before you normalize anything.

Questions people ask

How do I check for batch effects in scATAC-seq?: Color your LSI or UMAP embedding by known technical variables, lane, prep date, operator, not just by condition, and look for a dominant axis that tracks those instead of cell type. Back that up with DepthCor() to see which LSI components correlate with per-cell fragment counts, and cross-tabulate cluster membership by batch to catch skewed proportions numerically, not just visually.
Is a batch effect in scATAC-seq the same problem as in bulk or single-cell RNA-seq?: The underlying idea is the same, systematic technical variation masquerading as biology, but the mechanism differs. scATAC-seq batch effects are driven heavily by nuclei:Tn5 stoichiometry at tagmentation, and the peak-by-cell matrix is sparser and closer to binary than a gene-by-cell matrix, so depth swamps LSI component 1 more severely than it swamps PC1 in RNA-seq. Peak sets themselves can also disagree between batches, a failure mode fixed-annotation RNA-seq doesn't have.
Should I use Harmony or ComBat for scATAC-seq batch correction?: Harmony is the standard choice because it operates directly on the LSI embedding that Signac or ArchR already produce, iteratively adjusting coordinates per batch without needing a peak-level statistical model. ComBat-style methods built for continuous, roughly normal bulk data don't map cleanly onto sparse, near-binary peak counts. Note Harmony corrects only the embedding, not the peak matrix; for cases where linear methods like Harmony or MNN fail outright, epiConv is a documented alternative that skips the embedding step entirely.
When should I not correct for batch effects in scATAC-seq?: When batch is confounded with your biological condition, every control sample in one batch, every treatment sample in another, no algorithm, Harmony included, can separate the two. Forcing a correction hides the confound rather than resolving it. The honest fix is redesigning the experiment or explicitly reporting the confound as a limitation, not running a batch-correction tool and trusting the output.
What does a batch effect look like on a scATAC-seq PCA or LSI plot?: A dominant component, frequently LSI2, after you've already dropped the depth-driven LSI1, separates cells by sequencing date, lane, or operator rather than by cell type. DepthCor often shows that same axis still correlating with fragment counts, meaning what looks like 'batch' is frequently sequencing depth wearing a batch label.

Guide · How to Avoid Pseudoreplication in Single-Cell ATAC-seq
Guide · How to Call Peaks You Can Trust in Single-Cell ATAC-seq
Guide · How to Tell If You Overclustered in Single-Cell ATAC-seq
Guide · How to Detect Batch Effects in Spatial Transcriptomics
Glossary · Peak calling
Glossary · scATAC-seq

Sources

Joint analysis of scATAC-seq datasets using epiConv — Nuclei:Tn5 stoichiometry driving depth differences; limits of linear correction methods like Harmony on LSI
Signac Peak Calling Vignette — CallPeaks() usage for per-cell-type vs whole-dataset peak calling
GitHub Discussion: DepthCor and LSI component 1 correlation with depth — Standard practice of excluding LSI component 1 from downstream analysis
GitHub Discussion: Reference mapping on harmony-corrected LSI embedding — Harmony does not produce corrected loadings for projecting new cells onto the corrected space

Part of the Batch effects series.