/harmony-batch-correction

name: "harmony-batch-correction"
description: "Harmony batch correction for scRNA-seq and other omics. Removes batch effects from PCA embeddings while preserving biology. Run after PCA, before UMAP. Scales to millions of cells. Python (harmonypy, scanpy) and R (Seurat)."
license: "MIT"

Harmony Batch Correction

Overview

Harmony is a fast, scalable algorithm for batch integration in single-cell data. It takes a PCA embedding (cells × PCs) as input and returns a corrected embedding from which batch effects have been regressed out via iterative soft-clustering and per-cluster linear regression. The corrected embedding is then used to compute neighbors, UMAP, and downstream clustering — the raw count matrix is never modified. Harmony works for single-cell RNA-seq, ATAC-seq, and other omics modalities where a PCA-like embedding is available.

When to Use

• Integrating scRNA-seq datasets from different samples, donors, sequencing runs, or experimental batches that should contain the same cell types
• Removing technical variation (library preparation protocol, 10x chemistry version, sequencing depth, sequencing platform) while preserving biological differences between cell types and conditions
• Performing fast, scalable batch correction on datasets with millions of cells where deep generative model training would be prohibitively slow
• Correcting for multiple confounding variables simultaneously (batch, donor, sequencing platform, tissue processing protocol)
• Preparing a corrected embedding for UMAP visualization, Leiden clustering, or label transfer without modifying the gene expression count matrix
• Use **omics-plotting** SKILL for generic result-table figures; batch-evaluation UMAPs use scanpy `sc.pl.*`
• Use **scVI/scvi-tools** instead when you need probabilistic batch correction with a variational autoencoder (deep learning), differential expression with uncertainty estimates, or multi-modal integration (RNA + protein)
• Use **BBKNN** instead when you want graph-based integration that avoids constructing a corrected embedding altogether and directly builds a cross-batch nearest-neighbor graph
• Use **Seurat Integration / CCA** (R) instead when you are already in a Seurat workflow and prefer anchor-based integration methods

Prerequisites

• **Python packages**: `harmonypy`, `scanpy>=1.10`, `leidenalg`, `igraph`, `anndata`, `pandas`, `matplotlib`
• **R packages** (optional): `harmony`, `Seurat>=4.0` (for R workflow in Step 7)
• **Data requirements**: AnnData object with raw counts, batch/sample metadata in `adata.obs`, and PCA embedding already computed (`adata.obsm["X_pca"]`)
• **Environment**: Python 3.9+; 8 GB RAM sufficient for up to ~500k cells; linear memory scaling

pip install harmonypy "scanpy[leiden]" anndata pandas matplotlib

# R installation (for Step 7 / Seurat integration)
install.packages("harmony")
# If using Seurat:
install.packages("Seurat")

Pre-flight Interview

Settle these with the user before writing any analysis code.

decisions:
  - id: D1
    param: batchKey
    kind: required
    source: data
    ask: "Which variable separates samples that were processed apart - donor, run, chemistry, or a combination?"
    default: null

  - id: D2
    param: inputEmbedding
    kind: derived
    source: upstream
    ask: "Which reduced space should be corrected?"
    default: "the PCA embedding computed upstream"

  - id: D3
    param: diversityPenalty
    kind: required
    source: user
    depends_on: [D1]
    ask: "How hard should batches be pushed together? Too hard erases real differences between samples."
    default: "2.0; raise toward 4-6 when batch sizes are very uneven"

  - id: D4
    param: numClusters
    kind: optional_conditional
    source: data
    ask: "Does the dataset contain rare populations that need more soft clusters to survive correction?"
    default: "automatic, min(100, cells/30)"

  - id: D5
    param: maxIterations
    kind: optional_conditional
    source: data
    ask: "Did the correction converge, or does it need more iterations?"
    default: 10

  - id: D6
    param: smallClusterProtection
    kind: optional_conditional
    source: user
    ask: "Are small populations being over-split during correction?"
    default: "off"

  - id: D7
    param: kernelBandwidth
    kind: optional
    source: user
    ask: "Should Harmony use its standard soft-clustering bandwidth, or be tuned for unusually heterogeneous batches?"
    default: 0.1

D3 is where batch correction goes wrong in the direction nobody checks. Under-correction is visible - batches stay separated in the embedding. Over-correction is not: samples mix beautifully because a real treatment effect was removed along with the batch effect. Ask what the batch variable is confounded with before raising it.

Quick Start

Minimal pipeline — load preprocessed data, run Harmony via scanpy, produce UMAP:

import scanpy as sc

# Load preprocessed AnnData (counts normalized, HVGs selected, PCA run)
adata = sc.read_h5ad("preprocessed.h5ad")   # must have adata.obs["batch"]

# Run Harmony batch correction (corrects adata.obsm["X_pca"] → "X_pca_harmony")
sc.external.pp.harmony_integrate(adata, key="batch")

# Build neighborhood graph on corrected embedding, then UMAP
sc.pp.neighbors(adata, use_rep="X_pca_harmony")
sc.tl.umap(adata)
sc.tl.leiden(adata, resolution=0.5)

# Visualize
sc.pl.umap(adata, color=["batch", "leiden"], ncols=2)
print(f"Clusters: {adata.obs['leiden'].nunique()} | Cells: {adata.n_obs}")

Workflow

Step 1: Quality Control and Preprocessing

Filter, normalize, select highly variable genes (HVGs), and run PCA. Harmony is applied to the PCA embedding produced here.

import scanpy as sc
import numpy as np

sc.settings.set_figure_params(dpi=80, facecolor="white")

# Load multi-batch data — adata.obs must contain a batch column
adata = sc.read_h5ad("multi_batch_counts.h5ad")
# Alternatively, concate