Batch Correction & Data Integration

When to Use

• Combining scRNA-seq data from multiple samples/batches
• Removing technical batch effects while preserving biological variation
• Comparing integration methods quantitatively
• Preparing integrated datasets for joint clustering and analysis

Do NOT use when you need to:

• 只有单个样本/批次 → no integration needed, use scrna-qc-preprocess directly
• 需要差异表达 → (use scrna-marker-genes after integration)
• 数据未预处理 → (use scrna-qc-preprocess first)
• 需要空间对齐而非批次校正 → (use paste/spateo for spatial alignment)

Prerequisites

• Python >= 3.10
• Required packages: scanpy, anndata, harmonypy, scib, numpy, pandas, matplotlib, sklearn; optional: scvi-tools
• Input: h5ad file with batch/sample column in obs, PCA embeddings in obsm["X_pca"]

Core Rules

1. Always run QC preprocessing first — batch correction requires PCA embeddings and normalized data as input.
2. Benchmark multiple methods quantitatively (silhouette_batch, silhouette_bio, graph_connectivity) before selecting the best one.
3. Preserve biological variation: a good integration should have low silhouette_batch (mixed batches) but high silhouette_bio (separated cell types).
4. Use the same n_pcs across all integration methods for fair comparison.
5. Prefer Harmony for quick integration on large datasets; use scVI for complex batch effects with deep learning.

Common Mistakes

• Integrating unpreprocessed data: Always run QC + normalization + PCA before integration — raw count matrices will fail.
• Over-correcting: Too strong correction merges distinct cell types; reduce Harmony theta or scVI epochs.
• Using wrong batch column: Verify the batch column distinguishes technical batches, not biological groups (e.g., don't use condition as batch).

Troubleshooting

Output

写入 --output 指定的 JSON 文件：

{
  "skill": "scrna-batch-correction",
  "status": "success",
  "steps": ["load_data", "pre_integration_viz", "run_methods", "benchmark", "select_best", "save"],
  "result": {
    "input_file": "multi_batch.h5ad",
    "output_file": "integrated.h5ad",
    "n_cells": 80000,
    "n_batches": 4,
    "methods_tested": ["harmony", "scanorama", "bbknn", "scvi"],
    "best_method": "harmony",
    "metrics": {
      "silhouette_batch": 0.05,
      "silhouette_bio": 0.42,
      "composite_score": 0.37
    },
    "metrics_csv": "integration_metrics.csv",
    "plots": ["pre_integration_umap.png", "integration_comparison.png"]
  }
}

Workflow

Step 1: Load Multi-Batch Data

import scanpy as sc
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.metrics import silhouette_score, adjusted_rand_score

adata = sc.read_h5ad("multi_batch.h5ad")

# Identify batch column
batch_col = None
for col in ["batch", "sample", "donor_id", "sample_id", "dataset"]:
    if col in adata.obs.columns:
        batch_col = col
        break
assert batch_col is not None, "No batch column found"

print(f"Data: {adata.shape}")
print(f"Batches: {adata.obs[batch_col].value_counts().to_dict()}")

# Check if PCA exists
assert "X_pca" in adata.obsm, "PCA required. Run QC preprocessing first."

Step 2: Pre-Integration Visualization

# Compute UMAP on uncorrected data
sc.pp.neighbors(adata, n_pcs=30, use_rep="X_pca")
sc.tl.umap(adata)

fig, axes = plt.subplots(1, 2, figsize=(14, 6))
sc.pl.umap(adata, color=batch_col, ax=axes[0], show=False, title="Before Integration (Batch)")
if "leiden" in adata.obs.columns:
    sc.pl.umap(adata, color="leiden", ax=axes[1], show=False, title="Before Integration (Cluster)")
plt.savefig("pre_integration_umap.png", dpi=150, bbox_inches="tight")
plt.close()

Step 3: Method 1 - Harmony

def run_harmony(adata, batch_col):
    """Run Harmony integration."""
    import harmonypy as hm

    adata_harmony = adata.copy()
    sc.external.pp.harmony_integrate(
        adata_harmony,
        batch_col,
        basis="X_pca",
        adjusted_basis="X_pca_harmony",
    )
    sc.pp.neighbors(adata_harmony, n_pcs=30, use_rep="X_pca_harmony")
    sc.tl.umap(adata_harmony)
    return adata_harmony

Step 4: Method 2 - Scanorama

def run_scanorama(adata, batch_col):
    """Run Scanorama integration."""
    import scanpy.external as sce

    adata_scanorama = adata.copy()
    sce.pp.scanorama_integrate(adata_scanorama, batch_col, basis="X_pca", adjusted_basis="X_scanorama")
    sc.pp.neighbors(adata_scanorama, n_pcs=30, use_rep="X_scanorama")
    sc.tl.umap(adata_scanorama)
    return adata_scanorama

Step 5: Method 3 - BBKNN

def run_bbknn(adata, batch_col):
    """Run BBKNN integration."""
    import scanpy.external as sce

    adata_bbknn = adata.copy()
    sce.pp.bbknn(adata_bbknn, batch_key=batch_col, use_rep="X_pca")
    sc.tl.umap(adata_bbknn)
    return adata_bbknn

Step 6: Method 4 - scVI (if available)

def run_scvi(adata, batch_col):
    """Run scVI integration."""
    try:
        import scvi

        adata_scvi = adata.copy()
        # Need raw counts
        if "counts" in adata_scvi.layers:
            adata_scvi.X = adata_scvi.layers["counts"]

        scvi.model.SCVI.setup_anndata(adata_scvi, batch_key=batch_col)
        model = scvi.model.SCVI(adata_scvi)
        train(max_epochs=100, early_stopping=True)
        adata_scvi.obsm["X_scVI"] = model.get_latent_representation()

        sc.pp.neighbors(adata_scvi, n_pcs=30, use_rep="X_scVI")
        sc.tl.umap(adata_scvi)
        return adata_scvi
    except ImportError:
        print("scvi-tools not available")
        return None

Step 7: Quantitative Benchmarking

def compute_integration_metrics(adata, batch_col, label_col=None, embedding="X_pca"):
    """Compute integration quality metrics."""
    metrics = {}

    # Silhouette score on batch (lower = better batch mixing)
    X = adata.obsm[embedding]
    batch_labels = adata.obs[batch_col].values

    if len(set(batch_labels)) > 1:
        metrics["silhouette_batch"] = silhouette_score(X, batch_labels)

    # Silhouette score on cell type (higher = better biological preservation)
    if label_col and label_col in adata.obs.columns:
        labels = adata.obs[label_col].values
        if len(set(labels)) > 1:
            metrics["silhouette_bio"] = silhouette_score(X, labels)

    # Graph connectivity (from scib if available)
    try:
        from scib.metrics import graph_connectivity
        metrics["graph_connectivity"] = graph_connectivity(adata, batch_col)
    except:
        pass

    return metrics

# Run metrics on all methods
methods = {}
methods["uncorrected"] = {"adata": adata, "embedding": "X_pca"}

try:
    adata_harmony = run_harmony(adata, batch_col)
    methods["harmony"] = {"adata": adata_harmony, "embedding": "X_pca_harmony"}
except Exception as e:
    print(f"Harmony failed: {e}")

try:
    adata_scanorama = run_scanorama(adata, batch_col)
    methods["scanorama"] = {"adata": adata_scanorama, "embedding": "X_scanorama"}
except Exception as e:
    print(f"Scanorama failed: {e}")

try:
    adata_bbknn = run_bbknn(adata, batch_col)
    methods["bbknn"] = {"adata": adata_bbknn, "embedding": "X_pca"}
except Exception as e:
    print(f"BBKNN failed: {e}")

# Compute metrics
label_col = None
for col in ["cell_type", "celltype", "leiden", "cluster_annotations"]:
    if col in adata.obs.columns:
        label_col = col
        break

all_metrics = []
for name, info in methods.items():
    m = compute_integration_metrics(info["adata"], batch_col, label_col, info["embedding"])
    m["method"] = name
    all_metrics.append(m)

metrics_df = pd.DataFrame(all_metrics)
print("\nIntegration Benchmarking:")
print(metrics_df.to_string(index=False))
metrics_df.to_csv("integration_metrics.csv", index=False)

Step 8: Select Best Method and Visualize

# Select best method (highest silhouette_bio with low silhouette_batch)
if "silhouette_bio" in metrics_df.columns:
    # Composite score: bio preservation - batch mixing
    metrics_df["composite"] = metrics_df["silhouette_bio"] - metrics_df["silhouette_batch"]
    best_method = metrics_df.loc[metrics_df["composite"].idxmax(), "method"]
else:
    # Just pick lowest batch silhouette
    best_method = metrics_df.loc[metrics_df["silhouette_batch"].idxmin(), "method"]

print(f"\nBest method: {best_method}")
best_adata = methods[best_method]["adata"]

# Visualize all methods
n_methods = len(methods)
fig, axes = plt.subplots(2, n_methods, figsize=(5 * n_methods, 10))
if n_methods == 1:
    axes = axes.reshape(-1, 1)

for i, (name, info) in enumerate(methods.items()):
    ad = info["adata"]
    sc.pl.umap(ad, color=batch_col, ax=axes[0, i], show=False, title=f"{name} (Batch)")
    if label_col:
        sc.pl.umap(ad, color=label_col, ax=axes[1, i], show=False, title=f"{name} (Cell Type)")

plt.tight_layout()
plt.savefig("integration_comparison.png", dpi=150, bbox_inches="tight")
plt.close()

Step 9: Save Best Result

best_adata.write_h5ad("integrated.h5ad")
metrics_df.to_csv("integration_metrics.csv", index=False)
print(f"Saved integrated data: {best_adata.shape}")

Method Selection Guide

Key Parameters

Common Issues

• Over-correction: Biological variation lost; reduce correction strength
• Under-correction: Batch effects remain; increase correction or try different method
• Memory errors: Subsample for benchmarking; use Harmony for large datasets
• Missing cell types: Some batches may lack certain populations; check before integration