Chapter 8: Batch Effect Correction & Multi-Sample Integration

Chapter 8: Batch Effect Correction & Multi-Sample Integration

When analyzing single-cell datasets compiled across multiple sequencing runs, donors, clinical sites, or technology platforms (e.g. 10x vs. Drop-seq), technical differences inevitably mask genuine biological signals. Computational integration projects diverse datasets into a shared coordinate space where shared cell types align without erasing true condition-specific perturbations.

1. Integration Algorithmic Paradigms

1. Mutual Nearest Neighbors (MNN / fastMNN): Identifies pairs of cells across two batches that are mutual nearest neighbors in high-dimensional space, calculating a batch-correction vector to translate one dataset onto the other.
2. Seurat v5 CCA / RPCA: Uses Canonical Correlation Analysis (CCA) or Reciprocal PCA (RPCA) to identify shared correlation structures, identifying cross-batch "anchors" to project datasets into a shared reference.
3. Harmony: Operates directly on low-dimensional PCA embeddings. Harmony iteratively groups cells into soft clusters using Maximum Diversity clustering, computes batch correction offsets, and adjusts cell positions until clusters are evenly mixed.
4. scVI (single-cell Variational Inference): A deep generative model using variational autoencoders (VAEs) that models gene expression as a zero-inflated negative binomial distribution while explicitly conditioning on batch covariates.

2. Fast Integration with Harmony in Seurat

library(harmony)
library(Seurat)

# Assume 'combined_obj' contains merged datasets with a 'batch' metadata column
combined_obj <- NormalizeData(combined_obj)
combined_obj <- FindVariableFeatures(combined_obj, nfeatures = 2500)
combined_obj <- ScaleData(combined_obj)
combined_obj <- RunPCA(combined_obj, npcs = 30)

# Run Harmony batch correction on PCA coordinates
combined_obj <- RunHarmony(
    combined_obj, 
    group.by.vars = "batch", 
    reduction.use = "pca",
    dims.use = 1:30,
    plot_convergence = FALSE
)

# Run UMAP and clustering directly on the corrected Harmony embedding
combined_obj <- RunUMAP(combined_obj, reduction = "harmony", dims = 1:30)
combined_obj <- FindNeighbors(combined_obj, reduction = "harmony", dims = 1:30)
combined_obj <- FindClusters(combined_obj, resolution = 0.6)

# Verify mixing: cells from different batches should interleave within cell type clusters
DimPlot(combined_obj, reduction = "umap", group.by = c("batch", "seurat_clusters"))

3. Deep Generative Integration with scVI in Python

import scvi
import scanpy as sc

# Set up AnnData with raw counts
adata_raw = adata.copy()
adata_raw.X = adata.layers['counts'] # scVI requires unnormalized raw counts!

# Register AnnData with scVI
scvi.model.SCVI.setup_anndata(
    adata_raw, 
    layer="counts", 
    batch_key="batch_id", 
    categorical_covariate_keys=["donor_id"]
)

# Train the Variational Autoencoder
vae = scvi.model.SCVI(adata_raw, n_layers=2, n_latent=30, gene_likelihood="nb")
vae.train(max_epochs=250, early_stopping=True)

# Extract batch-corrected 30-dimensional latent space
adata.obsm["X_scVI"] = vae.get_latent_representation()

# Cluster and visualize on scVI representation
sc.pp.neighbors(adata, use_rep="X_scVI")
sc.tl.umap(adata)
sc.tl.leiden(adata, resolution=0.6)
sc.pl.umap(adata, color=["batch_id", "leiden"])