Batch Integration
Integrate multiple scRNA-seq datasets to remove batch effects while preserving biological variation.
Tool Comparison
Tool Speed Scalability Best For
Harmony Fast Good Quick integration, most use cases
scVI Moderate Excellent Large datasets, deep learning
Seurat CCA/RPCA Moderate Good Conserved biology across batches
fastMNN Fast Good MNN-based correction
Harmony (R/Python)
R with Seurat
library(Seurat) library(harmony)
Merge datasets first
merged <- merge(sample1, y = list(sample2, sample3), add.cell.ids = c('S1', 'S2', 'S3'))
Standard preprocessing
merged <- NormalizeData(merged) merged <- FindVariableFeatures(merged) merged <- ScaleData(merged) merged <- RunPCA(merged)
Run Harmony on PCA embeddings
merged <- RunHarmony(merged, group.by.vars = 'orig.ident', dims.use = 1:30)
Use harmony embeddings for downstream
merged <- RunUMAP(merged, reduction = 'harmony', dims = 1:30) merged <- FindNeighbors(merged, reduction = 'harmony', dims = 1:30) merged <- FindClusters(merged, resolution = 0.5)
Multiple Batch Variables
Correct for both sample and technology
merged <- RunHarmony(merged, group.by.vars = c('sample', 'technology'), dims.use = 1:30, max.iter.harmony = 20)
Python with Scanpy
import scanpy as sc import scanpy.external as sce
adata = sc.read_h5ad('merged.h5ad')
Standard preprocessing
sc.pp.normalize_total(adata, target_sum=1e4) sc.pp.log1p(adata) sc.pp.highly_variable_genes(adata, batch_key='batch') adata = adata[:, adata.var.highly_variable] sc.pp.scale(adata) sc.tl.pca(adata)
Run Harmony
sce.pp.harmony_integrate(adata, key='batch')
Use corrected embedding
sc.pp.neighbors(adata, use_rep='X_pca_harmony') sc.tl.umap(adata) sc.tl.leiden(adata)
scVI (Python)
Deep learning-based integration for large datasets.
import scvi import scanpy as sc
adata = sc.read_h5ad('merged.h5ad')
Setup for scVI
scvi.model.SCVI.setup_anndata(adata, batch_key='batch')
Train model
model = scvi.model.SCVI(adata, n_latent=30, n_layers=2) model.train(max_epochs=100, early_stopping=True)
Get latent representation
adata.obsm['X_scVI'] = model.get_latent_representation()
Use for downstream
sc.pp.neighbors(adata, use_rep='X_scVI') sc.tl.umap(adata) sc.tl.leiden(adata)
scVI with Covariates
Include continuous covariates
scvi.model.SCVI.setup_anndata(adata, batch_key='batch', continuous_covariate_keys=['percent_mito'])
model = scvi.model.SCVI(adata, n_latent=30) model.train()
scANVI (with cell type labels)
If you have reference labels for some cells
scvi.model.SCANVI.setup_anndata(adata, batch_key='batch', labels_key='cell_type', unlabeled_category='Unknown')
model = scvi.model.SCANVI(adata, n_latent=30) model.train(max_epochs=100)
Predict labels for unlabeled cells
adata.obs['predicted_type'] = model.predict()
Seurat Integration (R)
CCA-based Integration
library(Seurat)
Split by batch
obj_list <- SplitObject(merged, split.by = 'batch')
Normalize each
obj_list <- lapply(obj_list, function(x) { x <- NormalizeData(x) x <- FindVariableFeatures(x, selection.method = 'vst', nfeatures = 2000) return(x) })
Find integration anchors
anchors <- FindIntegrationAnchors(object.list = obj_list, dims = 1:30)
Integrate
integrated <- IntegrateData(anchorset = anchors, dims = 1:30)
Switch to integrated assay for downstream
DefaultAssay(integrated) <- 'integrated' integrated <- ScaleData(integrated) integrated <- RunPCA(integrated) integrated <- RunUMAP(integrated, dims = 1:30)
RPCA (Faster for Large Datasets)
Use reciprocal PCA for faster integration
anchors <- FindIntegrationAnchors(object.list = obj_list, dims = 1:30, reduction = 'rpca') integrated <- IntegrateData(anchorset = anchors, dims = 1:30)
Seurat v5 Integration
Seurat v5 uses layers
merged[['RNA']] <- split(merged[['RNA']], f = merged$batch) merged <- IntegrateLayers(merged, method = CCAIntegration, orig.reduction = 'pca', new.reduction = 'integrated.cca') merged <- JoinLayers(merged)
fastMNN (R)
library(batchelor) library(SingleCellExperiment)
Convert Seurat to SCE
sce <- as.SingleCellExperiment(merged)
Run fastMNN
corrected <- fastMNN(sce, batch = sce$batch, d = 30, k = 20)
Extract corrected values
reducedDim(sce, 'MNN') <- reducedDim(corrected, 'corrected')
Evaluate Integration
Mixing Metrics (R)
LISI score (lower = more mixed)
library(lisi) lisi_scores <- compute_lisi(Embeddings(merged, 'harmony'), merged@meta.data, c('batch', 'cell_type'))
Batch mixing should be high, cell type separation preserved
mean(lisi_scores$batch) # Want high mean(lisi_scores$cell_type) # Want low (preserved)
Visual Assessment
Before integration
DimPlot(merged, reduction = 'pca', group.by = 'batch') DimPlot(merged, reduction = 'pca', group.by = 'cell_type')
After integration
DimPlot(merged, reduction = 'harmony', group.by = 'batch') DimPlot(merged, reduction = 'harmony', group.by = 'cell_type')
Silhouette Score (Python)
from sklearn.metrics import silhouette_score
Batch silhouette (want low - batches mixed)
batch_sil = silhouette_score(adata.obsm['X_scVI'], adata.obs['batch'])
Cell type silhouette (want high - types separated)
celltype_sil = silhouette_score(adata.obsm['X_scVI'], adata.obs['cell_type'])
Complete Workflow
library(Seurat) library(harmony)
Load and merge samples
samples <- list.files('data/', pattern = '*.h5', full.names = TRUE) obj_list <- lapply(samples, Read10X_h5) names(obj_list) <- gsub('.h5', '', basename(samples))
merged <- merge(CreateSeuratObject(obj_list[[1]], project = names(obj_list)[1]), y = lapply(2:length(obj_list), function(i) CreateSeuratObject(obj_list[[i]], project = names(obj_list)[i])))
QC
merged[['percent.mt']] <- PercentageFeatureSet(merged, pattern = '^MT-') merged <- subset(merged, nFeature_RNA > 200 & nFeature_RNA < 5000 & percent.mt < 20)
Preprocess
merged <- NormalizeData(merged) merged <- FindVariableFeatures(merged, nfeatures = 2000) merged <- ScaleData(merged, vars.to.regress = 'percent.mt') merged <- RunPCA(merged, npcs = 50)
Integrate with Harmony
merged <- RunHarmony(merged, group.by.vars = 'orig.ident')
Downstream analysis on integrated data
merged <- RunUMAP(merged, reduction = 'harmony', dims = 1:30) merged <- FindNeighbors(merged, reduction = 'harmony', dims = 1:30) merged <- FindClusters(merged, resolution = 0.5)
DimPlot(merged, group.by = c('orig.ident', 'seurat_clusters'), ncol = 2)
When to Use Each Method
Scenario Recommended
Quick integration, most cases Harmony
Large datasets (>500k cells) scVI or Harmony
Strong batch effects scVI
Reference mapping Seurat anchors or scANVI
Preserving rare populations fastMNN
Related Skills
-
single-cell/preprocessing - QC before integration
-
single-cell/clustering - Clustering after integration
-
single-cell/cell-annotation - Annotation after integration
-
single-cell/multimodal-integration - Multi-omic integration (different from batch)