add harmony (#317)

* add harmony

* fix formatting

* update exception

* add test for harmony

* install harmony

Co-authored-by: Michaela Mueller <[email protected]>
Co-authored-by: Michaela Mueller <[email protected]>

.github/workflows/test.yml

@@ -93,7 +93,7 @@ jobs:
         -   name: Install dependencies
             run: |
                 python -m pip install --upgrade pip
-                pip install '.[test,scanorama,scvi]'
+                pip install '.[test,scanorama,scvi,harmony]'
 
         -   name: Test with pytest
             run: |

scib/integration.py

@@ -10,6 +10,29 @@
 from .exceptions import OptionalDependencyNotInstalled
 
 
+def harmony(adata, batch, hvg=None, **kwargs):
+    """Harmony wrapper function
+
+    Based on `harmony-pytorch <https://github.com/lilab-bcb/harmony-pytorch>`_ version 0.1.7
+
+    :param adata: preprocessed ``anndata`` object
+    :param batch: batch key in ``adata.obs``
+    :param hvg: list of highly variables to subset to. If ``None``, the full dataset will be used
+    :return: ``anndata`` object containing the corrected feature matrix as well as an embedding representation of the
+        corrected data
+    """
+    try:
+        from harmony import harmonize
+    except ModuleNotFoundError as e:
+        raise OptionalDependencyNotInstalled(e)
+
+    utils.check_sanity(adata, batch, hvg)
+    sc.tl.pca(adata)
+    adata.obsm["X_emb"] = harmonize(adata.obsm["X_pca"], adata.obs, batch_key=batch)
+
+    return adata
+
+
 def scanorama(adata, batch, hvg=None, **kwargs):
     """Scanorama wrapper function
 

setup.cfg

@@ -83,6 +83,7 @@ scvi = scvi-tools  >=0.16.1
 trvae = trvae  ==1.1.2
 trvaep = trvaep  ==0.1.0
 desc = desc  ==2.0.3
+harmony = harmony-pytorch
 
 [tool.isort]
 include_trailing_comma = true

tests/integration/test_harmony.py

@@ -0,0 +1,22 @@
+import scib
+from tests.common import LOGGER, assert_near_exact
+
+
+def test_harmony(adata_paul15_template):
+    adata = scib.ig.harmony(adata_paul15_template, batch="batch")
+
+    scib.pp.reduce_data(
+        adata, n_top_genes=200, neighbors=True, use_rep="X_emb", pca=True, umap=False
+    )
+
+    # check NMI after clustering
+    res_max, score_max, _ = scib.cl.opt_louvain(
+        adata,
+        label_key="celltype",
+        cluster_key="cluster",
+        plot=False,
+        inplace=True,
+    )
+    LOGGER.info(f"max resolution: {res_max}, max score: {score_max}")
+
+    assert_near_exact(score_max, 0.4954028837739288, 1e-2)