diff --git a/src/squidpy/_docs.py b/src/squidpy/_docs.py index 12138d914..d98c7e5b3 100644 --- a/src/squidpy/_docs.py +++ b/src/squidpy/_docs.py @@ -252,6 +252,14 @@ def decorator2(obj: Any) -> Any: {_niche_mask} {_library_key} {_copy}""" +_niche_leiden_params = f"""\ +flavor + Leiden backend passed to :func:`scanpy.tl.leiden`. Defaults to ``'igraph'`` + (the ``'leidenalg'`` backend is deprecated in scanpy). +n_iterations + Number of Leiden iterations. ``-1`` iterates until convergence. +{_rng} + Every resolution is clustered with an independent rng derived from it.""" # static plotting docs _plotting_kwargs_static = """\ @@ -509,6 +517,7 @@ def decorator2(obj: Any) -> Any: niche_mask=_niche_mask, niche_min_niche_size=_niche_min_niche_size, niche_common_params=_niche_common_params, + niche_leiden_params=_niche_leiden_params, sdata_params=_sdata_params, graph_common_params=_graph_common_params, n_jobs_libraries=_n_jobs_libraries, diff --git a/src/squidpy/gr/_niche.py b/src/squidpy/gr/_niche.py index 9c4ce85de..99ce66962 100644 --- a/src/squidpy/gr/_niche.py +++ b/src/squidpy/gr/_niche.py @@ -219,6 +219,8 @@ def calculate_niche( library_key=library_key, copy=not inplace, table_key=table_key, + n_iterations=n_iterations, + rng=rng, ) elif flavor == "utag": @@ -232,6 +234,8 @@ def calculate_niche( library_key=library_key, copy=not inplace, table_key=table_key, + n_iterations=n_iterations, + rng=rng, ) elif flavor == "cellcharter": @@ -287,6 +291,10 @@ def calculate_niche_neighborhood( library_key: str | None = None, copy: bool = False, table_key: str | None = None, + *, + flavor: Literal["igraph", "leidenalg"] = "igraph", + n_iterations: int = -1, + rng: SeedLike | RNGLike | None = None, ) -> AnnData | None: """Compute niche neighborhoods using a neighborhood profile embedding and Leiden clustering. @@ -315,6 +323,7 @@ def calculate_niche_neighborhood( Weights for combining neighborhood profiles across hops. %(niche_common_params)s %(table_key)s + %(niche_leiden_params)s Returns ------- @@ -334,7 +343,9 @@ def calculate_niche_neighborhood( ) # Create instance of _LeidenClusterer using provided inputs - clusterer = _LeidenClusterer(n_neighbors, resolutions, "nhood_niche") + clusterer = _LeidenClusterer( + n_neighbors, resolutions, "nhood_niche", flavor=flavor, n_iterations=n_iterations, rng=rng + ) return _calculate_niche_custom( data, @@ -359,6 +370,10 @@ def calculate_niche_utag( library_key: str | None = None, copy: bool = False, table_key: str | None = None, + *, + flavor: Literal["igraph", "leidenalg"] = "igraph", + n_iterations: int = -1, + rng: SeedLike | RNGLike | None = None, ) -> AnnData | None: """Compute niche assignments using a UTAG-style neighborhood embedding. @@ -375,6 +390,7 @@ def calculate_niche_utag( %(niche_spatial_conn_key)s %(niche_common_params)s %(table_key)s + %(niche_leiden_params)s Returns ------- @@ -385,7 +401,9 @@ def calculate_niche_utag( embedder = _UtagEmbedder(spatial_connectivities_key) - clusterer = _LeidenClusterer(n_neighbors, resolutions, "utag_niche") + clusterer = _LeidenClusterer( + n_neighbors, resolutions, "utag_niche", flavor=flavor, n_iterations=n_iterations, rng=rng + ) return _calculate_niche_custom( data, @@ -820,21 +838,21 @@ def _validate_niche_args( "abs_nhood", "distance", "n_hop_weights", + "rng", + "n_iterations", ], "unused": [ "aggregation", "n_components", - "rng", "latent_connectivities_key", "layer_ratio", - "n_iterations", "use_weights", "use_rep", ], }, "utag": { "required": ["n_neighbors", "resolutions", "spatial_connectivities_key"], - "optional": [], + "optional": ["rng", "n_iterations"], "unused": [ "groups", "min_niche_size", @@ -844,10 +862,8 @@ def _validate_niche_args( "n_hop_weights", "aggregation", "n_components", - "rng", "latent_connectivities_key", "layer_ratio", - "n_iterations", "use_weights", "use_rep", ], @@ -1391,6 +1407,7 @@ class _LeidenClusterer(_NicheClusterer): base_colname Base name for columns added to ``adata.obs``. Resolution is appended to this to unique identify columns for each resolution. + %(niche_leiden_params)s Notes ----- @@ -1403,10 +1420,17 @@ def __init__( n_neighbors: int, resolutions: float | list[float], base_colname: str = "niche_leiden", + *, + flavor: Literal["igraph", "leidenalg"] = "igraph", + n_iterations: int = -1, + rng: SeedLike | RNGLike | None = None, ): self.n_neighbors = n_neighbors self.resolutions = resolutions if isinstance(resolutions, list) else [resolutions] self.base_colname = base_colname + self.flavor = flavor + self.n_iterations = n_iterations + self.rng = np.random.default_rng(rng) def cluster(self, adata: AnnData, embedding: NDArrayA) -> list: # first create an adata object using the embedding provided @@ -1417,18 +1441,28 @@ def cluster(self, adata: AnnData, embedding: NDArrayA) -> list: # For each resolution, apply leiden on neighborhood profile. Each cluster label equals to a niche label niche_keys = [] - for res in self.resolutions: + # every resolution is a separate clustering run, so seed each one independently + resolution_rngs = self.rng.spawn(len(self.resolutions)) + for res, res_rng in zip(self.resolutions, resolution_rngs, strict=True): niche_key = f"{self.base_colname}_res={res}" niche_keys.append(niche_key) if niche_key in adata.obs.columns: logg.info(f"Overwriting existing column '{niche_key}'") - sc.tl.leiden( - adata_embedding, - resolution=res, - key_added=niche_key, - ) + # Default to the igraph backend so niche labels are reproducible across + # versions; leidenalg is deprecated in scanpy and unstable on small graphs. + # See scverse/squidpy#1260. + leiden_kwargs: dict[str, Any] = { + "flavor": self.flavor, + "n_iterations": self.n_iterations, + "random_state": legacy_random(res_rng), + } + # scanpy's igraph backend only supports undirected graphs and errors if + # ``directed`` is left at the leidenalg default of True, so pin it to False. + if self.flavor == "igraph": + leiden_kwargs["directed"] = False + sc.tl.leiden(adata_embedding, resolution=res, key_added=niche_key, **leiden_kwargs) adata.obs[niche_key] = list( adata_embedding.obs[niche_key] diff --git a/tests/graph/test_niche.py b/tests/graph/test_niche.py index 368f6f508..2ea09a533 100644 --- a/tests/graph/test_niche.py +++ b/tests/graph/test_niche.py @@ -1,8 +1,9 @@ from __future__ import annotations import numpy as np +import pytest from anndata import AnnData -from pandas import Categorical, Series +from pandas import Series from scanpy.pp import neighbors from scipy.sparse import csr_matrix from spatialdata import SpatialData @@ -19,31 +20,42 @@ N_NEIGHBORS = 20 GROUPS = "celltype_mapped_refined" -# test if calculate_niche() gives appropriate output for dummy_adata2 for the different flavors +# Niche labels come from Leiden clustering, whose exact partition is not stable across +# igraph/leidenalg versions on tiny toy graphs (multiple equal-modularity optima). These +# tests therefore assert squidpy's behavioural contract - every cell is assigned, the +# postprocessors and library stratification behave, and a fixed seed is reproducible - +# rather than a specific (arbitrary) partition. See scverse/squidpy#1260. + + +def _assert_all_assigned(adata: AnnData, column: str) -> Series: + """Every observation receives a niche label under ``column``.""" + assert column in adata.obs.columns + niches = adata.obs[column] + assert len(niches) == adata.n_obs + assert niches.notna().all() + return niches def test_niche_calc_nhood_dummy_adata(dummy_adata2: AnnData): "Check whether niche calculation using neighborhood profile approach works as intended for dummy_adata2." - calculate_niche(dummy_adata2, flavor="neighborhood", groups="celltype", n_neighbors=3, resolutions=1.0) - assert "nhood_niche_res=1.0" in dummy_adata2.obs.columns - expected_niches = Series( - ["0", "2", "0", "2", "1", "0", "0", "1", "0", "1"], - index=["a", "b", "c", "d", "e", "f", "g", "h", "i", "j"], - name="nhood_niche_res=1.0", - ) - assert (expected_niches == dummy_adata2.obs["nhood_niche_res=1.0"]).all() + rerun = dummy_adata2.copy() + calculate_niche(dummy_adata2, flavor="neighborhood", groups="celltype", n_neighbors=3, resolutions=1.0, rng=0) + niches = _assert_all_assigned(dummy_adata2, "nhood_niche_res=1.0") + + # a fixed rng gives reproducible niches + calculate_niche(rerun, flavor="neighborhood", groups="celltype", n_neighbors=3, resolutions=1.0, rng=0) + assert (niches.to_numpy() == rerun.obs["nhood_niche_res=1.0"].to_numpy()).all() def test_niche_calc_utag_dummy_adata(dummy_adata2: AnnData): "Check whether niche calculation using utag approach works as intended for dummy_adata2." - calculate_niche(dummy_adata2, flavor="utag", n_neighbors=3, resolutions=1.0) - assert "utag_niche_res=1.0" in dummy_adata2.obs.columns - expected_niches = Series( - Categorical(["1", "0", "0", "0", "1", "0", "0", "1", "1", "0"], categories=["0", "1"]), - index=["a", "b", "c", "d", "e", "f", "g", "h", "i", "j"], - name="utag_niche_res=1.0", - ) - assert (expected_niches == dummy_adata2.obs["utag_niche_res=1.0"]).all() + rerun = dummy_adata2.copy() + calculate_niche(dummy_adata2, flavor="utag", n_neighbors=3, resolutions=1.0, rng=0) + niches = _assert_all_assigned(dummy_adata2, "utag_niche_res=1.0") + + # a fixed rng gives reproducible niches + calculate_niche(rerun, flavor="utag", n_neighbors=3, resolutions=1.0, rng=0) + assert (niches.to_numpy() == rerun.obs["utag_niche_res=1.0"].to_numpy()).all() def test_niche_calc_cellcharter_dummy_adata(dummy_adata2: AnnData): @@ -54,18 +66,12 @@ def test_niche_calc_cellcharter_dummy_adata(dummy_adata2: AnnData): calculate_niche(dummy_adata2, flavor="cellcharter", distance=2, aggregation="mean", rng=np.random.default_rng(0)) - assert "cellcharter_niche" in dummy_adata2.obs.columns - - expected_niches = Series( - Categorical([2, 6, 4, 9, 3, 0, 1, 5, 8, 7], categories=[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]), - index=["a", "b", "c", "d", "e", "f", "g", "h", "i", "j"], - name="cellcharter_niche", - ) - assert (expected_niches == dummy_adata2.obs["cellcharter_niche"]).all() + _assert_all_assigned(dummy_adata2, "cellcharter_niche") def test_niche_calc_spatialleiden_dummy_adata(dummy_adata2: AnnData): "Check whether niche calculation using spatialleiden approach works as intended for dummy_adata2." + pytest.importorskip("spatialleiden") # need the latent_connectivities_key, meaning have to run the graph construction neighbors(dummy_adata2, n_neighbors=3, use_rep="X") @@ -79,14 +85,7 @@ def test_niche_calc_spatialleiden_dummy_adata(dummy_adata2: AnnData): rng=np.random.default_rng(0), ) - assert "spatialleiden_res=1.0" in dummy_adata2.obs.columns - expected_niches = Series( - Categorical([0, 0, 0, 0, 1, 1, 1, 2, 2, 2], categories=[0, 1, 2]), - index=["a", "b", "c", "d", "e", "f", "g", "h", "i", "j"], - name="spatialleiden_res=1.0", - ) - - assert (expected_niches == dummy_adata2.obs["spatialleiden_res=1.0"]).all() + _assert_all_assigned(dummy_adata2, "spatialleiden_res=1.0") # rng handling @@ -146,65 +145,27 @@ def test_niche_calc_library_key_dummy_adata(dummy_adata2: AnnData): "Check whether niche calculation when library_key is supplied works as intended for dummy_adata2." # add library_key information in dummy_adata - dummy_adata2.obs["batch"] = [ - "batch1", - "batch1", - "batch1", - "batch1", - "batch1", - "batch2", - "batch2", - "batch2", - "batch2", - "batch2", - ] + dummy_adata2.obs["batch"] = ["batch1"] * 5 + ["batch2"] * 5 calculate_niche( dummy_adata2, flavor="neighborhood", groups="celltype", n_neighbors=3, resolutions=1.5, library_key="batch" ) - assert "nhood_niche_res=1.5" in dummy_adata2.obs.columns - - expected_niches = Series( - [ - "lib=batch1_1", - "lib=batch1_0", - "lib=batch1_2", - "lib=batch1_0", - "lib=batch1_1", - "lib=batch2_2", - "lib=batch2_0", - "lib=batch2_1", - "lib=batch2_0", - "lib=batch2_1", - ], - index=["a", "b", "c", "d", "e", "f", "g", "h", "i", "j"], - name="nhood_niche_res=1.5", - dtype=str, - ) - - assert (expected_niches == dummy_adata2.obs["nhood_niche_res=1.5"]).all() + niches = _assert_all_assigned(dummy_adata2, "nhood_niche_res=1.5") + # niches are computed per library and prefixed with the originating library + for cell, label in niches.items(): + assert label.startswith(f"lib={dummy_adata2.obs['batch'][cell]}_") def test_niche_calc_spatialleiden_library_key_dummy_adata(dummy_adata2: AnnData): "Check whether niche calculation for spatialleiden works as intended for dummy_adata2 when library_key is supplied." + pytest.importorskip("spatialleiden") # need the latent_connectivities_key, meaning have to run the graph construction neighbors(dummy_adata2, n_neighbors=3, use_rep="X") # add library_key information in dummy_adata - dummy_adata2.obs["batch"] = [ - "batch1", - "batch1", - "batch1", - "batch1", - "batch1", - "batch2", - "batch2", - "batch2", - "batch2", - "batch2", - ] + dummy_adata2.obs["batch"] = ["batch1"] * 5 + ["batch2"] * 5 calculate_niche( dummy_adata2, @@ -216,27 +177,10 @@ def test_niche_calc_spatialleiden_library_key_dummy_adata(dummy_adata2: AnnData) rng=np.random.default_rng(0), ) - assert "spatialleiden_res=1.0" in dummy_adata2.obs.columns - - expected_niches = Series( - [ - "lib=batch1_1", - "lib=batch1_0", - "lib=batch1_0", - "lib=batch1_1", - "lib=batch1_0", - "lib=batch2_1", - "lib=batch2_1", - "lib=batch2_0", - "lib=batch2_0", - "lib=batch2_0", - ], - index=["a", "b", "c", "d", "e", "f", "g", "h", "i", "j"], - name="spatialleiden_res=1.0", - dtype=str, - ) - - assert (expected_niches == dummy_adata2.obs["spatialleiden_res=1.0"]).all() + niches = _assert_all_assigned(dummy_adata2, "spatialleiden_res=1.0") + # niches are computed per library and prefixed with the originating library + for cell, label in niches.items(): + assert label.startswith(f"lib={dummy_adata2.obs['batch'][cell]}_") def test_niche_calc_nhood_multipostprocessor_dummy_adata(dummy_adata2: AnnData): @@ -254,13 +198,12 @@ def test_niche_calc_nhood_multipostprocessor_dummy_adata(dummy_adata2: AnnData): mask=mask, min_niche_size=3, ) - assert "nhood_niche_res=1.0" in dummy_adata2.obs.columns - expected_niches = Series( - ["not_a_niche", "not_a_niche", "0", "not_a_niche", "1", "0", "0", "1", "0", "1"], - index=["a", "b", "c", "d", "e", "f", "g", "h", "i", "j"], - name="nhood_niche_res=1.0", - ) - assert (expected_niches == dummy_adata2.obs["nhood_niche_res=1.0"]).all() + niches = _assert_all_assigned(dummy_adata2, "nhood_niche_res=1.0") + # masked-out observations are never assigned to a real niche + assert (niches[["a", "b"]] == "not_a_niche").all() + # every real niche respects the requested minimum size + real = niches[niches != "not_a_niche"] + assert (real.value_counts() >= 3).all() def test_niche_calc_nhood_dummy_sdata(dummy_adata2: AnnData): @@ -276,16 +219,7 @@ def test_niche_calc_nhood_dummy_sdata(dummy_adata2: AnnData): calculate_niche(sdata, flavor="neighborhood", groups="celltype", n_neighbors=3, resolutions=1.0, table_key="adata") - assert "nhood_niche_res_1.0" in sdata["adata"].obs.columns - - expected_niches = Series( - ["0", "2", "0", "2", "1", "0", "0", "1", "0", "1"], - index=["a", "b", "c", "d", "e", "f", "g", "h", "i", "j"], - name="nhood_niche_res_1.0", - dtype=str, - ) - - assert (expected_niches == sdata["adata"].obs["nhood_niche_res_1.0"]).all() + _assert_all_assigned(sdata["adata"], "nhood_niche_res_1.0") # older tests