Harmony Single Cell Rna Seq
Harmony is an algorithm designed to integrate multiple single-cell RNA sequencing (scRNA-seq) datasets by removing batch effects while preserving biological variation. This guide explains how Harmony works, when to use it, and how to implement it in your analysis pipeline. It is written for computational biologists, bioinformatics trainees, and researchers who process scRNA-seq data and need a practical, source-bounded framework. If you are new to scRNA-seq integration, start by reviewing foundational resources like the NCBI Bookshelf for background on sequencing technology and the EMBL-EBI Training materials for general bioinformatics workflows.
Harmony operates by iteratively adjusting principal component (PC) embeddings to maximize the overlap in cluster assignments across batches. It does not require explicit annotation of cell types and scales well to millions of cells. The method is implemented in the Harmony R package available through Bioconductor. This guide will help you decide if Harmony fits your integration needs, walk you through a typical workflow, and alert you to common mistakes.
At a Glance
| Feature | Detail |
|---|---|
| Purpose | Correct batch effects in scRNA-seq data while preserving biological heterogeneity. |
| Input | Normalized, log-transformed expression matrix and a PC embedding. |
| Output | Corrected PC embeddings that can be used for clustering, visualization, and downstream analysis. |
| Key Assumption | Batches share at least some cell types or states, batch effects are linear or approximately linear in PC space. |
| Strengths | Fast, memory efficient, no need for reference dataset, handles large datasets. |
| Limitations | May not perform well with extreme batch effects or when batches have very different cell type compositions. |
| Software | Harmony R package (Bioconductor), also available in Python via harmony‑py. |
| Learning Resources | Galaxy Training Network offers tutorials on scRNA-seq integration. |
Core Concepts of Harmony Integration
Harmony belongs to the family of integration algorithms that work on low‑dimensional embeddings. After standard scRNA-seq preprocessing (quality control, normalization, highly variable gene selection, PCA), you obtain a matrix of cells by principal components. Harmony takes this matrix and a batch assignment vector as inputs.
The algorithm uses a maximum diversity clustering step to assign cells to soft clusters, then computes a correction term that shifts each cell’s PC coordinates to reduce the dispersion of clusters across batches. This correction is applied iteratively until convergence. A key innovation is that Harmony penalizes the variance of cluster centroids within each batch, effectively pushing batches to share similar cluster structures. The corrected embeddings can then be used for UMAP visualization, clustering, or differential expression analysis.
For a deeper understanding of the statistical underpinnings, the Bioconductor documentation includes detailed vignettes and references to the original methods paper. Several recent applications demonstrate Harmony’s utility. For example, a multi‑omics analysis of arginine metabolism in ovarian cancer used Harmony to integrate scRNA-seq data from multiple patient samples and identified prognostic signatures Multi‑omics analysis of arginine metabolism in ovarian cancer. Similarly, an integrated study of microglial subpopulations in Parkinson’s disease leveraged Harmony to combine datasets from different sources, revealing key therapeutic targets Integrated multi‑omics analysis identifies key microglial subpopulations and therapeutic targets in Parkinson’s disease.
Decision Criteria for Using Harmony
Not every integration problem is well suited to Harmony. You should consider the following criteria before adopting it:
- Number of cells: Harmony scales well to hundreds of thousands or millions of cells. It is a good choice for large atlas projects. For very small datasets (fewer than 500 cells per batch), simpler methods like simple linear regression may suffice.
- Batch effect severity: Harmony assumes batch effects are relatively mild and can be captured by linear shifts in PC space. If your data show strong non‑linear artifacts (e.g., one batch has extremely high dropout rates), Harmony may undercorrect. In such cases, consider methods based on mutual nearest neighbors or deep learning.
- Biological heterogeneity: Harmony preserves biology when batches share at least a few cell types. If your batches are completely disjoint (no overlapping cell types), Harmony may incorrectly force them to align. For such scenarios, you might need reference‑based integration.
- Computational resources: Harmony requires only modest memory and runs quickly. However, for datasets with more than 10 million cells, you may need to sample cells or use a more scalable implementation.
A practical workflow for evaluating these criteria is outlined in the NCBI Bookshelf sections on scRNA-seq analysis design. The Galaxy Training Network also provides decision trees for choosing integration methods.
Practical Workflow and Implementation Steps
Below is a typical implementation sequence using the Harmony R package. The steps assume you have already performed standard preprocessing up to PCA.
Step 1: Load and preprocess data
Start with a Seurat or SingleCellExperiment object. Ensure you have a normalized log‑counts matrix and have computed PCA (e.g., with RunPCA in Seurat). The number of PCs should be chosen based on elbow plot or jackstraw method. In Seurat you can run:
library(Seurat)
library(harmony)
data <- NormalizeData(data)
data <- FindVariableFeatures(data)
data <- ScaleData(data)
data <- RunPCA(data)
Step 2: Run Harmony
Supply the data object and the name of the batch column. For Seurat, use:
data <- RunHarmony(data, group.by.vars = "batch")
The function returns the object with a new reduction called harmony. You can adjust parameters like theta (diversity penalty), sigma (soft clustering bandwidth), and max.iter.harmony (convergence iterations). For most datasets the defaults work well.
Step 3: Use corrected embeddings
You can now cluster and visualize using the harmony reduction:
data <- FindNeighbors(data, reduction = "harmony", dims = 1:20)
data <- FindClusters(data, resolution = 0.5)
data <- RunUMAP(data, reduction = "harmony", dims = 1:20)
DimPlot(data, reduction = "umap", group.by = "batch")
Check that cells from different batches mix well in UMAP space.
Step 4: Downstream analysis
Proceed with cell type annotation, differential expression, or trajectory inference using the corrected embeddings. Harmony does not alter the raw expression matrix, so you can always revert to original counts for gene‑level analyses.
Detailed code examples are available in the Bioconductor package vignettes and in public workflows hosted on the Galaxy Training Network. Data for practice can be obtained from the NCBI Sequence Read Archive.
Quality Checks and Validation
After integration, you must verify that batch effects are reduced and biological variation is preserved. Common quality checks include:
- Visual inspection of UMAP or t‑SNE plots colored by batch. Mixing of batches in shared cell types is a good sign.
- Silhouette score on batch labels. Compute the average silhouette width for each cell using its batch identity and the corrected embeddings. A score near zero indicates good mixing.
- Comparison of cluster proportions per batch. For clusters that represent known cell types, the proportion of cells from each batch should be similar (unless the batch truly lacks that cell type).
- Differential expression analysis before and after integration. If a known marker gene shows batch‑correlated expression that disappears after integration, the correction is working.
- Preservation of biological signal. Check that known cell type markers remain differentially expressed between clusters after integration. For guidance on these metrics, refer to the [NC
BI Bookshelf](https://www.ncbi.nlm.nih.gov/books/) chapters on scRNA-seq quality control and the EMBL-EBI Training materials on integration validation.
A practical example comes from a study of gliomas that used Harmony to integrate primary, recurrent, and grade III tumor samples. The authors validated that the integrated data preserved known transcriptional dynamics of GFAP+ cells while removing donor‑specific batch effects GFAP-Dependent Transcriptional Dynamics and Cellular Heterogeneity in Primary, Recurrent, and Grade III Gliomas. Similarly, a single‑cell analysis of periodontitis used Harmony to combine oral mucosa samples and then confirmed neutrophil‑mast cell communication programs were reproducible across donors Single-cell Transcriptomics Inference of Neutrophil-mast Cell Communication Programs in Periodontitis.
Common Mistakes and Pitfalls
Even with a robust algorithm like Harmony, mistakes can derail your analysis.
- Incorrect preprocessing: Harmony works best on properly normalized log‑counts. Using raw counts or incorrect scaling can produce poor results.
- Too many or too few PCs: Using too few PCs can discard biological signal, too many can include noise. Always choose the number of PCs based on data variance (e.g., elbow plot) and not arbitrarily.
- Ignoring batch composition differences: If one batch contains a unique cell type not present in others, Harmony may attempt to erase that difference. Inspect the silhouette score per cluster to detect forced merging.
- Applying Harmony globally without subgrouping: For complex multi‑condition experiments, first integrate within conditions, then use Harmony again to integrate conditions.
- Running Harmony on highly variable genes only: While acceptable, ensure that the gene selection is not biased by batch effects. Consider using a harmonized HVG selection that accounts for batches.
- Not checking convergence: The default number of iterations is usually sufficient, but for very large datasets you may need to increase
max.iter.harmony. Monitor if the objective function stabilizes.
These pitfalls are discussed in depth in the Galaxy Training Network tutorials on scRNA-seq integration. Real‑world applications, such as the integrative analysis of TRIM31+ colorectal tumor cells, explicitly acknowledge the need to test multiple integration parameters to avoid over‑correction Integrative single‑cell analysis reveals TRIM31(+) colorectal tumor cells orchestrating macrophage crosstalk within the cancer‑immunity regulome.
Limits of Interpretation and Uncertainty
No integration method is perfect, and Harmony has well‑defined limitations.
- Biological assumption: Harmony assumes that cell type compositions are not vastly different across batches. If one batch has a cell type that is completely absent in another, Harmony may incorrectly force that cell type to align with a different cluster. This can lead to false biological conclusions. Always examine whether the mixing is biologically plausible.
- Linear correction: The algorithm operates in PC space and applies linear corrections. Strong non‑linear batch effects (e.g., caused by different library preparation protocols with vastly different sensitivity) may not be adequately corrected. In such cases, consider methods like scVI or scANVI.
- Preservation of subtle populations: Very rare cell types (less than 0.1% of cells) can be lost during integration if they are not consistently present across batches. Harmony may assign them to the nearest abundant cluster. Use dedicated cell‑type detection methods before integration.
- Statistical uncertainty: The corrected embeddings are point estimates. Downstream analyses that depend on these embeddings (e.g., pseudotime, differential abundance) should be interpreted with caution. Bootstrapping or repeated runs with different random seeds can help assess stability.
- No expression‑level correction: Harmony does not output corrected expression values. If you need to perform differential expression between conditions across batches, you must apply a separate normalization step (e.g., using ComBat-seq). See the Bioconductor page for guidance on combining methods.
A cautionary example is found in the analysis of a tumor‑specific T‑cell metabolic signature in ovarian cancer, the authors used Harmony but noted that some low‑abundance T‑cell subpopulations required additional validation with flow cytometry Single‑cell analysis identifies a tumor‑specific T‑cell metabolic signature: prognostic model and association with immunosuppressive microenvironment in ovarian cancer. Always cross‑validate integration results with independent experiments or external references.
Frequently Asked Questions
Q1: Can Harmony integrate data from different sequencing technologies (e.g., 10x Genomics vs. Smart-seq2)? Yes, but success depends on the severity of technical differences. Harmony works best when the biological variation outweighs technical variation. If the technologies produce vastly different gene detection rates, pre‑normalization steps (e.g., using SCTransform) may be necessary before applying Harmony.
Q2: How do I choose the number of principal components for Harmony input? Use an elbow plot or jackstraw method on the combined dataset. Typically 20 to 30 PCs capture biological variation without including excessive noise. For very heterogeneous datasets, you may need up to 50 PCs. You can test a range and compare silhouette scores.
Q3: Does Harmony work with multi‑modal data (e.g., CITE‑seq, scATAC‑seq)? Harmony is designed for single‑cell transcriptomics. For multi‑modal data, specialized integration methods like MOFA or Seurat’s weighted nearest neighbor approach are more appropriate. However, you can apply Harmony on the RNA modality alone if your primary focus is gene expression.
Q4: I have only two batches. Should I still use Harmony? Yes, Harmony works with two or more batches. For two batches, simpler methods like canonical correlation analysis (CCA) may also be effective. Harmony may still be preferred if the batches contain many cells and you want a fast, scalable solution.
References and Further Reading
- NCBI Bookshelf , Authoritative textbooks on molecular biology and sequencing technology.
- EMBL-EBI Training , Official online courses covering scRNA-seq data analysis and integration.
- Galaxy Training Network , Interactive tutorials for bioinformatics workflows, including Harmony.
- Bioconductor , Harmony package documentation and vignettes.
- NCBI Sequence Read Archive , Repository for raw sequencing data used in Harmony tutorials.
- Multi‑omics analysis of arginine metabolism in ovarian cancer: A prognostic signature and GTF2F2‑driven stromal remodeling , Example of Harmony integration in cancer research.
- Integrated multi‑omics analysis identifies key microglial subpopulations and therapeutic targets in Parkinson’s disease , Application of Harmony to neurological disease.
- GFAP‑Dependent Transcriptional Dynamics and Cellular Heterogeneity in Primary, Recurrent, and Grade III Gliomas , Validation of Harmony in glioma studies.
- Single‑cell Transcriptomics Inference of Neutrophil‑mast Cell Communication Programs in Periodontitis , Use of Harmony for oral tissue data.
- Integrative single‑cell analysis reveals TRIM31(+) colorectal tumor cells orchestrating macrophage crosstalk within the cancer‑immunity regulome , Workflow demonstrating parameter tuning for Harmony.