Bio Applied Sc Integration

Single-Cell Batch Correction and Dataset Integration

When to Use

• Combining scRNA-seq samples processed on different days, 10x chips, operators, or sequencing runs into one analysis
• A UMAP colored by batch shows batches separating into distinct blobs instead of mixing within cell types
• Deciding whether an observed separation is a real biological effect (disease vs. healthy) or a technical batch effect
• Transferring cell-type labels from a labeled reference atlas onto a new unlabeled query dataset
• Benchmarking multiple integration methods (Harmony/scVI/BBKNN) before picking one for a pipeline

Version Compatibility

scanpy ≥1.10, anndata ≥0.10, harmonypy ≥0.0.10, scvi-tools ≥1.1, scikit-learn ≥1.3, Python ≥3.10.

Prerequisites

• pip install scanpy harmonypy scvi-tools scikit-learn
• An AnnData object with raw counts in adata.layers['counts'], normalized/log data in adata.X, and a batch column in adata.obs
• Prior skill: bio-single-cell-preprocessing (QC, normalization, PCA already computed)

Diagnosing Batch Effects

Color the UMAP by batch first — a cluster made of cells from only one batch, or the same cell type forming two separate blobs, is the classic red flag. Quantify it with LISI (Local Inverse Simpson's Index): for each cell, look at its k nearest neighbors and compute the inverse Simpson diversity of their batch labels. High LISI (near the number of batches) = well mixed; low LISI (near 1) = segregated.

Goal: get an objective, per-cell mixing score before and after correction. Approach: build a k-NN graph in the embedding and compute inverse-Simpson diversity of a label column within each cell's neighborhood.

import numpy as np
from sklearn.neighbors import NearestNeighbors

def compute_lisi(X, labels, n_neighbors=30):
    """Simplified LISI: per-cell inverse Simpson diversity of `labels` among its k-NN.

    X: (n_cells, n_dims) embedding (e.g. adata.obsm['X_pca'])
    labels: array-like of length n_cells (batch IDs or cell types)
    Returns: array of per-cell LISI scores. High = well mixed w.r.t. labels.
    """
    labels = np.asarray(labels)
    nbrs = NearestNeighbors(n_neighbors=n_neighbors).fit(X)
    _, indices = nbrs.kneighbors(X)
    unique_labels = np.unique(labels)
    scores = np.empty(X.shape[0])
    for i, neighbor_idx in enumerate(indices):
        neighbor_labels = labels[neighbor_idx]
        counts = np.array([(neighbor_labels == lab).sum() for lab in unique_labels])
        p = counts / counts.sum()
        simpson = np.sum(p ** 2)
        scores[i] = 1.0 / simpson if simpson > 0 else 1.0
    return scores

# lisi_batch near n_batches = good mixing; lisi_celltype should stay high (structure preserved)
# lisi_batch = compute_lisi(adata.obsm['X_pca'][:, :20], adata.obs['batch'], n_neighbors=30)

Harmony Integration

Harmony (Korsunsky et al. 2019, Nature Methods) iteratively soft-clusters cells in PCA space and shifts each batch's cluster centroid onto the shared centroid, without touching the count matrix. It is fast (seconds–minutes, RAM only) and is the standard first choice. Recommended theta (diversity penalty) range: 1-3, default 2 — too high forces mixing even when batches truly differ biologically (e.g. disease vs. healthy).

Goal: produce a batch-corrected PCA embedding for clustering/UMAP. Approach: run PCA with scanpy, hand the embedding to harmonypy.run_harmony, store the corrected coordinates back in obsm.

import scanpy as sc
import harmonypy

def run_harmony_integration(adata, batch_key="batch", n_pcs=20):
    """Batch-correct an AnnData's PCA embedding with Harmony.

    Requires adata.X to be normalized/log-transformed and HVG-selected already.
    Returns adata with adata.obsm['X_pca_harmony'] added.
    """
    sc.pp.pca(adata, n_comps=n_pcs)
    ho = harmonypy.run_harmony(
        adata.obsm["X_pca"], adata.obs, [batch_key], theta=2.0, max_iter_harmony=10
    )
    adata.obsm["X_pca_harmony"] = ho.Z_corr.T
    return adata

# adata = run_harmony_integration(adata, batch_key="batch")
# sc.pp.neighbors(adata, use_rep="X_pca_harmony")
# sc.tl.umap(adata)

Harmony fails when: batches have no overlapping cell types (compositional imbalance — nothing to align), the batch effect is larger than the biological signal (needs a deeper model like scVI), or theta is too high and merges genuinely distinct populations.

scVI: Deep Generative Integration

scVI (Lopez et al. 2018, Nature Methods) is a variational autoencoder trained directly on raw counts using a negative-binomial likelihood, with batch fed as a covariate to both encoder and decoder. Unlike Harmony it works on counts (no log-normalization needed), gives probabilistic latent estimates, and supports batch-aware differential expression natively. scANVI extends it with semi-supervised cell-type labels for label transfer.

Goal: learn a batch-corrected latent space directly from raw counts, usable for both embedding and DE. Approach: register the AnnData with setup_anndata, train an SCVI model, pull out the latent representation.

import scvi

def run_scvi_integration(adata, batch_key="batch", counts_layer="counts", n_latent=10, max_epochs=400):
    """Train scVI and return the trained model; adds adata.obsm['X_scVI'].

    adata.layers[counts_layer] must hold raw integer counts.
    """
    scvi.model.SCVI.setup_anndata(adata, batch_key=batch_key, layer=counts_layer)
    model = scvi.model.SCVI(adata, n_layers=2, n_latent=n_latent, gene_likelihood="nb")
    model.train(max_epochs=max_epochs, early_stopping=True)
    adata.obsm["X_scVI"] = model.get_latent_representation()
    return model

# model = run_scvi_integration(adata, batch_key="batch")
# de_df = model.differential_expression(adata, groupby="cell_type",
#                                        group1="T_cell", group2="B_cell", delta=0.25)

BBKNN (Batch Balanced k-NN) is a faster, lighter alternative: for each cell it finds k neighbors from each batch separately, then merges these into one graph before UMAP — guaranteeing cross-batch edges without altering the PCA space.

Label Transfer

Transfer cell-type labels from a labeled reference (e.g. Azimuth PBMC atlas) onto an unlabeled query without re-annotating from scratch. Simplest approach: KNN in a shared embedding, confidence = fraction of neighbors sharing the predicted label. Score > 0.7 = reliable; 0.4-0.7 = inspect markers manually; < 0.4 = possibly a novel/transitional state.

Goal: assign cell-type labels to query cells and flag low-confidence calls. Approach: fit a distance-weighted KNN classifier on the reference embedding, predict and threshold on predict_proba.

from sklearn.neighbors import KNeighborsClassifier

def transfer_labels_knn(X_ref, ref_labels, X_query, n_neighbors=15, confidence_threshold=0.7):
    """KNN-based label transfer from a labeled reference to an unlabeled query.

    X_ref, X_query: embeddings in the SAME coordinate space (e.g. joint PCA/scVI latent).
    Returns (predicted_labels, confidence, high_confidence_mask).
    """
    clf = KNeighborsClassifier(n_neighbors=n_neighbors, weights="distance")
    clf.fit(X_ref, ref_labels)
    predicted = clf.predict(X_query)
    confidence = clf.predict_proba(X_query).max(axis=1)
    high_conf = confidence > confidence_threshold
    return predicted, confidence, high_conf

For anchor-based transfer across platforms use Seurat CCA (R); for probabilistic transfer with uncertainty use scANVI, which trains a semi-supervised VAE on scVI's latent space using the reference labels as weak supervision.

Choosing and Benchmarking a Method

Per the scib benchmark (Luecken et al. 2022), no single method wins everywhere — always validate that LISI-batch improves while LISI-celltype does not collapse.

Pitfalls

• Over-correction: aggressive Harmony theta or excessive scVI training can merge biologically distinct populations (e.g. disease vs. healthy immune composition) — always check LISI-celltype didn't drop
• DE on corrected embeddings: never run differential expression on Harmony/scVI-corrected coordinates directly for hypothesis testing; use raw counts (or scVI's own differential_expression, which accounts for batch properly)
• Compositional imbalance: if one batch is all T cells and another all B cells, no method can correctly align them — inspect cell-type composition per batch before integrating
• Correcting away real biology: only integrate to remove technical batch, never to erase condition/treatment effects you are trying to detect
• Multiple testing: apply FDR (Benjamini-Hochberg) when comparing integration metrics or DE genes across many features

See Also

• bio-single-cell-preprocessing — QC, normalization, and PCA required before integration
• bio-single-cell-clustering — Leiden/Louvain clustering on the corrected embedding
• bio-single-cell-cell-annotation — marker-based annotation to validate transferred labels
• bio-single-cell-multimodal-integration — integrating across modalities (RNA + ATAC/protein) rather than batches