Integrating Datasets Across Batches

Integrating Datasets Across Batches

Single-cell RNA-seq experiments are rarely generated in a single run. Samples are processed on different days, by different technicians, using different library preparation kits, or sequenced on different instruments. Each of these sources of variation introduces batch effects : systematic technical differences in gene expression measurements that are unrelated to true biological differences between cells. If left uncorrected, batch effects can dominate a UMAP plot, causing cells to cluster by sequencing run rather than by cell type, and leading to biologically meaningless conclusions. Batch correction is not simply a matter of subtracting out unwanted variation. Overly aggressive correction can erase genuine biological signal — for example, removing a real disease-associated transcriptional state because it happens to be enriched in one batch. The goal of dataset integration is therefore to remove technical variation while preserving biological variation . Modern integration methods accomplish this by first identifying pairs or groups of similar cells across datasets (often in a shared low-dimensional space), then learning a correction that aligns these matched populations without forcing unrelated cell types together. Two widely used approaches illustrate the concept. Mutual nearest neighbor (MNN)-based methods , such as the integration pipeline in Seurat, identify cells in one dataset whose nearest neighbors in another dataset are also mutually closest, and use these anchor pairs to estimate and remove batch-specific shifts in expression space. Graph- and factorization-based methods , such as Harmony, perform iterative clustering in which cells are assigned to clusters, cluster-level batch effects are estimated and subtracted, and the process repeats until the embedding is balanced across batches. Both operate on the top principal components rather than raw counts, which reduces noise and computational cost. A critical point for pharmacy and biotech applications: integration changes the expression matrix used for visualization and clustering , but downstream analyses such as differential expression should generally be performed on the original, uncorrected counts using batch (or sample) as a covariate. Corrected values are useful for computing neighborhoods and embeddings, but they are not reliable for quantitative gene-level statistics. The first example integrates two Seurat objects representing control and treated PBMC samples processed in separate batches: library(Seurat) # Load two separately processed objects ctrl <- readRDS("pbmc_control.rds") trt <- readRDS("pbmc_treated.rds") # Standard per-dataset preprocessing ctrl <- NormalizeData(ctrl) |> FindVariableFeatures() trt <- NormalizeData(trt) |> FindVariableFeatures() # Identify integration anchors and transfer the correction anchors <- FindIntegrationAnchors( object.list = list(ctrl, trt), dims = 1:30, reduction = "rpca" # robust to strong batch differences ) integrated <- IntegrateData( anchorset = anchors, normalization.method = "LogNormalization" ) # Cluster the integrated data integrated <- ScaleData(integrated, verbose = FALSE) |> RunPCA() |> RunUMAP(dims = 1:30) |> FindNeighbors() |> FindClusters(resolution = 0.5) # Check that batches are now mixed within clusters DimPlot(integrated, group.by = "batch") The second example performs the equivalent integration in Python with Scanpy and Harmony, a common choice for large-scale studies: import scanpy as sc import scanpy.external as sce import anndata as ad # Concatenate batches into a single AnnData object adata_ctrl = sc.read_h5ad("pbmc_control.h5ad") adata_trt = sc.read_h5ad("pbmc_treated.h5ad") adata = ad.concat([adata_ctrl, adata_trt], label="batch") # Standard preprocessing on the combined object 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") sc.pp.scale(adata, max_value=10) sc.tl.pca(adata, n_comps=30) # Harmony corrects batch effects in the PCA space sce.pp.harmony_integrate(adata, key="batch") # Build neighbors and embeddings using the corrected space sc.pp.neighbors(adata, use_rep="X_pca_harmony", n_neighbors=15) sc.tl.umap(adata) sc.tl.leiden(adata, resolution=0.5, key_added="clusters") # Inspect batch mixing on the UMAP sc.pl.umap(adata, color=["batch", "clusters"]) How do you know whether integration worked? Inspect the UMAP colored by batch: cells from different batches should be well mixed within each cluster. Quantitative diagnostics include the kBET statistic, which measures the local batch composition of neighborhoods, and silhouette scores computed on batch labels. If clusters remain batch-segregated, consider increasing the number of variable features, adjusting the number of principal components, or using a more conservative method such as reciprocal PCA anchors. Conversely, if known cell types that should be distinct collapse together, the correction may be too strong. Batch effects are systematic, not random: they arise from processing variables (kit, time, operator, sequencer) and can masquerade as biological signal if uncorrected. Integration aligns shared biology: methods like Seurat anchors and Harmony match similar cells across batches and correct in a reduced dimensional space, preserving cell-type structure. Use corrected data for clustering and visualization; use raw counts for differential expression, modeling batch or sample as a covariate to avoid inflated false positives. Validate integration: check batch mixing on UMAPs and use metrics such as kBET; both under-correction (batch-driven clusters) and over-correction (merged distinct cell types) are common failure modes. Design matters: integration cannot rescue poorly designed experiments — balance conditions across batches (e.g., avoid processing all controls in one run and all treated samples in another) whenever possible. Practice exercise: You receive two single-cell datasets from a drug-perturbation study: one batch of vehicle-treated PBMCs and one batch of compound-treated PBMCs, each with a "batch" metadata column. Using either Seurat or Scanpy/Harmony, integrate the two datasets, cluster the integrated data, and produce a UMAP colored by batch and by cluster. Then answer the following: (1) Are batches well mixed within clusters, or do any clusters appear batch-specific? (2) If a cluster contains cells almost exclusively from the treated batch, how would you distinguish a genuine drug-induced cell state from residual batch effect? (3) Write the differential expression comparison you would run to identify drug-responsive genes, and explain why you would not use the integrated expression matrix for that test.