harmonypy

harmonypy is a Python package for the Harmony algorithm for integrating multiple high-dimensional datasets. It uses a multithreaded C++ backend that follows the R harmony2 package step-by-step.

This animation shows Harmony aligning three single-cell RNA-seq datasets from different donors. → How to make this animation. Before Harmony, you can clearly distinguish cells from each of the three donors. After Harmony, the cells from different donors are mixed while preserving the overall shape of the data.

Installation

Install from PyPI (pre-built wheels for Linux and macOS):

pip install harmonypy

Building from source

Building from source requires a C++17 compiler and the Python development headers (python3-dev on Debian and Ubuntu, python3-devel on Fedora); harmonypy does not use BLAS or LAPACK. CMake is installed during the build if it is missing, and the Armadillo headers are downloaded (which needs network access) if Armadillo is not installed:

pip install .

Quick Start

import harmonypy as hm
import pandas as pd

# Load the principal components and metadata
pcs = pd.read_csv("data/pbmc_3500_pcs.tsv.gz", sep="\t")
meta = pd.read_csv("data/pbmc_3500_meta.tsv.gz", sep="\t")

# Run Harmony to correct for batch effects (donor)
harmony_out = hm.run_harmony(pcs, meta, "donor")

# Save corrected PCs (same shape as input)
result = pd.DataFrame(harmony_out.Z_corr, columns=pcs.columns)
result.to_csv("pbmc_3500_pcs_harmony.tsv", sep="\t", index=False)

Usage with Scanpy

import scanpy as sc
import harmonypy as hm

# Load and preprocess your data
adata = sc.read_h5ad("my_data.h5ad")
sc.pp.pca(adata)

# Get PCs from the AnnData object
pcs = adata.obsm['X_pca']
print(pcs.shape)  # (n_cells, n_pcs)

# Run Harmony on the PCA embedding
harmony_out = hm.run_harmony(pcs, adata.obs, "batch")

# Store corrected PCs back in the AnnData object
adata.obsm['X_pca_harmony'] = harmony_out.Z_corr

# Use harmonized PCs for downstream analysis
sc.pp.neighbors(adata, use_rep='X_pca_harmony')
sc.tl.umap(adata)
sc.tl.leiden(adata)

Parameters

run_harmony accepts the same parameters as the R package:

Parameter	Default	Description
theta	2	Diversity penalty per batch variable
sigma	0.1	Kernel bandwidth for soft clustering
nclust	min(N/30, 100)	Number of clusters
max_iter_harmony	10	Maximum Harmony iterations
max_iter_kmeans	4	K-means iterations per Harmony round
epsilon_harmony	1e-2	Convergence threshold
ncores	0	Threads (0 = one per physical core available to the process)
lamb	None	Ridge penalty (None = auto-estimate)

harmonypy runs on its own pool of threads (it uses no BLAS, LAPACK or OpenMP). ncores sets the number of threads; the default (0) uses one per physical core available to the process. A second thread per core (hyperthreading) made Harmony about 5% faster on a 6-core laptop and no faster on a 64-core server, where it used 70% more CPU time. On Linux the default respects CPU affinity (e.g. Slurm or taskset) but not container CPU quotas, so set ncores explicitly there. Results are identical for any ncores.

Performance

Run time with default settings on an Apple M1 Ultra (20 cores):

  Dataset                                   2.0.2    this version
  --------------------------------------- -------- ---------------
  Small (3.5k cells, 3 donors)               0.20s           0.06s
  Medium (69k cells, 11 batches)             4.92s           0.24s
  Large (858k cells, 120 batches)           68.34s           1.92s
  Large, by batch and sample (870 levels)  160.15s           3.05s

On one thread (ncores=1) the large dataset takes 18.0 s, 3.8x faster than 2.0.2; using all 20 cores brings it to 1.92 s. On a 6-core Linux laptop (AMD Ryzen 5 5560U), compared with the 2.0.2 wheel:

  Dataset                                   2.0.2    this version
  --------------------------------------- -------- ---------------
  Small (3.5k cells, 3 donors)               0.33s           0.06s
  Medium (69k cells, 11 batches)             6.86s           0.54s
  Large (858k cells, 120 batches)          113.78s           5.18s
  Large, by batch and sample (870 levels)  244.32s           9.15s

scripts/compare_outputs.py compares the results and run times of two builds.

On Tahoe-100M (95.6 million cells, 1,344 samples, 50 PCs), this version corrects all cells in 4.9 minutes on a 64-core server, where 2.0.2 took 4.5 hours; with the same number of iterations, it is 16-28x faster on 1-16 million cells and uses 26-28% less memory. See notebooks/tahoe-100m-benchmark/.

Citation

If you use Harmony in your work, please cite the original paper:

Korsunsky, I., Millard, N., Fan, J. et al. Fast, sensitive and accurate integration of single-cell data with Harmony. Nat Methods 16, 1289–1296 (2019). https://doi.org/10.1038/s41592-019-0619-0

The Supplementary Information PDF provides detailed mathematical descriptions and implementation notes.

To learn more about Harmony 2, please see the preprint here:

Patikas, Nikolaos, Hongcheng Yao, Roopa Madhu, Soumya Raychaudhuri, Martin Hemberg, and Ilya Korsunsky. 2026. Integration of Large, Complex Single-Cell Datasets with Harmony2. bioRxiv. https://doi.org/10.64898/2026.03.16.711825