Harmony: correcting batch effects across single-cell datasets

The problem Harmony addresses: several batches, one embedding

Single-cell experiments are routinely split into batches. A study may sequence samples on two different runs, combine a public dataset with a freshly generated one, or merge donors processed on different days. The result is that cells cluster by their technical origin rather than by biology, and a downstream differential expression test ends up reporting the batch instead of the cell type.

Harmony is an R package from the Broad Institute for exactly this situation. The project describes itself as integration of large, complex single-cell datasets, and the topics list algorithm, data-integration, r and scrna-seq. The homepage sits on the Broad Institute portals site, and the package is published on CRAN.

The key design decision is what the input is. Harmony does not take raw counts and does not run clustering. It takes a cell embedding you have already computed, in practice the principal component coordinates from your own pipeline, together with a metadata table that says which batch each cell came from. It returns a new embedding in which the batch structure has been removed while the biological structure is preserved as well as the method allows.

That boundary is what makes the tool usable in an existing analysis. If you already have a Seurat object with a pca reduction, you are already most of the way there. Harmony is a correction step inside a workflow rather than a workflow replacement, which is also why it composes with both Seurat and SingleCellExperiment pipelines according to the README.

The data flow: embeddings in, corrected embeddings out

The minimal call takes three things: your embedding matrix, a metadata table, and the name of the covariate that encodes the unwanted variation. In standalone mode that is RunHarmony with the embedding, the metadata and a string.

library(harmony)
my_harmony_embeddings <- RunHarmony(my_pca_embeddings, meta_data, "dataset")

The third argument is the whole point. If the cells were run in two experiments, the covariate is the experiment identifier, and Harmony uses it as the thing to correct away. You pass the metadata, not a numeric vector of labels, so the covariate names are looked up in that table.

When there is more than one thing to correct for, the argument accepts a vector of covariate names, and the README is explicit that Harmony can integrate over multiple covariates:

my_harmony_embeddings <- RunHarmony(
  my_pca_embeddings, meta_data, c("dataset", "donor", "batch_id")
)

In the Seurat path the same function is used, but the first argument is the Seurat object rather than a matrix:

seuratObj <- RunHarmony(seuratObj, "dataset")
seuratObj <- RunUMAP(seuratObj, reduction = "harmony")

Two details in that snippet matter for practice. The PCA cell embeddings have to be computed first through Seurat's own API before RunHarmony is called, and downstream steps should use the harmony embeddings rather than the pca reduction. The second line shows the intended usage: UMAP is computed on the corrected reduction, not on the original one.

The internals are documented elsewhere. The README points to a detailed walkthrough that explores the internal data structures and mechanics of the Harmony algorithm, which is the place to look if you want to know how the correction is computed rather than just how to call it.

Installing the package from CRAN or straight from GitHub

Harmony is on CRAN, and the README gives both installation routes. From CRAN:

install.packages("harmony")

The README says that this takes approximately 15 seconds. That speed is a hint about the dependency footprint: this is a compiled package but it does not drag in a large stack.

For the development version, the README gives a GitHub route through devtools:

devtools::install_github("immunogenomics/harmony", build_vignettes=TRUE)

The README says that this should take less than 5 minutes, and the build_vignettes argument explains why: the vignettes are where the tutorials live, and building them pulls in the documentation toolchain.

System requirements are short and worth reading before you start. The README states that Harmony has been tested on R versions 4.2 and newer, and that it has been tested on Linux, OS X and Windows. It also points to the DESCRIPTION file for details on the required R packages, which is where you should look if a dependency resolves to something unexpected on your platform.

For a first run, the README points to a quick start vignette that demonstrates standalone mode, described as taking about 4 seconds. The expectation the project sets for runtime is stated plainly as well: for a few samples under 100K cells, the integration should finish within seconds.

Standalone mode versus the Seurat and SingleCellExperiment routes

The README describes Harmony as user-friendly and states that it supports some SingleCellExperiment and Seurat R analysis pipelines, and that it can alternatively be used in standalone mode. That gives you two fairly different shapes of use.

In standalone mode you hold a matrix of embeddings and a metadata table, you call RunHarmony, and you get a matrix back. Nothing in your existing objects changes. This is the route to pick if your pipeline is built around a different object model, or if you want to run Harmony on embeddings you computed by some other means, since the function only requires a matrix and a covariate name.

In the Seurat route the correction is applied to the object itself, the result is stored as a reduction named harmony, and the rest of your analysis continues to use that object. The README is specific about the ordering: precompute the PCA embeddings through Seurat's API, then run Harmony, then use the harmony embeddings instead of pca for downstream work, including the UMAP step in the example.

The SingleCellExperiment support is mentioned in one line without a code example, and the vignettes directory exists in the repository. That asymmetry is worth flagging: the README example set is Seurat-heavy, so a SingleCellExperiment user should expect to consult the vignettes rather than copy a snippet.

The same split appears on the documentation side. There is a quick start article for standalone mode, a separate Seurat vignette, and a detailed walkthrough for the algorithm itself, all hosted on the korsunskylab package website linked from the README.

Performance notes: OpenBLAS first, then the ten million cell question

The README keeps performance to two numbered points, and both point at the same file.

The first is that OpenBLAS will make a substantial performance difference, and readers who are not using it are pointed to PERFORMANCE.md. That is a specific, checkable claim about the BLAS implementation behind the linear algebra, and it is actionable: on a machine where R links against a reference BLAS instead of OpenBLAS, the correction step can be substantially slower without anything in the algorithm changing.

The second point is about scale. For very large datasets, described as more than 10M cells, the README points to the OpenMP notes in PERFORMANCE.md and to the project channel. So the package's own documentation splits into a normal regime and a large-data regime, and the second regime is where the notes matter.

For context on the normal regime, the README says that for a few samples under 100K cells the integration should finish within seconds, and describes the standalone quick start as about 4 seconds. Those are the project's own expectations for its own example, not a measurement, and they are useful as a sanity check when you run it for the first time: if a small dataset takes minutes rather than seconds, the OpenBLAS note is the first place to look.

The repository also carries a PERFORMANCE.md at the top level, alongside src/ for the compiled code, which is consistent with a package where the implementation language matters to performance.

Licence, release history and what a CRAN package here signals

The licence is the first thing to check and the least conventional. The project is recorded as using a custom licence that GitHub cannot classify, with the terms in the LICENSE file at the repository root. A commercial user, or anyone redistributing the package, has to read that file rather than assume a standard permissive licence. This article does not interpret its terms.

The release history is thin. There is a single release recorded, version 0.1, tagged on 2020-02-19, which sits oddly against a repository whose last push was on 2026-10-15 and whose README now describes Harmony2 rather than the original version.

The explanation is visible in the repository layout. There is a CRAN-SUBMISSION file and a cran-comments.md file, both of which are artefacts of a CRAN submission, and a CRAN status badge in the README pointing at the CRAN package page. In other words, the versioned releases on GitHub and the version that CRAN serves are not the same numbering scheme, and for most users the CRAN artefact is the one that matters. The README also cites a preprint for Harmony2 on biorxiv, alongside the earlier Harmony v1 manuscript in Nature Methods, so the method has a publication trail and a newer write-up to read.

Other files in the root support the usual R package shape: a DESCRIPTION file that lists the dependencies, a NAMESPACE file, a NEWS file, a man/ directory for documentation, tests/, vignettes/, and CI configuration through .travis.yml and appveyor.yml. A Python user is pointed to the separate harmonypy package by Kamil Slowikowski, which is the same idea in a different language rather than a port maintained in this repository.

Where Harmony is the wrong tool

The clearest boundary is data type. Everything documented here is about single-cell expression data corrected in embedding space. If your problem is integrating imaging data, proteomics, or spatial measurements, nothing in the README suggests Harmony addresses that, and forcing a single-cell correction onto another data type will not be defensible.

The second boundary is the input requirement. Harmony needs a precomputed embedding and a metadata table naming the covariate. If you have not run dimensionality reduction yet, or if you cannot describe which column of the metadata encodes the batch you want removed, the package cannot help. The covariate choice is also scientifically consequential, and the documentation points to the walkthrough rather than giving a decision procedure, so choosing the right covariate is on you.

The third is scale, and it is a documented limit rather than a vague one. Above 10M cells you are outside the regime the README describes, in the territory of the OpenMP notes. Teams planning single-cell studies at that size should read PERFORMANCE.md before committing to this pipeline.

The fourth is the licence and the release scheme. A custom licence needs reading, and a project whose public version history stops at 0.1 while CRAN carries the current build means your version pinning has to come from CRAN rather than from GitHub tags.

Editorial conclusion

Adopt Harmony if you are merging single-cell datasets that were processed in separate batches and you already work in R with Seurat or SingleCellExperiment, since the corrected embeddings replace the pca reduction in an existing pipeline. Do not adopt it expecting it to integrate modalities other than single-cell expression, or as a substitute for reading the algorithm walkthrough. Verify first that the artefact on CRAN matches the Harmony2 code you need, that your R version is 4.2 or newer, and that you have read PERFORMANCE.md, because the project flags OpenBLAS as a substantial performance factor and sends very large datasets above ten million cells to a separate OpenMP note. Note that the repository uses a custom licence, with the terms in the LICENSE file, and the last push was on 2026-10-15.

Frequently asked questions

Official sources

1. immunogenomics/harmony on GitHub
2. Issues
3. Project website
4. README
5. Releases