Bio Single Cell Batch Integration

Version Compatibility

Reference examples tested with: scanpy 1.10+, Seurat 5.0+, scvi-tools 1.1+, harmonypy 0.0.10+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures
• R: packageVersion('<pkg>') then ?function_name to verify parameters

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Batch Integration

"Integrate my batches" -> Learn a shared low-dimensional representation that mixes technical batches while preserving biological cell states, then cluster and visualize on it.

• Python: sce.pp.harmony_integrate, scvi.model.SCVI, sce.pp.bbknn, scanorama
• R: RunHarmony, IntegrateLayers (Seurat v5), fastMNN (batchelor)

Governing Principle

Integration trades batch-mixing against biological-signal preservation, and the two cannot be jointly maximized. The algorithm removes variance along directions where batches differ, on the assumption that cell-type composition is shared across batches; it does not "know" which variance is technical. When batch correlates with a real biological axis, the method cannot distinguish them and, by construction, erases biology - this is information-theoretic, not a tuning problem. Over-correction is the silent failure: it absorbs rare cell types into common neighbors, snaps continuous gradients toward shared anchors, and deletes condition-specific populations, all while batch-mixing metrics improve. The cells most worth finding (rare, transitional, novel) are exactly the ones integration most endangers. Batch and biology are unidentifiable under confounding - the only tie-breaker is external information (shared controls, multiplexed designs, known shared types), and the fix for a confounded design is experimental, not computational. Always keep the uncorrected embedding for before/after comparison, and never run differential expression on batch-corrected expression.

When NOT to Integrate

Visualize the uncorrected data first; integration is a bias-variance trade and removing batch variance risks removing biology correlated with batch.

• Confounded design (each condition is its own batch, e.g. all controls day 1, all treated day 2): no algorithm can separate batch from biology. The diagnostic: cluster the uncorrected data and cross-tabulate clusters x batch x condition; pure-by-batch clusters that are also condition-aligned mean integration is unsafe. The fix is experimental - multiplex conditions across batches (cell hashing; genetic demux via souporcell/vireo; split each condition across capture days).
• Technical replicates of the same tissue that already mix well: over-correction risk outweighs benefit.
• Per-sample analyses (CNV/tumor-clone inference): integration would erase the signal of interest.

Over-correction signatures: rare types collapsing into neighbors, lost known gradients, disappearing condition-specific populations, markers no longer separating known cell types.

Method Selection

No method wins universally - the scIB benchmark (Luecken 2022, 68 method/preprocessing combos) scores integration as overall = 0.6 x bio-conservation + 0.4 x batch-removal, deliberately weighting biology higher because erasing it is worse than imperfect mixing. Methodology evolves; verify current best practice and APIs against the installed package docs, and in practice run 2-3 candidates and score them (see Evaluating Integration).

scIB headline: top combined performers were scANVI, scVI, Scanorama, scGen; Harmony and Seurat were strong on simpler tasks with the best usability; BBKNN sits at the batch-removal end. "Deep methods are always best" is not supported - Harmony/Seurat win simple/small tasks; deep methods win complex/large/label-rich tasks.

Strength Parameters

Aggressive settings increase mixing and over-correction risk in lockstep - raise correction strength only after confirming under-correction, and re-check rare populations after each change.

Integrate with Harmony

Goal: Correct batch in PCA space and run downstream steps on the corrected embedding. Approach: Joint preprocessing -> PCA -> Harmony -> neighbors/UMAP/clustering on X_pca_harmony.

import scanpy as sc
import scanpy.external as sce

adata = sc.read_h5ad('merged.h5ad')
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=2000, batch_key='batch')
adata.raw = adata
adata = adata[:, adata.var.highly_variable]
sc.pp.scale(adata, max_value=10)
sc.tl.pca(adata, n_comps=50)

sce.pp.harmony_integrate(adata, key='batch')  # writes adata.obsm['X_pca_harmony']
sc.pp.neighbors(adata, use_rep='X_pca_harmony')
sc.tl.umap(adata)
sc.tl.leiden(adata, flavor='igraph', n_iterations=2, directed=False)

In Seurat: RunHarmony(obj, group.by.vars = 'orig.ident', reduction.use = 'pca') writes a harmony reduction; group.by.vars takes a vector to correct multiple covariates.

Integrate with scVI / scANVI

Goal: Learn a batch-invariant latent space from raw counts, optionally protecting known labels. Approach: Put raw counts in a layer, register batch (and labels for scANVI), train, and use the latent embedding downstream.

import scvi
import scanpy as sc

adata = sc.read_h5ad('merged.h5ad')
adata.layers['counts'] = adata.X.copy()  # scVI needs raw counts
sc.pp.highly_variable_genes(adata, n_top_genes=2000, flavor='seurat_v3',
                            layer='counts', batch_key='batch')
adata = adata[:, adata.var.highly_variable].copy()

scvi.model.SCVI.setup_anndata(adata, layer='counts', batch_key='batch')
model = scvi.model.SCVI(adata, n_latent=10, gene_likelihood='zinb')
model.train()  # default max_epochs heuristic scales down for large data
adata.obsm['X_scVI'] = model.get_latent_representation()

scanvi = scvi.model.SCANVI.from_scvi_model(model, 'Unknown', labels_key='cell_type')
scanvi.train(max_epochs=20)
adata.obs['scanvi_label'] = scanvi.predict()

The scVI latent space is not "biology with batch removed": it is a learned nonlinear embedding optimized to reconstruct counts while being marginally independent of batch. Its dimensions are entangled, individually uninterpretable, and carry no guaranteed correspondence to any biological quantity - treat it as a coordinate system for neighbors/clustering, not a measurement. Note unlabeled_category ('Unknown') is the second positional argument to from_scvi_model, before labels_key.

Integrate with Seurat v5

Goal: Use Seurat v5's modular layer-based integration with a chosen method. Approach: Split layers by batch, run the standard pipeline, call IntegrateLayers, rejoin.

library(Seurat)

merged[['RNA']] <- split(merged[['RNA']], f = merged$batch)
merged <- NormalizeData(merged)
merged <- FindVariableFeatures(merged)
merged <- ScaleData(merged)
merged <- RunPCA(merged)

merged <- IntegrateLayers(merged, method = RPCAIntegration,
                          orig.reduction = 'pca', new.reduction = 'integrated.rpca')
merged <- JoinLayers(merged)
merged <- FindNeighbors(merged, reduction = 'integrated.rpca', dims = 1:30)
merged <- FindClusters(merged, resolution = 0.5)
merged <- RunUMAP(merged, reduction = 'integrated.rpca', dims = 1:30)

Methods are passed as bare symbols: CCAIntegration, RPCAIntegration, HarmonyIntegration, FastMNNIntegration, scVIIntegration. For graph-only correction with BBKNN in Python: sce.pp.bbknn(adata, batch_key='batch') rewrites the neighbor graph in place (very fast, feeds Leiden/UMAP only).