Repository: theislab/scib (language: python)

Issue:
Update louvain call to a more recent method
Current version raises `BaseException: Could not construct partition: vector` error when calculating metrics that use Louvain clustering algorithm. As discussed in https://github.com/vtraag/louvain-igraph/issues/61 and [scanpy docs](https://scanpy-tutorials.readthedocs.io/en/latest/paga-paul15.html#Clustering-and-PAGA), the code should use `sc.tl.leiden` instead.




Discussion / hints:
Thanks for this PR @johnarevalo. I agree this is a good idea. The tests would also need to be adapted though. Is this the only place we have louvain in the code @mumichae?

Fix patch:
diff --git a/scib/metrics/__init__.py b/scib/metrics/__init__.py
index 13454854..bc4ed44e 100644
--- a/scib/metrics/__init__.py
+++ b/scib/metrics/__init__.py
@@ -1,6 +1,8 @@
+# fmt: off
 from .ari import ari
 from .cell_cycle import cell_cycle
-from .clustering import opt_louvain
+from .clustering import (cluster_optimal_resolution, get_resolutions,
+                         opt_louvain)
 from .graph_connectivity import graph_connectivity
 from .highly_variable_genes import hvg_overlap
 from .isolated_labels import isolated_labels
diff --git a/scib/metrics/clustering.py b/scib/metrics/clustering.py
index cef1dade..18200fa3 100644
--- a/scib/metrics/clustering.py
+++ b/scib/metrics/clustering.py
@@ -1,11 +1,127 @@
 import matplotlib.pyplot as plt
+import numpy as np
 import pandas as pd
 import scanpy as sc
 import seaborn as sns
+from deprecated import deprecated
 
 from .nmi import nmi
 
 
+def get_resolutions(n=20, min=0.1, max=2):
+    min = np.max([1, int(min * 10)])
+    max = np.max([min, max * 10])
+    frac = n / 10
+    return [frac * x / n for x in range(min, max + 1)]
+
+
+def cluster_optimal_resolution(
+    adata,
+    label_key,
+    cluster_key,
+    cluster_function=None,
+    metric=None,
+    resolutions=None,
+    use_rep=None,
+    force=True,
+    verbose=True,
+    return_all=False,
+    **kwargs,
+):
+    """Optimised clustering
+
+    Leiden, louvain or any custom clustering algorithm with resolution optimised against a metric
+
+    :param adata: anndata object
+    :param label_key: name of column in adata.obs containing biological labels to be
+        optimised against
+    :param cluster_key: name of column to be added to adata.obs during clustering.
+        Will be overwritten if exists and ``force=True``
+    :param cluster_function:
+    :param metric: function that computes the cost to be optimised over. Must take as
+        arguments ``(adata, group1, group2, **kwargs)`` and returns a number for maximising
+        Default is :func:`~scib.metrics.nmi()`
+    :param resolutions: list of resolutions to be optimised over. If ``resolutions=None``,
+        default resolutions of 10 values ranging between 0.1 and 2 will be used
+    :param use_rep: key of embedding to use only if ``adata.uns['neighbors']`` is not
+        defined, otherwise will be ignored
+    :param force: whether to overwrite the cluster assignments in the ``.obs[cluster_key]``
+    :param verbose: whether to print out intermediate results
+    :param return_all: whether to results for all resolutions
+    :param kwargs: arguments to pass to clustering function
+    :returns:
+        Only if ``return_all=True``, return tuple of ``(res_max, score_max, score_all)``
+        ``res_max``: resolution of maximum score;
+        ``score_max``: maximum score;
+        ``score_all``: ``pd.DataFrame`` containing all scores at resolutions. Can be used to plot the score profile.
+    """
+    if cluster_key in adata.obs.columns:
+        if force:
+            print(
+                f"WARNING: cluster key {cluster_key} already exists in adata.obs and will be overwritten because "
+                "force=True "
+            )
+        else:
+            raise ValueError(
+                f"cluster key {cluster_key} already exists in adata, please remove the key or choose a different "
+                "name. If you want to force overwriting the key, specify `force=True` "
+            )
+
+    if cluster_function is None:
+        cluster_function = sc.tl.leiden
+
+    if metric is None:
+        metric = nmi
+
+    if resolutions is None:
+        resolutions = get_resolutions(n=10, max=2)
+
+    score_max = 0
+    res_max = resolutions[0]
+    clustering = None
+    score_all = []
+
+    if use_rep is None:
+        try:
+            adata.uns["neighbors"]
+        except KeyError:
+            raise RuntimeError(
+                "Neighbours must be computed when setting use_rep to None"
+            )
+    else:
+        print(f"Compute neighbors on rep {use_rep}")
+        sc.pp.neighbors(adata, use_rep=use_rep)
+
+    for res in resolutions:
+        cluster_function(adata, resolution=res, key_added=cluster_key, **kwargs)
+        score = metric(adata, label_key, cluster_key)
+        if verbose:
+            print(f"resolution: {res}, {metric.__name__}: {score}")
+        score_all.append(score)
+
+        # optimise score
+        if score_max < score:
+            score_max = score
+            res_max = res
+            clustering = adata.obs[cluster_key]
+        del adata.obs[cluster_key]
+
+    if verbose:
+        print(f"optimised clustering against {label_key}")
+        print(f"optimal cluster resolution: {res_max}")
+        print(f"optimal score: {score_max}")
+
+    score_all = pd.DataFrame(
+        zip(resolutions, score_all), columns=("resolution", "score")
+    )
+
+    adata.obs[cluster_key] = clustering
+
+    if return_all:
+        return res_max, score_max, score_all
+
+
+@deprecated
 def opt_louvain(
     adata,
     label_key,


Test patch:
diff --git a/tests/preprocessing/test_clustering.py b/tests/preprocessing/test_clustering.py
index e3a84fa5..004a4743 100644
--- a/tests/preprocessing/test_clustering.py
+++ b/tests/preprocessing/test_clustering.py
@@ -1,4 +1,5 @@
 import pandas as pd
+import scanpy as sc
 
 import scib
 from tests.common import LOGGER, assert_near_exact
@@ -18,3 +19,34 @@ def test_opt_louvain(adata_neighbors):
     LOGGER.info(f"max resolution: {res_max}, max score: {score_max}")
     assert res_max == 0.7
     assert_near_exact(score_max, 0.7432787576640969, diff=1e-3)
+
+
+def test_cluster_optimal_resolution_louvain(adata_neighbors):
+    res_max, score_max, score_all = scib.cl.cluster_optimal_resolution(
+        adata_neighbors,
+        label_key="celltype",
+        cluster_key="cluster",
+        cluster_function=sc.tl.louvain,
+        resolutions=scib.cl.get_resolutions(n=20, min=0.1, max=2),
+        return_all=True,
+    )
+    assert isinstance(score_all, pd.DataFrame)
+
+    LOGGER.info(f"max resolution: {res_max}, max score: {score_max}")
+    assert res_max == 0.7
+    assert_near_exact(score_max, 0.7432787576640969, diff=1e-3)
+
+
+def test_cluster_optimal_resolution_leiden(adata_neighbors):
+    res_max, score_max, score_all = scib.cl.cluster_optimal_resolution(
+        adata_neighbors,
+        label_key="celltype",
+        cluster_key="cluster",
+        cluster_function=sc.tl.leiden,
+        return_all=True,
+    )
+    assert isinstance(score_all, pd.DataFrame)
+
+    LOGGER.info(f"max resolution: {res_max}, max score: {score_max}")
+    assert res_max == 0.5
+    assert_near_exact(score_max, 0.7424614219722735, diff=1e-3)