henhao
Sign in
BioinformaticsNature Methods
Peer-reviewed journal article

Fast, sensitive, and accurate integration of single cell data with Harmony

Ilya KorsunskyiD, Nghia MillardiD, Jean FaniD, Kamil SlowikowskiiD, Fan ZhangiD, Kevin WeiiD, Yuriy BaglaenkoiD, Michael B. BrenneriD, Po‐Ru LohiD, Soumya RaychaudhuriiD

Abstract

The rapidly emerging diversity of single cell RNAseq datasets allows us to characterize the transcriptional behavior of cell types across a wide variety of biological and clinical conditions. With this comprehensive breadth comes a major analytical challenge. The same cell type across tissues, from different donors, or in different disease states, may appear to express different genes. A joint analysis of multiple datasets requires the integration of cells across diverse conditions. This is particularly challenging when datasets are assayed with different technologies in which real biological differences are interspersed with technical differences. We present Harmony, an algorithm that projects cells into a shared embedding in which cells group by cell type rather than dataset-specific conditions. Unlike available single-cell integration methods, Harmony can simultaneously account for multiple experimental and biological factors. We develop objective metrics to evaluate the quality of data integration. In four separate analyses, we demonstrate the superior performance of Harmony to four single-cell-specific integration algorithms. Moreover, we show that Harmony requires dramatically fewer computational resources. It is the only available algorithm that makes the integration of[~] 106 cells feasible on a personal computer. We demonstrate that Harmony identifies both broad populations and fine-grained subpopulations of PBMCs from datasets with large experimental differences. In a meta-analysis of 14,746 cells from 5 studies of human pancreatic islet cells, Harmony accounts for variation among technologies and donors to successfully align several rare subpopulations. In the resulting integrated embedding, we identify a previously unidentified population of potentially dysfunctional alpha islet cells, enriched for genes active in the Endoplasmic Reticulum (ER) stress response. The abundance of these alpha cells correlates across donors with the proportion of dysfunctional beta cells also enriched in ER stress response genes. Harmony is a fast and flexible general purpose integration algorithm that enables the identification of shared fine-grained subpopulations across a variety of experimental and biological conditions.

PDFDOICode
Full-text digest ↓

02 /Specification & Classification

Click any tag to browse matching papers
Citations
10,490
FWCI
131.88
Field percentile
top 1%
OpenAlex Topics

In-depth summary

✦ Multi-agent verified full-text synthesis

Claude Opus 5.5 · max effort · 8 Nov 2026 · from the full text of version 2 (5 Nov 2018) · CC BY 4.0

Harmony alternates diversity-penalized soft clustering with cluster-specific linear correction of PCA embeddings, integrating single-cell datasets over multiple covariates and handling 500,000 cells in 7.2GB of memory.

Joint scRNA-seq analysis across donors, tissues, and platforms is confounded by technical variation; prior integrators (Seurat MultiCCA, MNN Correct, Scanorama, BBKNN) integrate over only one variable, and several run out of memory beyond 125,000 cells. Harmony offers cluster-aware linear correction without quadratic cell-to-cell comparisons, models several categorical covariates jointly, and adds local inverse Simpson's index (LISI) metrics for mixing and accuracy.

Key results

  1. ShownHarmony mixed datasets more completely (higher integration LISI, iLISI) than PCA or four competing integrators, both in a controlled Jurkat/HEK293T design where cell-type LISI (cLISI) stayed near 1 and across tissues and donors in the Human Cell Atlas (HCA) benchmark.

    Cell lines: iLISI median 1.59 (95% [1.27, 1.97]) vs 1.01 for PCA (expected maxima 1.8 for Jurkat, 1.5 for 293T), cLISI 1.00 (95% [1.00, 1.02]); MultiCCA lowered iLISI to 1.00 and raised the cLISI 95% upper limit to 1.23. HCA: tissue iLISI 1.40 (95% [1.04, 1.97]) vs other methods' medians 1.00 to 1.12; donor iLISI 3.93 (95% [2.46, 4.95]) vs 1.07 to 2.82. · Quantification of Error and Accuracy in Cell…; Harmony Scales to Enable Analysis of Large Data; Footnotes

  2. ShownHarmony's runtime and memory scaled well from 30,000 to 500,000 cells, keeping atlas-scale integration within typical personal-computer memory, whereas Scanorama, MNN Correct, and Seurat MultiCCA exceeded server memory beyond 125,000 cells.

    HCA subsets: Harmony ran 4 minutes (30,000 cells) to 68 minutes (500,000), 30 to 200 times faster than MNN Correct and MultiCCA, using 0.9GB to 7.2GB; at 125,000 cells it used 30 to 50 times less memory than MultiCCA, MNN Correct, Scanorama. Only BBKNN also finished 500,000 cells (44 minutes, 45GB); Harmony, BBKNN, Scanorama runtimes were similar up to 125,000 cells. · Harmony Scales to Enable Analysis of Large Data; Discussion

  3. ShownAcross three human peripheral blood mononuclear cell (PBMC) libraries from 10x Chromium 3pV1, 3pV2, and 5p chemistries, Harmony aligned the datasets and exposed shared fine-grained T and B cell subtypes that remained dataset-specific in other methods' embeddings.

    Dataset iLISI median 1.96 (95% [1.36, 2.56]) after Harmony vs 1.00 (95% [1.00, 1.00]) before; no other method exceeded a median of 1.1; Harmony cLISI 1.00 (95% [1.00, 1.02]). Joint clustering yielded 5 shared T cell subtypes (naive and effector memory CD4, Treg, memory and effector CD8) and 2 shared B cell subtypes (naive, memory). · Deep Integration Enables Identification of Broad…; Figure legends

  4. ShownIn a meta-analysis of human islet cells from 5 technologies and 36 donors, Harmony integrated over donor and technology simultaneously, recovering rare activated and quiescent stellate cells, mast cells, macrophages, and beta cells showing endoplasmic reticulum (ER) stress.

    Across 14,746 cells, technology iLISI median rose from 1.00 (95% [1.00, 1.06]) under PCA to 2.27 (95% [1.31, 3.27]) and donor iLISI from 1.42 (95% [1.00, 5.50]) to 4.71 (95% [1.81, 6.36]). Rare types except activated stellate cells appeared in all 5 technology datasets. No competing method was run, as the authors state none integrates over more than one variable. · Abstract; Simultaneous Integration Across Donors and…; Figure legends

  5. SuggestedClustering the integrated islet embedding uncovered a previously undescribed alpha-cell subset with ER-stress gene enrichment and lower expression of alpha-cell function genes, whose donor-level abundance correlated with that of ER-stressed beta cells.

    Alpha cluster over-expressed DDIT3, ATF3, ATF4, HSPA5 and had lower GCG, ISL1, ARX, MAFB (limma with technology and library-complexity covariates; Gene Ontology enrichment via enrichR); beta ER-stress cells had lower PDX1, MAFA, INSM1, NEUROD1. Proportions of the two populations correlated across donors (Spearman r=0.46, p=0.004). · Simultaneous Integration Across Donors and…; Discussion; Identification of alpha and beta ER stress…

Methods, data and reproducibility

Starting from a PCA embedding, Harmony iterates cosine soft k-means clustering with a Kullback-Leibler penalty against batch-dominated clusters, then a Gaussian mixture-model linear correction giving each cell a cluster-weighted batch offset, until convergence; iLISI and cLISI score dataset mixing and cell-type separation.

  • Cell-line scRNA-seq: pure Jurkat, pure 293T, and a 50:50 Jurkat/293T mixture, used as a controlled integration benchmark · human (Jurkat and HEK293T cell lines) · n = 3255 Jurkat, 2859 293T (pure); 1799 Jurkat, 1565 293T (mixture)
  • Human Cell Atlas (HCA) cord blood and bone marrow scRNA-seq, downsampled to 500,000, 250,000, 125,000, 60,000, and 30,000 cells · human · n = 528,688 cells; 16 donors; 2 tissues
  • Three PBMC datasets from 10x Chromium with 3-prime v1 (3pV1), 3-prime v2 (3pV2), and 5-prime chemistries; 3pV1 requantified from FASTQ · human
  • Pancreatic islet scRNA-seq from five studies on CEL-seq, CEL-seq2, Fluidigm C1, SmartSeq2, and inDrop platforms · human · n = 14,746 cells; 36 donors
Methods & tools
Harmony (R package; maximum diversity clustering plus mixture-model linear correction), LISI, iLISI, cLISI (local inverse Simpson's index; perplexity 30), MNN Correct (scran mnncorrect, version 1.9.4), Seurat MultiCCA (Seurat version 2.3.4), Scanorama (assemble function, knn=30, sigma=1), BBKNN with scanpy UMAP, Cell Ranger v2.1.0 (GRCh38), Truncated SVD/PCA (top 30 components), UMAP (k = 30, correlation distance, min dist = 0.1), Louvain clustering (Seurat BuildSNN, RunModularityClustering), limma (differential expression with covariates), enrichR (Gene Ontology enrichment)
Code
https://github.com/immunogenomics/harmonyhttps://github.com/immunogenomics/LISI
Reagents & environment
Harmony R package; Cell Ranger v2.1.0 (GRCh38); scran 1.9.4 mnncorrect; Seurat 2.3.4 MultiCCA; Scanorama and BBKNN from GitHub; 30 PCs input; Harmony defaults θ = 2, K = 100, σ = 0.1, α = 0.1, block_size 0.05; islets τ = 5, technology θ = 4.
Assay conditions & compute
Linux servers with Intel Xeon E5-2690 v3 CPUs, each job allotted 6 cores and 120GB RAM (Discussion cites 128GB servers); wall time and peak resident memory from /usr/bin/time -v. Harmony: 4 to 68 minutes and 0.9GB to 7.2GB for 30,000 to 500,000 cells.
Reproducibility gaps
Data source links are deferred to supplementary Table S6, with no accessions in the text; versions of Harmony itself, Scanorama, and BBKNN are not stated; centroid initialization is random and no seed is reported.
Limitations (6)
  • The clustering penalty accepts several covariates, but the linear correction (regression) step currently handles only one batch variable; in the islet analysis donor was the primary covariate (θ = 2) and technology secondary (θ = 4). · stated by the authors
  • Only additive batch offsets in the reduced embedding are modeled (scaling effects are not estimated), and the Gaussian, orthogonal-data assumptions of the correction do not hold for the L2-normalized space used in clustering, so correction runs on unnormalized coordinates. · stated by the authors
  • Only categorical covariates are supported; continuous nuisance factors such as read depth or cell quality are not yet handled by Harmony, unlike global regression approaches. · stated by the authors
  • Harmony returns corrected embedding coordinates rather than adjusted gene expression, so differential expression needs batch-aware models; a very large diversity penalty θ could overmix and overcorrect datasets, so users are urged to test several values. · stated by the authors
  • Competitors shared the authors' preprocessing (pooled instead of per-dataset gene scaling, Scanorama assemble on precomputed PCs rather than its correct function), departing from some tools' own workflows, which may understate their integration performance. · noted in this summary
  • The islet analysis has no baseline comparison (not even the stepwise global-regression route the authors describe); the alpha ER-stress subtype and its donor-level correlation (r=0.46) lack experimental validation and could reflect shared tissue-processing stress. · noted in this summary

Who should read:Read this if you integrate scRNA-seq data over several covariates (donor, technology, tissue) at 10^5 to 10^6 cells, or need LISI metrics scoring dataset mixing separately from cell-type preservation.

Reproducibility ledger

Nobody has reported trying to reproduce this paper yet. Did you replicate an experiment or run the code? Your report saves the next lab a week.

Sign in to file a reproducibility report.

Related papers

Similarity Map
SelectedSame fieldCross-field
THISPatikas '26Li '22Hie '19Ming '22DePasquale '19Stuart '19
82% semantic matchBioinformatics
Nikolaos Patikas, Hongcheng Yao, Roopa Madhu+3 more·bioRxiv·Mar 19, 2026· 8 citations