(feat): igraph leiden implementation now included as an option in sc.tl.leiden

docs/release-notes/1.10.0.md

@@ -14,6 +14,7 @@
 * {func}`scanpy.pp.pca`, {func}`scanpy.pp.scale`, {func}`scanpy.pl.embedding`, and {func}`scanpy.experimental.pp.normalize_pearson_residuals_pca`
   now support a `mask` parameter {pr}`2272` {smaller}`C Bright, T Marcella, & P Angerer`
 * {func}`scanpy.tl.rank_genes_groups` no longer warns that it's default was changed from t-test_overestim_var to t-test {pr}`2798` {smaller}`L Heumos`
+* {func}`scanpy.tl.leiden` now relies by default on `igraph` fully.  `leidenalg`'s implementation is still available via `flavor` when set to `leidenalg`, although other defaults have now changed.  {pr}`2815` {smaller}`I Gold`
 
 ```{rubric} Docs
 ```

scanpy/_utils/__init__.py

@@ -7,9 +7,11 @@
 
 import importlib.util
 import inspect
+import random
 import sys
 import warnings
 from collections import namedtuple
+from contextlib import contextmanager
 from enum import Enum
 from functools import partial, singledispatch, wraps
 from textwrap import dedent
@@ -20,10 +22,10 @@
 import numpy as np
 from anndata import AnnData
 from anndata import __version__ as anndata_version
-from numpy import random
 from numpy.typing import NDArray
 from packaging import version
 from scipy import sparse
+from sklearn.utils import check_random_state
 
 from .. import logging as logg
 from .._compat import DaskArray
@@ -45,7 +47,38 @@ def __repr__(self) -> str:
 _empty = Empty.token
 
 # e.g. https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html
-AnyRandom = Union[int, random.RandomState, None]  # maybe in the future random.Generator
+# maybe in the future random.Generator
+AnyRandom = Union[int, np.random.RandomState, None]
+
+
+class RNGIgraph:
+    """
+    Random number generator for ipgraph so global seed is not changed.
+    See :func:`igraph.set_random_number_generator` for the requirements.
+    """
+
+    def __init__(self, random_state: int = 0) -> None:
+        self._rng = check_random_state(random_state)
+
+    def __getattr__(self, attr: str):
+        return getattr(self._rng, "normal" if attr == "gauss" else attr)
+
+
+@contextmanager
+def set_igraph_random_state(random_state: int):
+    try:
+        import igraph
+    except ImportError:
+        raise ImportError(
+            "Please install igraph: `conda install -c conda-forge igraph` or `pip3 install igraph`."
+        )
+    rng = RNGIgraph(random_state)
+    try:
+        igraph.set_random_number_generator(rng)
+        yield None
+    finally:
+        igraph.set_random_number_generator(random)
+
 
 EPS = 1e-15
 

scanpy/tests/notebooks/test_pbmc3k.py

@@ -28,11 +28,10 @@
 
 @needs.leidenalg
 def test_pbmc3k(image_comparer):
+    # ensure violin plots and other non-determinstic plots have deterministic behavior
+    np.random.seed(0)
     save_and_compare_images = partial(image_comparer, ROOT, tol=20)
-
-    adata = sc.read(
-        "./data/pbmc3k_raw.h5ad", backup_url="https://falexwolf.de/data/pbmc3k_raw.h5ad"
-    )
+    adata = sc.datasets.pbmc3k()
 
     # Preprocessing
 
@@ -105,13 +104,41 @@ def test_pbmc3k(image_comparer):
 
     # Clustering the graph
 
-    sc.tl.leiden(adata, resolution=0.9)
-    # sc.pl.umap(adata, color=['leiden', 'CST3', 'NKG7'], show=False)
-    # save_and_compare_images('umap_2')
+    sc.tl.leiden(adata, resolution=0.9, random_state=0)
+
+    # sc.pl.umap(adata, color=["leiden", "CST3", "NKG7"], show=False)
+    # save_and_compare_images("umap_2")
     sc.pl.scatter(adata, "CST3", "NKG7", color="leiden", show=False)
     save_and_compare_images("scatter_3")
 
     # Finding marker genes
+    # Due to incosistency with our test runner vs local, these clusters need to
+    # be pre-annotated as the numbers for each cluster are not consistent.
+    marker_genes = [
+        "RP11-18H21.1",
+        "GZMK",
+        "CD79A",
+        "FCGR3A",
+        "GNLY",
+        "S100A8",
+        "FCER1A",
+        "PPBP",
+    ]
+    new_labels = ["0", "1", "2", "3", "4", "5", "6", "7"]
+    data_df = adata[:, marker_genes].to_df()
+    data_df["leiden"] = adata.obs["leiden"]
+    max_idxs = data_df.groupby("leiden", observed=True).mean().idxmax()
+    leiden_relabel = {}
+    for marker_gene, new_label in zip(marker_genes, new_labels):
+        leiden_relabel[max_idxs[marker_gene]] = new_label
+    adata.obs["leiden_old"] = adata.obs["leiden"].copy()
+    adata.rename_categories(
+        "leiden", [leiden_relabel[key] for key in sorted(leiden_relabel.keys())]
+    )
+    # ensure that the column can be sorted for consistent plotting since it is by default unordered
+    adata.obs["leiden"] = adata.obs["leiden"].cat.reorder_categories(
+        list(map(str, range(len(adata.obs["leiden"].cat.categories)))), ordered=True
+    )
 
     sc.tl.rank_genes_groups(adata, "leiden")
     sc.pl.rank_genes_groups(adata, n_genes=20, sharey=False, show=False)
@@ -132,23 +159,22 @@ def test_pbmc3k(image_comparer):
     if adata[adata.obs["leiden"] == "4", "CST3"].X.mean() < 1:
         (  # switch clusters
             adata.obs["leiden"][adata.obs["leiden"] == "4"],
-            adata.obs["leiden"][adata.obs["leiden"] == "5"],
-        ) = ("5", "4")
+            adata.obs["leiden"][adata.obs["leiden"] == "3"],
+        ) = ("3", "4")
     new_cluster_names = [
         "CD4 T cells",
-        "CD14+ Monocytes",
-        "B cells",
         "CD8 T cells",
-        "NK cells",
+        "B cells",
         "FCGR3A+ Monocytes",
+        "NK cells",
+        "CD14+ Monocytes",
         "Dendritic cells",
         "Megakaryocytes",
     ]
     adata.rename_categories("leiden", new_cluster_names)
 
     # sc.pl.umap(adata, color='leiden', legend_loc='on data', title='', frameon=False, show=False)
     # save_and_compare_images('umap_3')
-
     sc.pl.violin(
         adata, ["CST3", "NKG7", "PPBP"], groupby="leiden", rotation=90, show=False
     )

scanpy/tests/test_clustering.py

@@ -1,6 +1,7 @@
 from __future__ import annotations
 
 import pytest
+from sklearn.metrics.cluster import normalized_mutual_info_score
 
 import scanpy as sc
 from scanpy.testing._helpers.data import pbmc68k_reduced
@@ -12,11 +13,80 @@ def adata_neighbors():
     return pbmc68k_reduced()
 
 
+FLAVORS = [
+    pytest.param("igraph", marks=needs.igraph),
+    pytest.param("leidenalg", marks=needs.leidenalg),
+]
+
+
+@pytest.mark.parametrize("flavor", FLAVORS)
+@pytest.mark.parametrize("resolution", [1, 2])
+@pytest.mark.parametrize("n_iterations", [-1, 3])
+def test_leiden_basic(adata_neighbors, flavor, resolution, n_iterations):
+    sc.tl.leiden(
+        adata_neighbors,
+        flavor=flavor,
+        resolution=resolution,
+        n_iterations=n_iterations,
+    )
+    assert adata_neighbors.uns["leiden"]["params"]["resolution"] == resolution
+    assert adata_neighbors.uns["leiden"]["params"]["n_iterations"] == n_iterations
+
+
+@pytest.mark.parametrize("flavor", FLAVORS)
+def test_leiden_random_state(adata_neighbors, flavor):
+    adata_1 = sc.tl.leiden(adata_neighbors, flavor=flavor, random_state=1, copy=True)
+    adata_1_again = sc.tl.leiden(
+        adata_neighbors, flavor=flavor, random_state=1, copy=True
+    )
+    adata_2 = sc.tl.leiden(adata_neighbors, flavor=flavor, random_state=50, copy=True)
+    assert (adata_1.obs["leiden"] == adata_1_again.obs["leiden"]).all()
+    assert (adata_2.obs["leiden"] != adata_1_again.obs["leiden"]).any()
+
+
+@needs.igraph
+def test_leiden_igraph_directed(adata_neighbors):
+    with pytest.raises(ValueError):
+        sc.tl.leiden(adata_neighbors, directed=True)
+
+
 @needs.leidenalg
-def test_leiden_basic(adata_neighbors):
-    sc.tl.leiden(adata_neighbors)
+@needs.igraph
+def test_leiden_equal_defaults(adata_neighbors):
+    """Ensure the two implementations are the same for the same args."""
+    leiden_alg_clustered = sc.tl.leiden(adata_neighbors, flavor="leidenalg", copy=True)
+    igraph_clustered = sc.tl.leiden(adata_neighbors, copy=True)
+    assert (
+        normalized_mutual_info_score(
+            leiden_alg_clustered.obs["leiden"], igraph_clustered.obs["leiden"]
+        )
+        > 0.9
+    )
+
+
+@needs.leidenalg
+@needs.igraph
+def test_leiden_equal_old_defaults(adata_neighbors):
+    """Ensure that the old leidenalg defaults are close enough to the current default outputs."""
+    leiden_alg_clustered = sc.tl.leiden(
+        adata_neighbors, flavor="leidenalg", directed=True, n_iterations=-1, copy=True
+    )
+    igraph_clustered = sc.tl.leiden(adata_neighbors, copy=True)
+    assert (
+        normalized_mutual_info_score(
+            leiden_alg_clustered.obs["leiden"], igraph_clustered.obs["leiden"]
+        )
+        > 0.9
+    )
+
+
+@needs.igraph
+def test_leiden_objective_function(adata_neighbors):
+    """Ensure that popping this as a `clustering_kwargs` and using it does not error out."""
+    sc.tl.leiden(adata_neighbors, objective_function="modularity")
 
 
+@needs.igraph
 @pytest.mark.parametrize(
     "clustering,key",
     [
@@ -52,6 +122,7 @@ def test_clustering_subset(adata_neighbors, clustering, key):
 
 
 @needs.louvain
+@needs.igraph
 def test_louvain_basic(adata_neighbors):
     sc.tl.louvain(adata_neighbors)
     sc.tl.louvain(adata_neighbors, use_weights=True)
@@ -60,6 +131,7 @@ def test_louvain_basic(adata_neighbors):
 
 
 @needs.louvain
+@needs.igraph
 def test_partition_type(adata_neighbors):
     import louvain
 

scanpy/tools/_leiden.py

@@ -1,6 +1,6 @@
 from __future__ import annotations
 
-from typing import TYPE_CHECKING
+from typing import TYPE_CHECKING, Literal
 
 import numpy as np
 import pandas as pd
@@ -34,14 +34,15 @@ def leiden(
     random_state: _utils.AnyRandom = 0,
     key_added: str = "leiden",
     adjacency: sparse.spmatrix | None = None,
-    directed: bool = True,
+    directed: bool = False,
     use_weights: bool = True,
-    n_iterations: int = -1,
+    n_iterations: int = 2,
     partition_type: type[MutableVertexPartition] | None = None,
     neighbors_key: str | None = None,
     obsp: str | None = None,
     copy: bool = False,
-    **partition_kwargs,
+    flavor: Literal["leidenalg", "ipgraph"] = "igraph",
+    **clustering_args,
 ) -> AnnData | None:
     """\
     Cluster cells into subgroups [Traag18]_.
@@ -96,9 +97,11 @@ def leiden(
         `obsp` and `neighbors_key` at the same time.
     copy
         Whether to copy `adata` or modify it inplace.
-    **partition_kwargs
-        Any further arguments to pass to `~leidenalg.find_partition`
-        (which in turn passes arguments to the `partition_type`).
+    flavor
+        Which package's implementation to use.
+    **clustering_args
+        Any further arguments to pass to :func:`~leidenalg.find_partition` (which in turn passes arguments to the `partition_type`)
+        or :method:`igraph.Graph.community_leiden` from `igraph`.
 
     Returns
     -------
@@ -112,13 +115,14 @@ def leiden(
         A dict with the values for the parameters `resolution`, `random_state`,
         and `n_iterations`.
     """
-    try:
-        import leidenalg
-    except ImportError:
-        raise ImportError(
-            "Please install the leiden algorithm: `conda install -c conda-forge leidenalg` or `pip3 install leidenalg`."
-        )
-    partition_kwargs = dict(partition_kwargs)
+    if flavor == "leidenalg":
+        try:
+            import leidenalg
+        except ImportError:
+            raise ImportError(
+                "Please install the leiden algorithm: `conda install -c conda-forge leidenalg` or `pip3 install leidenalg`."
+            )
+    clustering_args = dict(clustering_args)
 
     start = logg.info("running Leiden clustering")
     adata = adata.copy() if copy else adata
@@ -134,22 +138,38 @@ def leiden(
             adjacency=adjacency,
         )
     # convert it to igraph
+    if not flavor == "leidenalg" and directed:
+        raise ValueError(
+            "Cannot use igraph's leiden implemntation with a directed graph."
+        )
     g = _utils.get_igraph_from_adjacency(adjacency, directed=directed)
     # flip to the default partition type if not overriden by the user
-    if partition_type is None:
+    if partition_type is None and flavor == "leidenalg":
         partition_type = leidenalg.RBConfigurationVertexPartition
+    elif not flavor == "leidenalg" and partition_type is not None:
+        raise ValueError("Do not pass in partition_type argument when using igraph.")
     # Prepare find_partition arguments as a dictionary,
     # appending to whatever the user provided. It needs to be this way
     # as this allows for the accounting of a None resolution
     # (in the case of a partition variant that doesn't take it on input)
     if use_weights:
-        partition_kwargs["weights"] = np.array(g.es["weight"]).astype(np.float64)
-    partition_kwargs["n_iterations"] = n_iterations
-    partition_kwargs["seed"] = random_state
+        clustering_args["weights"] = (
+            "weight"
+            if not flavor == "leidenalg"
+            else np.array(g.es["weight"]).astype(np.float64)
+        )
+    clustering_args["n_iterations"] = n_iterations
+    if flavor == "leidenalg":
+        clustering_args["seed"] = random_state
     if resolution is not None:
-        partition_kwargs["resolution_parameter"] = resolution
+        clustering_args["resolution_parameter"] = resolution
     # clustering proper
-    part = leidenalg.find_partition(g, partition_type, **partition_kwargs)
+    if not flavor == "leidenalg":
+        clustering_args.setdefault("objective_function", "modularity")
+        with _utils.set_igraph_random_state(random_state):
+            part = g.community_leiden(**clustering_args)
+    else:
+        part = leidenalg.find_partition(g, partition_type, **clustering_args)
     # store output into adata.obs
     groups = np.array(part.membership)
     if restrict_to is not None: