How do I check for a batch effect in snRNA-seq?

Cross-tabulate batch against condition, then color a PCA by freezing date, prep day, lane and condition. If PC1 separates a technical variable instead of the phenotype, batch is likely driving variance. Add per-sample QC, intronic read ratio and a sample-to-sample distance heatmap to confirm.

When should I correct a batch effect in snRNA-seq?

Correct for clustering and annotation when batch is crossed with condition and you want a shared map. For testing expression, model batch in a pseudobulk design rather than correcting the matrix. If batch and condition are fully confounded, do not correct: the answer is a redesign.

Can I use ComBat for snRNA-seq batch correction?

It is not a good default for nuclei data. A benchmark found ComBat inflated gene expression, and it does not address ambient RNA differences between preps. Harmony on the embedding, or batch as a covariate in a pseudobulk model, are safer choices.

Does Harmony change my expression values?

No. Harmony iteratively adjusts the PCA embedding of each cell with respect to batch covariates and leaves the original expression untouched. Use the corrected embedding for neighbors, clusters and UMAP, and run DE on the original counts with batch in the model.

Why is batch more of a problem in snRNA-seq than scRNA-seq?

Nuclear isolation releases cytoplasmic material, so snRNA-seq carries more ambient RNA, and the amount varies by prep and freezing date. That adds a technical axis that is mostly absent in whole-cell data, and integration tools may align it away together with real biology.

Sanity check · Single-Nucleus RNA-seq

How to Detect Batch Effects in Single-Nucleus RNA-seq

Before you run Harmony on your frozen-tissue nuclei, find out whether the batch axis is technical noise or the condition you came to study.

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

You have nuclei from frozen tissue, processed on several dates, and the UMAP shows islands. Some islands are cell types. Others look suspiciously like freezing date or prep day. Now you need to decide whether to run Harmony, scVI or Seurat integration, and whether doing so will rescue the comparison or erase it.

What is at stake: nuclei preps from different freezing dates differ in ambient RNA fraction and intronic read ratio. Integration tools align those differences away, and if your condition tracks the freezing date, they align away the condition effect with it. You end up with a clean-looking UMAP and no biology left to test.

This page gives you an order of checks for the next hour. Start with the cheap ones (PCA colored by batch, a sample table). Then use the snRNA-seq-specific ones (ambient and intronic signal). End with a decision: model batch, correct it, or admit the design cannot answer the question.

What it looks like when it's happening

PC1 or PC2 in the PCA separates samples by freezing date, prep day, lane or operator instead of by condition.
UMAP clusters that contain nuclei from only one or two samples, while the same cell type from other samples sits somewhere else.
Glia (for example oligodendrocytes or astrocytes) express highly expressed neuronal genes, and a cluster annotated as 'immature oligodendrocytes' is mostly ambient-contaminated glia.
Per-sample violin plots of genes detected, UMI counts or intronic read fraction show one batch sitting clearly higher or lower than the others.
The sample-to-sample distance heatmap groups samples by batch label, not by condition label.
A cross-tab of batch by condition shows each condition living in only one batch.
Top differentially expressed genes between conditions are ribosomal, mitochondrial or synaptic genes typical of cytoplasmic ambient RNA.

Why it happens

A batch effect is any systematic technical difference between groups of samples: processing date, library prep, sequencing lane, operator, reagent lot, or here, freezing date. In snRNA-seq the largest source is ambient RNA. Nuclear isolation lyses cytoplasm, and that material ends up in every droplet. snRNA-seq is more exposed to this than scRNA-seq, so the ambient fraction can differ a lot between preps even when the tissue is the same.

Ambient RNA comes in two flavors. Extra-nuclear ambient RNA has a lower intronic read ratio and is enriched in ribosomal, mitochondrial and synaptic genes. Nuclear ambient RNA has a higher intronic read ratio and comes from highly expressed neuronal genes. Because the contamination reflects whatever is abundant in the tissue, it is cell-type-specific in its effect: glia in a brain prep can all carry neuronal transcripts, and a contaminated glial population can be mistaken for a new cell type.

Depth and gene properties add to this. Batch effects depend on expression level, transcript length and nucleotide composition, and cell types differ in total mRNA content, so a depth-biased batch can look like a cell-type difference. Nuclei also detect fewer genes and reads than whole cells and are biased toward longer genes, so the technical signal is not the same as in scRNA-seq.

The decision turns on confounding, not on how big the batch effect looks. If every sample of condition A was processed in batch 1 and every sample of condition B in batch 2, batch and biology are inseparable. No method can pull them apart. Integration will either leave both in place or remove both, and you cannot tell which from the UMAP.

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

01Cross-tabulate batch against condition
Before any plot, build a sample table with condition, freezing date, prep day, lane, operator and reagent lot, and cross-tabulate each technical variable against condition. This takes a minute and decides whether correction is even possible.
r
```
meta <- seu@meta.data
# one row per sample, then cross-tab each technical variable vs condition
samples <- unique(meta[, c("sample", "condition", "freeze_date", "prep_day", "lane")])
table(samples$freeze_date, samples$condition)
table(samples$prep_day, samples$condition)
```
Healthy
Each batch contains samples from more than one condition, ideally in comparable numbers, so batch and condition can be estimated separately.
Red flag
A condition appears in only one batch, or batch and condition columns line up perfectly. Stop here: no correction method can fix this and the honest answer is a redesign.
02Color PCA by batch variables and by condition
Run PCA on the normalized, variable-gene data and color PC1 vs PC2 (and the next few PCs) by freezing date, prep day, lane, operator, reagent lot and condition, one plot per variable. Compare which variable separates samples along the leading PCs. For pseudobulk-level PCA, aggregate counts per sample first.
r
```
seu <- NormalizeData(seu) |> FindVariableFeatures() |> ScaleData() |> RunPCA()
for (v in c("condition", "freeze_date", "prep_day", "lane")) {
  print(DimPlot(seu, reduction = "pca", group.by = v) + ggtitle(v))
}
ElbowPlot(seu)
```
Healthy
Leading PCs separate cell types and, where relevant, condition. Technical variables are mixed within each cell type.
Red flag
PC1 separates freezing date or prep day rather than phenotype. That means batch likely dominates the variance and may confound the comparison.
03Check per-sample QC distributions, with nuclei-appropriate expectations
Plot nFeature_RNA, nCount_RNA and percent mitochondrial reads per sample as violins. In nuclei, mitochondrial content should be near absent, so low mito% is the norm and an elevated value points to a nuclei isolation problem, not dying cells. Do not apply scRNA-seq thresholds unchanged: nuclei detect fewer genes and reads, so cutoffs need to sit lower. The sources do not give a numeric cutoff, so set it from your own distributions.
r
```
seu[["percent.mt"]] <- PercentageFeatureSet(seu, pattern = "^MT-")
VlnPlot(seu, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"),
        group.by = "sample", pt.size = 0, ncol = 3)
```
Healthy
Samples overlap in genes detected and UMI counts, with mito% low everywhere.
Red flag
One batch has systematically different depth or a mito% well above the rest. Depth-imbalanced batches produce batch clusters whose main axis is depth, and a high-mito sample points to a bad isolation.
04Measure the intronic read ratio per sample
Compute the fraction of intronic (unspliced) reads per nucleus, for example from spliced and unspliced layers produced by your counting pipeline, then plot it by sample and batch. Use it to tell the two ambient types apart: low intronic ratio points to extra-nuclear ambient RNA, high intronic ratio to nuclear ambient RNA. The sources give no cutoff, so compare samples against each other.
r
```
# assumes 'unspliced' and 'spliced' count assays from your quantification
seu$intron_frac <- Matrix::colSums(GetAssayData(seu, assay = "unspliced", layer = "counts")) /
  (Matrix::colSums(GetAssayData(seu, assay = "unspliced", layer = "counts")) +
   Matrix::colSums(GetAssayData(seu, assay = "spliced", layer = "counts")))
VlnPlot(seu, "intron_frac", group.by = "sample", pt.size = 0)
```
Healthy
Intronic fraction is similar across samples and batches, with spread within each sample.
Red flag
Intronic fraction clusters by freezing date or prep day. If it also tracks condition, the batch and the condition are entangled in a way integration will blur.
05Look for ambient signatures in the wrong cell types
Check whether highly expressed neuronal genes (or the dominant genes of your tissue) show up in glia or other non-neuronal clusters, and whether ribosomal, mitochondrial or synaptic genes rank high in cluster markers. Plot the top ambient genes per cluster, split by sample. Run a decontamination tool such as CellBender, DecontX, SoupX or scAR if the signal is strong, then re-check.
Healthy
Marker genes for each cluster are specific to that cell type, and ambient genes show similar low levels across clusters.
Red flag
Glia carry neuronal transcripts, an odd 'immature' cluster is mostly ambient-contaminated cells, or the ambient signal varies by batch.
06Draw a sample-to-sample distance heatmap
Aggregate counts per sample (pseudobulk, ideally within one cell type), normalize, compute sample-to-sample distances and plot a heatmap annotated by batch variables and condition.
r
```
pb <- AggregateExpression(seu, group.by = "sample", assays = "RNA")$RNA
logpb <- log1p(sweep(pb, 2, colSums(pb), "/") * 1e6)
d <- dist(t(logpb))
ann <- unique(seu@meta.data[, c("sample", "condition", "freeze_date")])
rownames(ann) <- ann$sample
pheatmap::pheatmap(as.matrix(d), annotation_row = ann[, c("condition", "freeze_date")],
                   annotation_col = ann[, c("condition", "freeze_date")])
```
Healthy
Samples group by condition, or by a mix, with batch labels scattered across the tree.
Red flag
Samples group by batch. If batch explains the major axes you need to model it, and if it is confounded with condition the answer is a redesign.
07Check UMAP mixing within cell types, before and after correction
Color the UMAP by sample and batch and judge mixing within each annotated cell type, not across the whole map. After correcting, compare against the uncorrected result using the published metrics: k-nearest-neighbor rank displacement for neighborhood preservation, Adjusted Rand Index for cluster stability, and differential expression fidelity (do the same genes stay significant?). Large changes on a dataset with little batch effect signal overcorrection.
Healthy
Each cell type mixes samples after correction, while neighborhoods and clusters stay mostly stable and DE results change modestly.
Red flag
Neighborhoods rearrange heavily, clusters merge or split, or condition-specific states disappear after correction. The method has likely removed biology along with batch.

What to do about it

Redesign or add samples when batch and condition are confounded

When: The cross-tab shows each condition entirely in its own batch, or every sample of one condition shares a freezing date.

Re-run a subset of samples so that each batch contains both conditions, randomizing sample assignment to prep day, lane and operator. If you cannot re-run, report the comparison as exploratory and say plainly that condition and batch cannot be separated.

Caveat: Costs tissue, money and time, and some precious frozen samples cannot be repeated. A candid limitation is still better than a clean UMAP that hides the problem.

Model batch as a covariate in the test, not in the matrix

When: Batch is partly crossed with condition and you want differential expression between conditions.

Aggregate counts to pseudobulk per sample within each cell type and include batch in the design formula of your DE model. The sample, not the nucleus, is the unit of replication. Do not pre-correct the count matrix and then test.

Caveat: Needs enough samples per condition per batch to estimate both terms. With few samples the model will have little power.

Run Harmony on the embedding for clustering and annotation

When: Batch is crossed with condition and your goal is a shared map to annotate cell types, not to test expression.

Run Harmony on the PCA embedding with the technical variables as covariates, for example RunHarmony(seuratObject, c('dataset', 'donor', 'batch_id')), then build the neighbor graph and UMAP from the Harmony embedding. Harmony adjusts the embedding only and never the expression values, so run DE on the original counts. In one benchmark it was the only method that performed consistently well on neighborhood preservation, cluster stability and DE fidelity.

Caveat: The sources do not state the default theta, so check the documentation for your version. Do not include condition as a covariate. Correcting for a variable that tracks condition removes the effect you want.

Use scVI or Seurat v5 CCA integration for larger shifts

When: Harmony leaves visible batch structure within cell types, for example across very different preps.

With scVI, pass batch as the covariate and use the corrected latent space for clustering; the sysVI extension is designed for substantial effects including cell versus nucleus data. In Seurat v5, use IntegrateLayers(object, method = CCAIntegration, orig.reduction = 'pca', new.reduction = 'integrated.cca'). Compare the result against the uncorrected map with the neighborhood and cluster stability metrics.

Caveat: Benchmarks report that MNN, scVI and LIGER substantially altered neighborhood structure and that ComBat inflated expression, even with minimal batch effects. Stronger methods overcorrect more easily.

Remove ambient RNA before integrating

When: Batches differ mainly in ambient fraction or intronic ratio rather than in cell state.

Run a decontamination tool such as CellBender (recommended in the ambient RNA study), DecontX, SoupX, scAR or others on each sample's raw droplet data, then redo QC and the PCA check. Confirm contaminating genes drop in non-neuronal clusters.

Caveat: Decontamination is itself a model and can strip real low-level expression. Check marker genes before and after.

When not to "fix" it

Do not integrate when your question is about a shifted or new cell state, or when batch is confounded with condition. Integration smooths exactly the variation you are trying to find, and for confounded designs it cannot tell the two apart. Counting B cells across donors is an integration job, finding activated B cells in disease is not. Also leave the data alone if the leading PCs track cell type and condition and batch is mixed within each cell type: a small, well-mixed batch term is better handled by including batch in the DE model than by reshaping the embedding. Finally, a sample that differs because the biology differs (a real change in cell-type proportions between conditions, for instance) should not be 'fixed' just because it separates on a UMAP.

Five things experienced analysts do here

Record freezing date, prep day, lane, operator and reagent lot for every sample the day you generate them. You cannot test for a batch variable you did not write down.
Fill in the batch-by-condition table before you look at a single UMAP. If the table shows full confounding, nothing else on this page changes the verdict.
Judge mixing within each cell type, not across the whole UMAP. A mixed global map can hide a batch-specific subcluster.
Run DE on pseudobulk with batch in the design and on uncorrected counts. Treat embedding-level correction as a tool for clustering and annotation only.
Compute the intronic read ratio and an ambient gene score per sample early. They tell you whether a batch is a prep-quality problem before you spend time tuning integration parameters.

Questions people ask

How do I check for a batch effect in snRNA-seq?: Cross-tabulate batch against condition, then color a PCA by freezing date, prep day, lane and condition. If PC1 separates a technical variable instead of the phenotype, batch is likely driving variance. Add per-sample QC, intronic read ratio and a sample-to-sample distance heatmap to confirm.
When should I correct a batch effect in snRNA-seq?: Correct for clustering and annotation when batch is crossed with condition and you want a shared map. For testing expression, model batch in a pseudobulk design rather than correcting the matrix. If batch and condition are fully confounded, do not correct: the answer is a redesign.
Can I use ComBat for snRNA-seq batch correction?: It is not a good default for nuclei data. A benchmark found ComBat inflated gene expression, and it does not address ambient RNA differences between preps. Harmony on the embedding, or batch as a covariate in a pseudobulk model, are safer choices.
Does Harmony change my expression values?: No. Harmony iteratively adjusts the PCA embedding of each cell with respect to batch covariates and leaves the original expression untouched. Use the corrected embedding for neighbors, clusters and UMAP, and run DE on the original counts with batch in the model.
Why is batch more of a problem in snRNA-seq than scRNA-seq?: Nuclear isolation releases cytoplasmic material, so snRNA-seq carries more ambient RNA, and the amount varies by prep and freezing date. That adds a technical axis that is mostly absent in whole-cell data, and integration tools may align it away together with real biology.

Guide · How to Detect Batch Effects in Single-Cell RNA-seq
Guide · How to Detect Batch Effects in Spatial Transcriptomics
Guide · How to Avoid Pseudoreplication in Single-Nucleus RNA-seq
Guide · How to Choose Cell QC Thresholds in Single-Nucleus RNA-seq
Guide · How to Detect Batch Effects in Single-Cell ATAC-seq
Glossary · Batch correction (ComBat, limma removeBatchEffect)
Glossary · Principal component analysis (PCA)
Glossary · Count matrix

Sources

Batch correction methods used in single-cell RNA sequencing analyses are often poorly calibrated — Harmony performed consistently; ComBat inflation and overcorrection by MNN, scVI and LIGER; evaluation metrics.
Ambient RNA analysis reveals misinterpreted and masked cell types in brain single-nuclei datasets — Two types of ambient RNA by intronic ratio, false cell types, decontamination tools.
Comparative Analysis of Single-Nucleus and Single-Cell RNA Sequencing in Human Bone Marrow Mononuclear Cells: Methodological Insights and Trade-offs — snRNA-seq is more susceptible to ambient RNA, lower genes per nucleus, HVG stability.
Batch alignment of single-cell transcriptomics data using deep metric learning — Harmony adjusts the embedding, not expression values.
Harmony on GitHub — RunHarmony with multiple batch covariates.
Seurat Integration Introduction — IntegrateLayers with CCAIntegration in Seurat v5.

Part of the Batch effects series.