diff --git a/hatch.toml b/hatch.toml index 1614be834..f7028839e 100644 --- a/hatch.toml +++ b/hatch.toml @@ -13,6 +13,10 @@ scripts.download = "python ./.scripts/ci/download_data.py {args}" [envs.hatch-test] dependency-groups = ["test"] +# to pass spatialleiden-using tests in test_niche.py +features = [ + "leiden" +] extra-dependencies = ["diff-cover"] matrix = [ { deps = ["stable"], python = ["3.12", "3.13", "3.14"] }, diff --git a/src/squidpy/gr/__init__.py b/src/squidpy/gr/__init__.py index 6216159d3..50d3a1046 100644 --- a/src/squidpy/gr/__init__.py +++ b/src/squidpy/gr/__init__.py @@ -20,7 +20,22 @@ interaction_matrix, nhood_enrichment, ) -from squidpy.gr._niche import calculate_niche +from squidpy.gr._niche import ( + CellcharterEmbedder, + GMMClusterer, + LeidenClusterer, + MaskPostprocessor, + MinNicheSizePostprocessor, + NhoodProfileEmbedder, + RenamePostprocessor, + UtagEmbedder, + calculate_niche, + calculate_niche_cellcharter, + calculate_niche_custom, + calculate_niche_neighborhood, + calculate_niche_spatialleiden, + calculate_niche_utag, +) from squidpy.gr._ppatterns import co_occurrence, spatial_autocorr from squidpy.gr._ripley import ripley from squidpy.gr._sepal import sepal @@ -38,11 +53,24 @@ "spatial_neighbors_radius", "spatial_neighbors_delaunay", "spatial_neighbors_grid", + "calculate_niche", + "calculate_niche_neighborhood", + "calculate_niche_utag", + "calculate_niche_cellcharter", + "calculate_niche_spatialleiden", + "calculate_niche_custom", + "NhoodProfileEmbedder", + "UtagEmbedder", + "CellcharterEmbedder", + "LeidenClusterer", + "GMMClusterer", + "MinNicheSizePostprocessor", + "MaskPostprocessor", + "RenamePostprocessor", "ligrec", "centrality_scores", "interaction_matrix", "nhood_enrichment", - "calculate_niche", "co_occurrence", "spatial_autocorr", "ripley", diff --git a/src/squidpy/gr/_niche.py b/src/squidpy/gr/_niche.py index 07db78246..d0b902390 100644 --- a/src/squidpy/gr/_niche.py +++ b/src/squidpy/gr/_niche.py @@ -1,6 +1,7 @@ from __future__ import annotations -import contextlib +import warnings +from abc import abstractmethod from typing import Any, Literal import anndata as ad @@ -9,21 +10,37 @@ import scanpy as sc import scipy.sparse as sps from anndata import AnnData -from numpy.typing import NDArray -from scipy.sparse import coo_matrix, hstack, issparse, spdiags -from scipy.spatial import distance -from sklearn.metrics import f1_score +from scipy.sparse import coo_matrix, hstack, issparse, lil_matrix, spdiags from sklearn.mixture import GaussianMixture from sklearn.preprocessing import normalize -from spatialdata import SpatialData +from spatialdata import SpatialData, sanitize_table from spatialdata._logging import logger as logg from squidpy._constants._constants import NicheDefinitions from squidpy._docs import d, inject_docs +from squidpy._utils import NDArrayA from squidpy._validators import assert_isinstance, assert_key_in_adata, assert_one_of from squidpy.gr._utils import extract_adata_if_sdata -__all__ = ["calculate_niche"] +__all__ = [ + "calculate_niche", + "calculate_niche_neighborhood", + "calculate_niche_utag", + "calculate_niche_cellcharter", + "calculate_niche_spatialleiden", + "calculate_niche_custom", + "NicheEmbedder", + "NhoodProfileEmbedder", + "UtagEmbedder", + "CellcharterEmbedder", + "NicheClusterer", + "LeidenClusterer", + "GMMClusterer", + "NichePostprocessor", + "MinNicheSizePostprocessor", + "MaskPostprocessor", + "RenamePostprocessor", +] @d.dedent @@ -58,6 +75,16 @@ def calculate_niche( Calculate niches (spatial clusters) based on a user-defined method in 'flavor'. The resulting niche labels with be stored in 'adata.obs'. + .. deprecated:: 1.8.3 + ``calculate_niche`` is deprecated and will be removed in squidpy + v1.9.0. Use one of the flavor-specific functions instead: + + - ``calculate_niche_neighborhood`` + - ``calculate_niche_utag`` + - ``calculate_niche_cellcharter`` + - ``calculate_niche_spatialleiden`` + - ``calculate_niche_custom`` + Parameters ---------- %(adata)s @@ -74,6 +101,9 @@ def calculate_niche( mask Boolean array to filter cells which won't get assigned to a niche. Note that if you want to exclude these cells during neighborhood calculation already, you should subset your AnnData table before running 'sq.gr.spatial_neigbors'. + Mask can look like the following. Here, the index values would correspond to adata.obs.index. + The entries that are False are the ones ignored. + mask = Series([False, False, True], index = ["a", "b", "c"]) groups Groups based on which to calculate neighborhood profile (E.g. columns of cell type annotations in adata.obs). Required if flavor == `{fla.NEIGHBORHOOD.s!r}`. @@ -170,635 +200,630 @@ def calculate_niche( if resolutions is None: resolutions = [0.5] - orig_adata = extract_adata_if_sdata(data, table_key=table_key) - adata = orig_adata.copy() - - assert_key_in_adata( - adata, - spatial_connectivities_key, - attr="obsp", - extra_msg="If you haven't computed a spatial neighborhood graph yet, use `sq.gr.spatial_neighbors`.", - ) - - if flavor == "spatialleiden": - assert_key_in_adata( - adata, - latent_connectivities_key, - attr="obsp", - extra_msg="If you haven't computed a latent neighborhood graph yet, use `sc.pp.neighbors`.", - ) - - result_columns = _get_result_columns( - flavor=flavor, - resolutions=resolutions, - library_key=None, - libraries=None, - ) - - if library_key is not None: - assert_key_in_adata(adata, library_key, attr="obs") - - logg.info(f"Stratifying by library_key '{library_key}'") - - for col in result_columns: - adata.obs[col] = "not_a_niche" - - for lib_id in adata.obs[library_key].unique(): - logg.info(f"Processing library '{lib_id}'") - - lib_indices = adata.obs[adata.obs[library_key] == lib_id].index - - if len(lib_indices) == 0: - logg.warning(f"Library '{lib_id}' contains no cells, skipping") - continue - - lib_adata = adata[lib_indices].copy() - - lib_mask = None - if mask is not None: - lib_mask = mask[mask.index.isin(lib_indices)] - - lib_result = calculate_niche( - lib_adata, - flavor=flavor, - library_key=None, - mask=lib_mask, - groups=groups, - n_neighbors=n_neighbors, - resolutions=None if flavor == "cellcharter" else resolutions, - min_niche_size=min_niche_size, - scale=scale, - abs_nhood=abs_nhood, - distance=None if flavor == "utag" else distance, - n_hop_weights=n_hop_weights, - aggregation=aggregation, - n_components=n_components, - random_state=random_state, - spatial_connectivities_key=spatial_connectivities_key, - latent_connectivities_key=latent_connectivities_key, - layer_ratio=layer_ratio, - n_iterations=n_iterations, - use_weights=use_weights, - inplace=False, - ) - - for col in result_columns: - if col in lib_result.obs.columns: - prefixed_values = lib_result.obs[col].apply( - lambda x, lib=lib_id: f"lib={lib}_{x}" if x != "not_a_niche" else x - ) - - adata.obs.loc[lib_indices, col] = prefixed_values.values - - else: - _calculate_niches( - adata, - mask, - flavor, + if flavor == "neighborhood": + return calculate_niche_neighborhood( + data, groups, - n_neighbors, resolutions, - min_niche_size, + n_neighbors, + spatial_connectivities_key, scale, - abs_nhood, distance, + abs_nhood, n_hop_weights, - aggregation, - n_components, - random_state, - spatial_connectivities_key, - latent_connectivities_key, - layer_ratio, - n_iterations, - use_weights, - use_rep, + min_niche_size, + mask, + library_key, + inplace, + table_key, ) - if not inplace: - return adata - # For SpatialData, update the table directly - if isinstance(data, SpatialData): - data.tables[table_key] = adata - else: - # For AnnData, copy results back to original object - for col in result_columns: - if col in orig_adata.obs.columns: - logg.info(f"Overwriting existing column '{col}'") - with contextlib.suppress(KeyError): - del orig_adata.obs[col] - if f"{col}_colors" in orig_adata.uns.keys(): - with contextlib.suppress(KeyError): - del orig_adata.uns[f"{col}_colors"] - - orig_adata.obs[col] = adata.obs[col] - - return None - - -def _get_result_columns( - flavor: str, - resolutions: float | tuple[float, float] | list[float | tuple[float, float]], - library_key: str | None, - libraries: list[str] | None, -) -> list[str]: - """Get the column names that will be populated based on flavor and resolutions.""" - - library_str = f"_{library_key}" if library_key is not None else "" - - if flavor == "cellcharter": - base_column = "cellcharter_niche" - if library_key is None: - return [base_column] - elif libraries is not None and len(libraries) > 0: - return [f"{base_column}_{lib}" for lib in libraries] - - # For neighborhood, utag and spatialleiden, we need to handle resolutions - if not isinstance(resolutions, list): - resolutions = [resolutions] - - if flavor == "neighborhood": - prefix = f"nhood_niche{library_str}" elif flavor == "utag": - prefix = f"utag_niche{library_str}" - elif flavor == "spatialleiden": - prefix = f"spatialleiden{library_str}" - - if library_key is None: - return [f"{prefix}_res={res}" for res in resolutions] - else: - assert isinstance(libraries, list) # for mypy - return [f"{prefix}_{lib}_res={res}" for lib in libraries for res in resolutions] - - -def _calculate_niches( - adata: AnnData, - mask: pd.core.series.Series | None, - flavor: str, - groups: str | None, - n_neighbors: int | None, - resolutions: float | tuple[float, float] | list[float | tuple[float, float]], - min_niche_size: int | None, - scale: bool, - abs_nhood: bool, - distance: int, - n_hop_weights: list[float] | None, - aggregation: str | None, - n_components: int | None, - random_state: int, - spatial_connectivities_key: str, - latent_connectivities_key: str, - layer_ratio: float, - n_iterations: int, - use_weights: bool | tuple[bool, bool], - use_rep: str | None, -) -> None: - """Calculate niches using the specified flavor and parameters.""" - if flavor == "neighborhood": - assert isinstance(resolutions, float | list) - _get_nhood_profile_niches( - adata, - mask, - groups, - n_neighbors, + return calculate_niche_utag( + data, resolutions, - min_niche_size, - scale, - abs_nhood, - distance, - n_hop_weights, + n_neighbors, spatial_connectivities_key, + min_niche_size, + mask, + library_key, + inplace, + table_key, ) - elif flavor == "utag": - assert isinstance(resolutions, float | list) - _get_utag_niches(adata, n_neighbors, resolutions, spatial_connectivities_key) + elif flavor == "cellcharter": - assert isinstance(aggregation, str) # for mypy - assert isinstance(n_components, int) # for mypy - _get_cellcharter_niches( - adata, + return calculate_niche_cellcharter( + data, distance, aggregation, - n_components, random_state, spatial_connectivities_key, + n_components, use_rep, + min_niche_size, + mask, + library_key, + inplace, + table_key, ) + elif flavor == "spatialleiden": - _get_spatialleiden_domains( - adata, - spatial_connectivities_key, - latent_connectivities_key, + return calculate_niche_spatialleiden( + data, resolutions, + latent_connectivities_key, + spatial_connectivities_key, layer_ratio, - use_weights, n_iterations, + use_weights, random_state, + min_niche_size, + mask, + prefix=None, + library_key=library_key, + inplace=inplace, + table_key=table_key, ) + return -def _get_nhood_profile_niches( - adata: AnnData, - mask: pd.core.series.Series | None, - groups: str | None, - n_neighbors: int | None, - resolutions: float | tuple[float, float] | list[float | tuple[float, float]], - min_niche_size: int | None, - scale: bool, - abs_nhood: bool, - distance: int, - n_hop_weights: list[float] | None, - spatial_connectivities_key: str, -) -> None: - """ - adapted from https://github.com/immunitastx/monkeybread/blob/main/src/monkeybread/calc/_neighborhood_profile.py - """ - - adata_masked = adata - # get obs x neighbor matrix from sparse matrix - matrix = adata_masked.obsp[spatial_connectivities_key].tocoo() +@d.dedent +def calculate_niche_neighborhood( + data: AnnData | SpatialData, + groups: str, + resolutions: float | list[float], + n_neighbors: int = 15, + spatial_connectivities_key: str = "spatial_connectivities", + scale: bool = True, + distance: int = 1, + abs_nhood: bool = False, + n_hop_weights: list[float] | None = None, + min_niche_size: int | None = None, + mask: pd.Series | None = None, + library_key: str | None = None, + inplace: bool = True, + table_key: str | None = None, +) -> AnnData | None: + """Compute niche neighborhoods using a neighborhood profile embedding and Leiden clustering. - # get obs x category matrix where each column is the absolute/relative frequency of a category in the neighborhood - nhood_profile = _calculate_neighborhood_profile(adata_masked, groups, matrix, abs_nhood) + This is a high-level convenience wrapper that constructs a + :class:`NhoodProfileEmbedder`, a :class:`LeidenClusterer`, and optional + postprocessors to compute niche assignments for each observation. - # Additionally use n-hop neighbors if distance > 1. This sums up the (weighted) neighborhood profiles of all n-hop neighbors. - if distance > 1: - n_hop_adjacency_matrix = adata_masked.obsp[spatial_connectivities_key].copy() - # if no weights are provided, use 1 for all n_hop neighbors - if n_hop_weights is None: - n_hop_weights = [1] * distance - # if weights are provided, start with applying weight to the original neighborhood profile - elif len(n_hop_weights) < distance: - # Extend weights if too few provided - n_hop_weights = n_hop_weights + [n_hop_weights[-1]] * (distance - len(n_hop_weights)) - logg.debug(f"Extended weights to match distance: {n_hop_weights}") + Parameters + ---------- + %(adata)s + groups + Column in ``adata.obs`` defining categorical groups (e.g. cell types) + used to compute neighborhood composition profiles. + n_neighbors + Number of neighbors used when constructing the graph for Leiden clustering. + resolutions + Resolution parameter(s) for Leiden clustering. Can be a single float or a list. + spatial_connectivities_key + Key in ``adata.obsp`` containing the spatial connectivity matrix. + scale + Whether to z-score the neighborhood profile prior to clustering. + distance + Number of hops to consider when constructing neighborhood profiles. + Values greater than ``1`` incorporate higher-order neighbors. + abs_nhood + If ``True``, use absolute counts; otherwise normalize to proportions. + n_hop_weights + Weights for combining neighborhood profiles across hops. + min_niche_size + Minimum number of observations required for a niche; smaller niches are filtered. + mask + Boolean mask or index specifying observations to exclude from niche assignment. + %(library_key)s + inplace + Whether to modify ``adata`` in place. + %(table_key)s - # Apply first weight to base profile - weighted_profile = n_hop_weights[0] * nhood_profile + Returns + ------- + If ``inplace = True``, modifies ``adata`` in place and returns ``None``. + Otherwise, returns a copy of ``adata`` with niche annotations added to ``.obs``. - # Calculate higher-order hop profiles - n_hop_adjacency_matrix = adata_masked.obsp[spatial_connectivities_key].copy() + See Also + -------- + calculate_niche_custom : Lower-level API for custom embedding, clustering, and postprocessing. + NhoodProfileEmbedder : Default embedding strategy based on neighborhood composition. + LeidenClusterer : Default clustering strategy. + """ - # get n_hop neighbor adjacency matrices by multiplying the original adjacency matrix with itself n times and get corresponding neighborhood profiles. - for n_hop in range(1, distance): - logg.debug(f"Calculating {n_hop + 1}-hop neighbors") - # Multiply adjacency matrix by itself to get n+1 hop adjacency - n_hop_adjacency_matrix = n_hop_adjacency_matrix @ adata_masked.obsp[spatial_connectivities_key] - matrix = n_hop_adjacency_matrix.tocoo() + # Create instance of NhoodProfileEmbedder using provided inputs + embedder = NhoodProfileEmbedder( + groups, + spatial_connectivities_key, + scale, + distance, + abs_nhood, + n_hop_weights, + ) - # Calculate and add weighted profile - hop_profile = _calculate_neighborhood_profile(adata_masked, groups, matrix, abs_nhood) - weighted_profile += n_hop_weights[n_hop] * hop_profile + # Create instance of LeidenClusterer using provided inputs + clusterer = LeidenClusterer(n_neighbors, resolutions, "nhood_niche") - if not abs_nhood: - weighted_profile = weighted_profile / sum(n_hop_weights) + # generate the list of postprocessor objects using the supplied args + postprocessors_list = [] + if mask is not None: + mask_postprocessor = MaskPostprocessor(mask) + postprocessors_list.append(mask_postprocessor) + if min_niche_size is not None: + min_niche_size_postprocessor = MinNicheSizePostprocessor(min_niche_size) + postprocessors_list.append(min_niche_size_postprocessor) - nhood_profile = weighted_profile + return calculate_niche_custom(data, embedder, clusterer, postprocessors_list, library_key, inplace, table_key) - # create AnnData object from neighborhood profile to perform scanpy functions - # Use .to_numpy(copy=True) to ensure the array is writeable (required for pandas CoW compatibility) - # Preserve the DataFrame index for later matching with adata_masked - adata_neighborhood = ad.AnnData(X=nhood_profile.to_numpy(copy=True), obs=pd.DataFrame(index=nhood_profile.index)) - # reason for scaling see https://monkeybread.readthedocs.io/en/latest/notebooks/tutorial.html#niche-analysis - if scale: - sc.pp.scale(adata_neighborhood, zero_center=True) +@d.dedent +def calculate_niche_utag( + data: AnnData | SpatialData, + resolutions: float | list[float], + n_neighbors: int = 15, + spatial_connectivities_key: str = "spatial_connectivities", + min_niche_size: int | None = None, + mask: pd.Series | None = None, + library_key: str | None = None, + inplace: bool = True, + table_key: str | None = None, +) -> AnnData | None: + """Compute niche assignments using a UTAG-style neighborhood embedding. - # mask obs to exclude cells for which no niche shall be assigned - if mask is not None: - mask = mask[mask.index.isin(adata_neighborhood.obs.index)] - adata_neighborhood = adata_neighborhood[mask] + This wrapper constructs a :class:`UtagEmbedder`, a + :class:`LeidenClusterer`, and optional postprocessors to generate niche + labels from spatial neighborhoods. - # required for leiden clustering (note: no dim reduction performed in original implementation) - sc.pp.neighbors(adata_neighborhood, n_neighbors=n_neighbors, use_rep="X") + Parameters + ---------- + %(adata)s + n_neighbors + Number of neighbors used when constructing the graph for Leiden clustering. + resolutions + Resolution parameter(s) for Leiden clustering. Can be a single float or a list. + spatial_connectivities_key + Key in ``adata.obsp`` containing the spatial connectivity matrix. + min_niche_size + Minimum number of observations required for a niche; smaller niches are filtered. + mask + Boolean mask or index specifying observations to exclude from niche assignment. + %(library_key)s + inplace + Whether to modify ``adata`` in place. + %(table_key)s - resolutions = resolutions if isinstance(resolutions, list) else [resolutions] + Returns + ------- + If ``inplace = True``, modifies ``adata`` in place and returns ``None``. + Otherwise, returns a copy of ``adata`` with niche annotations added to ``.obs``. - # For each resolution, apply leiden on neighborhood profile. Each cluster label equals to a niche label - for res in resolutions: - niche_key = f"nhood_niche_res={res}" + See Also + -------- + UtagEmbedder : Embedding strategy based on UTAG neighborhood feature propagation. + LeidenClusterer : Leiden clustering backend used to assign niches. + calculate_niche_custom : Lower-level API for custom niche pipelines. + """ - if niche_key in adata_masked.obs.columns: - del adata_masked.obs[niche_key] + embedder = UtagEmbedder(spatial_connectivities_key) - if f"{niche_key}_colors" in adata_masked.uns.keys(): - del adata_masked.uns[f"{niche_key}_colors"] - # print(adata_masked.obs[niche_key]) + clusterer = LeidenClusterer(n_neighbors, resolutions, "utag_niche") - sc.tl.leiden( - adata_neighborhood, - resolution=res, - key_added=niche_key, - ) + # generate the list of postprocessor objects using the supplied args + postprocessors_list = [] + if mask is not None: + mask_postprocessor = MaskPostprocessor(mask) + postprocessors_list.append(mask_postprocessor) + if min_niche_size is not None: + min_niche_size_postprocessor = MinNicheSizePostprocessor(min_niche_size) + postprocessors_list.append(min_niche_size_postprocessor) - adata_masked.obs[niche_key] = "not_a_niche" + return calculate_niche_custom(data, embedder, clusterer, postprocessors_list, library_key, inplace, table_key) - neighborhood_clusters = dict(zip(adata_neighborhood.obs.index, adata_neighborhood.obs[niche_key], strict=False)) - mask_indices = adata_masked.obs.index - adata_masked.obs.loc[mask_indices, niche_key] = [ - neighborhood_clusters.get(idx, "not_a_niche") for idx in mask_indices - ] +@d.dedent +def calculate_niche_cellcharter( + data: AnnData | SpatialData, + distance: int = 2, + aggregation: str = "mean", + random_state: int = 0, + spatial_connectivities_key: str = "spatial_connectivities", + n_components: int = 10, + use_rep: str | None = None, + min_niche_size: int | None = None, + mask: pd.Series | None = None, + library_key: str | None = None, + inplace: bool = True, + table_key: str | None = None, +) -> AnnData | None: + """Compute niche assignments using a CellCharter-style aggregation embedding. - # filter niches with n_cells < min_niche_size - if min_niche_size is not None: - counts_by_niche = adata_masked.obs[niche_key].value_counts() - to_filter = counts_by_niche[counts_by_niche < min_niche_size].index - adata_masked.obs[niche_key] = adata_masked.obs[niche_key].apply( - lambda x, to_filter=to_filter: "not_a_niche" if x in to_filter else x - ) - adata_masked.obs[niche_key] = adata_masked.obs.index.map(adata_masked.obs[niche_key]).fillna("not_a_niche") + This wrapper constructs a :class:`CellcharterEmbedder`, a + :class:`GMMClusterer`, and optional postprocessors to generate niche labels + from spatial neighborhoods. - return + Parameters + ---------- + %(adata)s + distance + Number of neighborhood hops to aggregate when building the embedding. + aggregation + Aggregation mode used for neighborhood features, typically ``"mean"`` or + ``"variance"``. + random_state + Random seed used by the Gaussian mixture clustering step. + spatial_connectivities_key + Key in ``adata.obsp`` containing the spatial connectivity matrix. + n_components + Number of embedding components to retain when ``use_rep`` is provided, + or number of mixture components used by the clusterer. + use_rep + Key in ``adata.obsm`` pointing to a precomputed representation to use + instead of deriving a spatially aggregated embedding. + min_niche_size + Minimum number of observations required for a niche; smaller niches are filtered. + mask + Boolean mask or index specifying observations to exclude from niche assignment. + %(library_key)s + inplace + Whether to modify ``adata`` in place. + %(table_key)s + Returns + ------- + If ``inplace = True``, modifies ``adata`` in place and returns ``None``. + Otherwise, returns a copy of ``adata`` with niche annotations added to ``.obs``. -def _get_utag_niches( - adata: AnnData, - n_neighbors: int | None, - resolutions: float | tuple[float, float] | list[float | tuple[float, float]], - spatial_connectivities_key: str, -) -> None: - """ - Adapted from https://github.com/ElementoLab/utag/blob/main/utag/segmentation.py + See Also + -------- + CellcharterEmbedder : Embedding strategy inspired by CellCharter. + GMMClusterer : Gaussian-mixture clustering backend used to assign niches. + calculate_niche_custom : Lower-level API for custom niche pipelines. """ - new_feature_matrix = _utag(adata, normalize_adj=True, spatial_connectivity_key=spatial_connectivities_key) - adata_utag = ad.AnnData(X=new_feature_matrix) - sc.tl.pca(adata_utag) # note: unlike with flavor 'neighborhood' dim reduction is performed here - sc.pp.neighbors(adata_utag, n_neighbors=n_neighbors, use_rep="X_pca") + embedder = CellcharterEmbedder(distance, aggregation, spatial_connectivities_key, n_components, use_rep) - if not isinstance(resolutions, list): - resolutions = [resolutions] - # For each resolution, apply leiden on neighborhood profile. Each cluster label equals to a niche label - for res in resolutions: - sc.tl.leiden(adata_utag, resolution=res, key_added=f"utag_niche_res={res}") - adata.obs[f"utag_niche_res={res}"] = adata_utag.obs[f"utag_niche_res={res}"].values + clusterer = GMMClusterer(n_components, random_state, base_colname="cellcharter_niche") - return + # generate the list of postprocessor objects using the supplied args + postprocessors_list = [] + if mask is not None: + mask_postprocessor = MaskPostprocessor(mask) + postprocessors_list.append(mask_postprocessor) + if min_niche_size is not None: + min_niche_size_postprocessor = MinNicheSizePostprocessor(min_niche_size) + postprocessors_list.append(min_niche_size_postprocessor) + return calculate_niche_custom(data, embedder, clusterer, postprocessors_list, library_key, inplace, table_key) -def _get_cellcharter_niches( - adata: AnnData, - distance: int, - aggregation: str, - n_components: int, - random_state: int, - spatial_connectivities_key: str, - use_rep: str | None = None, -) -> None: - """adapted from https://github.com/CSOgroup/cellcharter/blob/main/src/cellcharter/gr/_aggr.py - and https://github.com/CSOgroup/cellcharter/blob/main/src/cellcharter/tl/_gmm.py""" - - adjacency_matrix = adata.obsp[spatial_connectivities_key] - layers = list(range(distance + 1)) - - aggregated_matrices = [] - adj_hop = _setdiag(adjacency_matrix, 0) # Remove self-loops, set diagonal to 0 - adj_visited = _setdiag(adjacency_matrix.copy(), 1) # Track visited neighbors - for k in layers: - if k == 0: - # get original count matrix (not aggregated) - aggregated_matrices.append(adata.X) - else: - # get count and adjacency matrix for k-hop (neighbor of neighbor of neighbor ...) and aggregate them - if k > 1: - adj_hop, adj_visited = _hop(adj_hop, adjacency_matrix, adj_visited) - adj_hop_norm = _normalize(adj_hop) - aggregated_matrix = _aggregate(adata, adj_hop_norm, aggregation) - aggregated_matrices.append(aggregated_matrix) - - concatenated_matrix = hstack(aggregated_matrices) # Stack all matrices horizontally - arr = concatenated_matrix.toarray() # Densify - - if use_rep is not None: - # Use provided embedding from adata.obsm - assert_key_in_adata(adata, use_rep, attr="obsm") - embedding = adata.obsm[use_rep] - # Ensure embedding has the right number of components - if embedding.shape[1] < n_components: - raise ValueError( - f"Embedding has {embedding.shape[1]} components, but n_components={n_components}. Please provide an embedding with at least {n_components} components." - ) - # Use only the first n_components - embedding = embedding[:, :n_components] - else: - logg.warning( - "CellCharter recommends to use a dimensionality reduced embedding of the data, e.g. a scVI embedding. Since 'use_rep' is not provided, PCA will be used as proxy - performance may be suboptimal." - ) - arr_ad = ad.AnnData(X=arr) - sc.tl.pca(arr_ad) - embedding = arr_ad.obsm["X_pca"] +@d.dedent +def calculate_niche_spatialleiden( + data: AnnData | SpatialData, + resolutions: float | tuple[float, float] | list[float | tuple[float, float]], + latent_connectivities_key: str = "connectivities", + spatial_connectivities_key: str = "spatial_connectivities", + layer_ratio: float = 1.0, + n_iterations: int = -1, + use_weights: bool | tuple[bool, bool] = True, + random_state: int = 42, + min_niche_size: int | None = None, + mask: pd.Series | None = None, + prefix: str | None = None, + library_key: str | None = None, + inplace: bool = True, + table_key: str | None = None, +) -> AnnData | None: + """Compute niche assignments using the SpatialLeiden algorithm. + + This is a wrapper around :func:`spatialleiden.multiplex_leiden` that + uses :class:`AnnData` as input and works with two layers; one latent + space and one spatial layer. + Adapted from https://github.com/HiDiHlabs/SpatialLeiden/. - # cluster concatenated matrix with GMM, each cluster label equals to a niche label - niches = _get_GMM_clusters(embedding, n_components, random_state) + Parameters + ---------- + %(adata)s + latent_connectivities_key + Key in ``adata.obsp`` containing the latent-space connectivity matrix. + spatial_connectivities_key + Key in ``adata.obsp`` containing the spatial connectivity matrix. + resolutions + Resolution parameter(s) for the Leiden optimization. Can be a single + float or a list of floats. + layer_ratio + Relative weight assigned to the latent and spatial layers. + n_iterations + Number of optimization iterations used by SpatialLeiden. + use_weights + Whether to use edge weights during clustering. + random_state + Random seed passed to the SpatialLeiden routine. + min_niche_size + Minimum number of observations required for a niche; smaller niches are filtered. + mask + Boolean mask or index specifying observations to exclude from niche assignment. + prefix + Prefix added to niche labels produced by SpatialLeiden. + When stratifying by ``library_key``, a library-specific prefix is added + automatically (something like "lib="). + %(library_key)s + inplace + Whether to modify ``adata`` in place. + %(table_key)s - adata.obs["cellcharter_niche"] = pd.Categorical(niches) - return + Returns + ------- + If ``inplace = True``, modifies ``adata`` in place and returns ``None``. + Otherwise, returns a copy of ``adata`` with niche annotations added to ``.obs``. + Notes + ----- + If ``library_key`` is provided, clustering is performed independently for + each library and the results are merged back into the parent object. -def _calculate_neighborhood_profile( - adata: AnnData, - groups: str | None, - matrix: coo_matrix, - abs_nhood: bool, -) -> pd.DataFrame: - """ - Returns an obs x category matrix where each column is the absolute/relative frequency of a category in the neighborhood + See Also + -------- + spatialleiden.multiplex_leiden : SpatialLeiden implementation used internally. """ - nonzero_indices = np.split(matrix.col, matrix.row.searchsorted(np.arange(1, matrix.shape[0]))) - neighbor_matrix = pd.DataFrame(nonzero_indices) + try: + import spatialleiden as sl + except ImportError as e: + msg = "Please install the spatialleiden algorithm: `pip install squidpy[leiden]` or `conda install bioconda::spatialleiden` or `pip install spatialleiden`." + raise ImportError(msg) from e + + # obtain adata if data was of sdata type + orig_adata = extract_adata_if_sdata(data, table_key=table_key) - # get unique categories - unique_categories = np.unique(adata.obs[groups].values) + if inplace: + adata = orig_adata + else: + adata = orig_adata.copy() - # get obs x k matrix where each column is the category of the k-th neighbor - indices_with_nan = neighbor_matrix.to_numpy() - valid_indices = neighbor_matrix.fillna(-1).astype(int).to_numpy() - cat_by_id = adata.obs[groups].values[valid_indices] - cat_by_id[indices_with_nan == -1] = np.nan - # cat_by_id = np.take(category_arr, neighbor_matrix) + if library_key is not None: + # first assert that library_key was there in adata.obs, and then, stratify the object according to that library_key and + # then re-call calculate_niche_spatialleiden for each subpart, with library_key = None and prefix with appropriate information like "lib=" + assert_key_in_adata(adata, library_key, attr="obs") + logg.info(f"Stratifying by library_key '{library_key}'") - # in obs x k matrix convert categorical values to numerical values - cat_indices = {category: index for index, category in enumerate(unique_categories)} - cat_values = np.vectorize(cat_indices.get)(cat_by_id) + # go through each library_id and process the corresponding adata subset + for itr, lib_id in enumerate(adata.obs[library_key].unique()): + logg.info(f"Processing library '{lib_id}'") - # get obx x category matrix where each column is the absolute amount of a category in the neighborhood - m, k = cat_by_id.shape - abs_freq = np.zeros((m, len(unique_categories)), dtype=int) - np.add.at(abs_freq, (np.arange(m)[:, None], cat_values), 1) + lib_indices = adata.obs[adata.obs[library_key] == lib_id].index - # normalize by n_neighbors to get relative frequency of each category - rel_freq = abs_freq / k + if len(lib_indices) == 0: + logg.warning(f"Library '{lib_id}' contains no cells, skipping") + continue - if abs_nhood: - return pd.DataFrame(abs_freq, index=adata.obs.index) - else: - return pd.DataFrame(rel_freq, index=adata.obs.index) + lib_adata = adata[lib_indices].copy() + # give prefix appropriate value so that the niche values indicate lib id. + calculate_niche_spatialleiden( + lib_adata, + resolutions, + latent_connectivities_key, + spatial_connectivities_key, + layer_ratio, + n_iterations, + use_weights, + random_state, + min_niche_size, + mask, + prefix=f"lib={lib_id}_", + library_key=None, + inplace=True, # to save memory + table_key=table_key, + ) -def _utag(adata: AnnData, normalize_adj: bool, spatial_connectivity_key: str) -> AnnData: - """ - Performs inner product of adjacency matrix and feature matrix, - such that each observation inherits features from its immediate neighbors as described in UTAG paper. - """ + # from itr==1 onwards, adata will hold the columns that are being added hence, + # added_columns will be empty. Hence only obtain added_columns when itr==0 + if itr == 0: + added_columns = list(set(lib_adata.obs.columns) - set(adata.obs.columns)) - adjacency_matrix = adata.obsp[spatial_connectivity_key] + for col in added_columns: + # ensure that adata has the columns in which we are adding the information + if col not in adata.obs: + adata.obs[col] = "not_a_niche" + adata.obs.loc[lib_indices, col] = list(lib_adata.obs[col].astype("str")) - if normalize_adj: - return normalize(adjacency_matrix, norm="l1", axis=1) @ adata.X else: - return adjacency_matrix @ adata.X - - -def _setdiag(adjacency_matrix: sps.spmatrix, value: int) -> sps.spmatrix: - """remove self-loops""" + # Simply call sl.spatialleiden with the provided arguments + if not isinstance(resolutions, list): + resolutions = [resolutions] + + for res in resolutions: + sl.spatialleiden( + adata, + resolution=res, + use_weights=use_weights, + n_iterations=n_iterations, + layer_ratio=layer_ratio, + latent_neighbors_key=latent_connectivities_key, + spatial_neighbors_key=spatial_connectivities_key, + random_state=random_state, + directed=False, + key_added=f"spatialleiden_res={res}", + ) - if issparse(adjacency_matrix): - adjacency_matrix = adjacency_matrix.tolil() - adjacency_matrix.setdiag(value) - adjacency_matrix = adjacency_matrix.tocsr() - if value == 0: - adjacency_matrix.eliminate_zeros() - return adjacency_matrix + # obtain the result_columns, which are basically the difference in columns in orig_adata and adata + result_columns = [f"spatialleiden_res={res}" for res in resolutions] + # generate the list of postprocessor objects using the supplied args + postprocessors_list = [] + if mask is not None: + mask_postprocessor = MaskPostprocessor(mask) + postprocessors_list.append(mask_postprocessor) + if min_niche_size is not None: + min_niche_size_postprocessor = MinNicheSizePostprocessor(min_niche_size) + postprocessors_list.append(min_niche_size_postprocessor) + if prefix is not None: + renaming_postprocessor = RenamePostprocessor(prefix) + postprocessors_list.append(renaming_postprocessor) -def _hop( - adj_hop: sps.spmatrix, adj: sps.spmatrix, adj_visited: sps.spmatrix = None -) -> tuple[sps.spmatrix, sps.spmatrix]: - """get nearest neighbor of neighbors""" + _postprocess_niche_results(adata, result_columns, postprocessors_list) - adj_hop = adj_hop @ adj + # For SpatialData, the column names shouldn't have = sign. Hence, run sanitize_table. + # TODO: In future, change the naming standard of any niche columns added to not have '=' to be compatible with spatialdata naming + if isinstance(data, SpatialData): + sanitize_table(adata) - if adj_visited is not None: - adj_hop = adj_hop > adj_visited - adj_visited = adj_visited + adj_hop + if inplace: + return None + else: + return adata - return adj_hop, adj_visited +@d.dedent +def calculate_niche_custom( + data: AnnData | SpatialData, + embedder: NicheEmbedder, + clusterer: NicheClusterer, + postprocessors_list: list[NichePostprocessor], + library_key: str | None, + inplace: bool, + table_key: str | None, +) -> AnnData | None: + """Compute niche assignments using user-defined embedding, clustering, and postprocessing. -def _normalize(adj: sps.spmatrix) -> sps.spmatrix: - """normalize adjacency matrix such that nodes with high degree don't disproportionately affect aggregation""" + This function provides a flexible pipeline where embedding, clustering, + and postprocessing are decoupled and customizable. - deg = np.array(np.sum(adj, axis=1)).squeeze() - with np.errstate(divide="ignore"): - deg_inv = 1 / deg - deg_inv[deg_inv == float("inf")] = 0 + Parameters + ---------- + %(adata)s + embedder + Instance of :class:`NicheEmbedder` used to compute an embedding from ``adata``. + clusterer + Instance of :class:`NicheClusterer` used to assign niches based on the embedding. + postprocessors_list + List of :class:`NichePostprocessor` objects applied sequentially to refine results. + %(library_key)s + inplace + Whether to modify ``adata`` in place. + %(table_key)s - return spdiags(deg_inv, 0, len(deg_inv), len(deg_inv)) * adj + Returns + ------- + If ``inplace = True``, modifies ``adata`` in place and returns ``None``. + Otherwise, returns a copy of ``adata`` with niche annotations added to ``.obs``. + + Notes + ----- + If ``library_key`` is provided, the computation is performed independently + for each library and results are merged back into ``adata``. + + See Also + -------- + calculate_niche_neighborhood : Convenience wrapper for neighborhood flavor niche analysis. + calculate_niche_utag : Convenience wrapper for utag flavor niche analysis. + calculate_niche_cellcharter : Convenience wrapper for cellcharter flavor niche analysis. + calculate_niche_spatialleiden : Convenience wrapper for spatialleiden flavor niche analysis. + NicheEmbedder : Base class for embedding strategies. + NicheClusterer : Base class for clustering strategies. + NichePostprocessor : Base class for postprocessing steps. + """ + # obtain adata if data was of sdata type + orig_adata = extract_adata_if_sdata(data, table_key=table_key) -def _aggregate(adata: AnnData, normalized_adjacency_matrix: sps.spmatrix, aggregation: str = "mean") -> Any: - """aggregate count and adjacency matrix either by mean or variance""" - # TODO: add support for other aggregation methods - if aggregation == "mean": - aggregated_matrix = normalized_adjacency_matrix @ adata.X - elif aggregation == "variance": - mean_matrix = (normalized_adjacency_matrix @ adata.X).toarray() - X_to_arr = adata.X.toarray() - mean_squared_matrix = normalized_adjacency_matrix @ (X_to_arr * X_to_arr) - aggregated_matrix = mean_squared_matrix - mean_matrix * mean_matrix + if inplace: + adata = orig_adata else: - raise ValueError(f"Invalid aggregation method '{aggregation}'. Please choose either 'mean' or 'variance'.") - - return aggregated_matrix + adata = orig_adata.copy() + if library_key is not None: + assert_key_in_adata(adata, library_key, attr="obs") + logg.info(f"Stratifying by library_key '{library_key}'") -def _get_GMM_clusters(A: NDArray[np.float64], n_components: int, random_state: int) -> Any: - """Returns niche labels generated by GMM clustering. - Compared to cellcharter this approach is simplified by using sklearn's GaussianMixture model without stability analysis. - """ + # go through each library_id and process the corresponding adata subset + for itr, lib_id in enumerate(adata.obs[library_key].unique()): + logg.info(f"Processing library '{lib_id}'") - gmm = GaussianMixture( - n_components=n_components, - random_state=random_state, - init_params="random_from_data", - ) - gmm.fit(A) - labels = gmm.predict(A) + lib_indices = adata.obs[adata.obs[library_key] == lib_id].index - return labels + if len(lib_indices) == 0: + logg.warning(f"Library '{lib_id}' contains no cells, skipping") + continue + lib_adata = adata[lib_indices].copy() -def _get_spatialleiden_domains( - adata: AnnData, - spatial_connectivities_key: str, - latent_connectivities_key: str, - resolutions: float | tuple[float, float] | list[float | tuple[float, float]], - layer_ratio: float, - use_weights: bool | tuple[bool, bool], - n_iterations: int, - random_state: int, -) -> None: - """ - Perform SpatialLeiden clustering. + # append a renaming postprocessor to postprocessors_list_lib + renaming_postprocessor = RenamePostprocessor(prefix_for_niches=f"lib={lib_id}_") + postprocessors_list_lib = postprocessors_list + [renaming_postprocessor] - This is a wrapper around :py:func:`spatialleiden.multiplex_leiden` that uses :py:class:`anndata.AnnData` as input and works with two layers; one latent space and one spatial layer. + calculate_niche_custom( + lib_adata, + embedder, + clusterer, + postprocessors_list_lib, + library_key=None, + inplace=True, # to save memory + table_key=None, + ) - Adapted from https://github.com/HiDiHlabs/SpatialLeiden/. - """ - try: - import spatialleiden as sl - except ImportError as e: - msg = "Please install the spatialleiden algorithm: `pip install squidpy[leiden]` or `conda install bioconda::spatialleiden` or `pip install spatialleiden`." - raise ImportError(msg) from e + # from itr==1 onwards, adata will hold the columns that are being added hence, + # added_columns will be empty. Hence only obtain added_columns when itr==0 + if itr == 0: + added_columns = list(set(lib_adata.obs.columns) - set(adata.obs.columns)) - if not isinstance(resolutions, list): - resolutions = [resolutions] - - for res in resolutions: - sl.spatialleiden( - adata, - resolution=res, - use_weights=use_weights, - n_iterations=n_iterations, - layer_ratio=layer_ratio, - latent_neighbors_key=latent_connectivities_key, - spatial_neighbors_key=spatial_connectivities_key, - random_state=random_state, - directed=False, - key_added=f"spatialleiden_res={res}", - ) + for col in added_columns: + # ensure that adata has the columns in which we are adding the information + if col not in adata.obs: + adata.obs[col] = "not_a_niche" + adata.obs.loc[lib_indices, col] = list(lib_adata.obs[col].astype("str")) - return + else: + # supply the adata object to the embedder object, and obtain appropriate embedding matrix + embedding = embedder.get_embedding(adata) + # Supply to the clusterer object, the embedding matrix just obtained, and get the appropriate clustering. + result_columns = clusterer.cluster(adata, embedding) -def _fide_score(adata: AnnData, niche_key: str, average: bool) -> Any: - """ - F1-score of intra-domain edges (FIDE). A high score indicates a great domain continuity. + # do postprocessing + _postprocess_niche_results(adata, result_columns, postprocessors_list) - The F1-score is computed for every class, then all F1-scores are averaged. If some classes - are not predicted, the `n_classes` argument allows to pad with zeros before averaging the F1-scores. - """ - i, j = adata.obsp["spatial_connectivities"].nonzero() # get row and column indices of non-zero elements - niche_labels, neighbor_niche_labels = ( - adata.obs.iloc[i][niche_key], - adata.obs.iloc[j][niche_key], - ) + # For SpatialData, the column names shouldn't have = sign. Hence, run sanitize_table. + # TODO: In future, change the naming standard of any niche columns added to not have '=' to be compatible with spatialdata naming + if isinstance(data, SpatialData): + sanitize_table(adata) - if not average: - fide = f1_score(niche_labels, neighbor_niche_labels, average=None) + if inplace: + return None else: - fide = f1_score(niche_labels, neighbor_niche_labels, average="macro") + return adata - return fide +def _postprocess_niche_results( + adata: AnnData, + result_columns: list[str], + postprocessors_list: list[NichePostprocessor], +) -> None: + """Apply a sequence of postprocessors to niche assignment results. -def _jensen_shannon_divergence(adata: AnnData, niche_key: str, library_key: str) -> Any: - """ - Calculate Jensen-Shannon divergence (JSD) over all slides. - This metric measures how well niche label distributions match across different slides. + Parameters + ---------- + adata + Annotated data matrix. + result_columns + List of column names in ``adata.obs`` containing initial niche assignments. + postprocessors_list + List of :class:`NichePostprocessor` objects applied sequentially. + + Notes + ----- + Each postprocessor may create new columns in ``adata.obs`` and returns + the updated list of result column names, which are passed to the next step. """ - niche_labels = sorted(adata.obs[niche_key].unique()) - label_distributions = [] - - for _, slide in adata.obs.groupby(library_key): - counts = slide[niche_key].value_counts(normalize=True) - relative_freq = [counts.get(label, 0) for label in niche_labels] - label_distributions.append(relative_freq) + # go through each postprocessor object, and apply it to the adata, and store + # results in the form of new columns in adata + for postprocessor in postprocessors_list: + # obtain the new columns in this process + result_columns = postprocessor.postprocess(adata, result_columns) - return distance.jensenshannon(np.array(label_distributions)) + return def _validate_niche_args( @@ -878,7 +903,7 @@ def _validate_niche_args( # Define parameters used by each flavor flavor_param_specs = { "neighborhood": { - "required": ["groups", "n_neighbors", "resolutions"], + "required": ["groups", "n_neighbors", "resolutions", "spatial_connectivities_key"], "optional": [ "min_niche_size", "scale", @@ -886,10 +911,19 @@ def _validate_niche_args( "distance", "n_hop_weights", ], - "unused": ["aggregation", "n_components", "random_state"], + "unused": [ + "aggregation", + "n_components", + "random_state", + "latent_connectivities_key", + "layer_ratio", + "n_iterations", + "use_weights", + "use_rep", + ], }, "utag": { - "required": ["n_neighbors", "resolutions"], + "required": ["n_neighbors", "resolutions", "spatial_connectivities_key"], "optional": [], "unused": [ "groups", @@ -901,10 +935,15 @@ def _validate_niche_args( "aggregation", "n_components", "random_state", + "latent_connectivities_key", + "layer_ratio", + "n_iterations", + "use_weights", + "use_rep", ], }, "cellcharter": { - "required": ["distance", "aggregation", "random_state"], + "required": ["distance", "aggregation", "random_state", "spatial_connectivities_key"], "optional": ["n_components", "use_rep"], "unused": [ "groups", @@ -914,6 +953,10 @@ def _validate_niche_args( "n_neighbors", "resolutions", "n_hop_weights", + "latent_connectivities_key", + "layer_ratio", + "n_iterations", + "use_weights", ], }, "spatialleiden": { @@ -925,14 +968,7 @@ def _validate_niche_args( "use_weights", "random_state", ], - "unused": [ - "groups", - "min_niche_size", - "scale", - "abs_nhood", - "n_neighbors", - "n_hop_weights", - ], + "unused": ["groups", "min_niche_size", "scale", "abs_nhood", "n_neighbors", "n_hop_weights", "use_rep"], }, } @@ -1049,3 +1085,647 @@ def _check_unnecessary_args(flavor: str, param_dict: dict[str, Any], param_specs logg.warning( f"Parameters {', '.join([f'{arg}' for arg in unnecessary_args])} are not used for flavor '{flavor}'.", ) + + +############ +### embedder classes +############ + + +class NicheEmbedder: + """Base class for computing embeddings used in niche analysis. + + Subclasses must implement :meth:`get_embedding`, which transforms an + :class:`AnnData` object into a feature matrix suitable for clustering. + The 0-index dimension of returned embedding (embedding.shape[0]) + should correspond to the rows in adata.obs (and adata.X), meaning in + the same order and having same length. + """ + + @abstractmethod + def get_embedding(self, adata: AnnData) -> NDArrayA: + """return an embedding matrix, with cells as rows""" + + +@d.dedent +class NhoodProfileEmbedder(NicheEmbedder): + """Compute neighborhood composition profiles as embeddings. + + Each observation is represented by the frequency of categorical labels + (e.g. cell types) in its spatial neighborhood. Optionally, higher-order + neighborhoods (multi-hop) can be incorporated. + + Parameters + ---------- + groups + Column in ``adata.obs`` defining categorical labels. + spatial_connectivities_key + Key in ``adata.obsp`` containing the spatial connectivity matrix. + scale + Whether to z-score the resulting embedding. + distance + Number of hops to consider for neighborhood aggregation. + abs_nhood + If ``True``, use absolute counts of categories in neighborhood; otherwise + normalize to proportions. + n_hop_weights + Weights for combining profiles across neighborhood hops. + + Notes + ----- + For ``distance > 1``, neighborhood profiles are iteratively aggregated using + powers of the adjacency matrix, optionally weighted per hop. + """ + + def __init__( + self, + groups: str, + spatial_connectivities_key: str, + scale: bool, + distance: int, + abs_nhood: bool, + n_hop_weights: list[float] | None, + ): + super().__init__() + self.groups = groups + self.spatial_connectivities_key = spatial_connectivities_key + self.scale = scale + self.distance = distance + self.abs_nhood = abs_nhood + self.n_hop_weights = n_hop_weights + + def _calculate_neighborhood_profile( + self, + adata: AnnData, + matrix: coo_matrix, + ) -> pd.DataFrame: + """ + Returns an obs x category matrix where each column is the absolute/relative frequency of a category in the neighborhood + """ + + # ensure that adata.obs[group] is of categorical type, as that makes it explicit, which cols of the returned profile_df + # correspond to which categories in group + if adata.obs[self.groups].dtype.name != "category": + warnings.warn( + "Since adata.obs[groups] does not already have categorical dtype, converting it into categorical type.", + stacklevel=2, + ) + adata.obs[self.groups] = adata.obs[self.groups].astype("category") + + # ensure matrix is in csc format for efficient column slicing + if matrix.format != "csc": + matrix = matrix.tocsc() + + # get cell categories in order + categories_order = adata.obs[self.groups].cat.categories + n_categories = len(categories_order) + + # map category to column index + category_to_idx = {ct: i for i, ct in enumerate(categories_order)} + + # pre allocate sparse LIL matrix for efficient assignment (n_cells x n_categories) + profile_sparse = lil_matrix((matrix.shape[0], n_categories), dtype=np.float64) + + # for each category, sum over cells of that category + for ct in categories_order: + ct_mask = adata.obs[self.groups] == ct # boolean mask for cells of this category + col_indices = np.where(ct_mask)[0] # indices of those cells + if len(col_indices) > 0: + col_slice = matrix[:, col_indices] # sparse submatrix + profile_sparse[:, category_to_idx[ct]] = col_slice.sum(axis=1).A1 + + # convert to dataframe (csr for final storage, dense for pandas) + profile_df = pd.DataFrame( + profile_sparse.tocsr().todense(), index=adata.obs[self.groups].index, columns=categories_order + ) + + # now according to parameter abs_nhood, make raw counts into proportions or not + if not self.abs_nhood: + total_neighs = profile_df.sum(axis=1) + profile_df = profile_df.div(total_neighs, axis=0) + # this may lead to some values being nan, as some cells might have had no neighbors. Make those values as 0 + profile_df = profile_df.fillna(0.0) + + return profile_df + + def get_embedding(self, adata: AnnData) -> NDArrayA: + """ + adapted from https://github.com/immunitastx/monkeybread/blob/main/src/monkeybread/calc/_neighborhood_profile.py + """ + + # get obs x neighbor matrix from sparse matrix + matrix = adata.obsp[self.spatial_connectivities_key].tocoo() + + # get obs x category matrix where each column is the absolute/relative frequency of a category in the neighborhood + nhood_profile = self._calculate_neighborhood_profile(adata, matrix) + + # Additionally use n-hop neighbors if distance > 1. This sums up the (weighted) neighborhood profiles of all n-hop neighbors. + if self.distance > 1: + n_hop_adjacency_matrix = adata.obsp[self.spatial_connectivities_key].copy() + # if no weights are provided, use 1 for all n_hop neighbors + if self.n_hop_weights is None: + self.n_hop_weights = [1] * self.distance + # if weights are provided, start with applying weight to the original neighborhood profile + elif len(self.n_hop_weights) < self.distance: + # Extend weights if too few provided + self.n_hop_weights = self.n_hop_weights + [self.n_hop_weights[-1]] * ( + self.distance - len(self.n_hop_weights) + ) + logg.debug(f"Extended weights to match distance: {self.n_hop_weights}") + + # Apply first weight to base profile + weighted_profile = self.n_hop_weights[0] * nhood_profile + + # Calculate higher-order hop profiles + n_hop_adjacency_matrix = adata.obsp[self.spatial_connectivities_key].copy() + + # get n_hop neighbor adjacency matrices by multiplying the original adjacency matrix with itself n times and get corresponding neighborhood profiles. + for n_hop in range(1, self.distance): + logg.debug(f"Calculating {n_hop + 1}-hop neighbors") + # Multiply adjacency matrix by itself to get n+1 hop adjacency + n_hop_adjacency_matrix = n_hop_adjacency_matrix @ adata.obsp[self.spatial_connectivities_key] + matrix = n_hop_adjacency_matrix.tocoo() + + # Calculate and add weighted profile + hop_profile = self._calculate_neighborhood_profile(adata, matrix) + weighted_profile += self.n_hop_weights[n_hop] * hop_profile + + if not self.abs_nhood: + weighted_profile = weighted_profile / sum(self.n_hop_weights) + + nhood_profile = weighted_profile + + # create AnnData object from neighborhood profile to perform scanpy functions + # Use .to_numpy(copy=True) to ensure the array is writeable (required for pandas CoW compatibility) + # Preserve the DataFrame index for later matching with adata_masked + adata_neighborhood = ad.AnnData( + X=nhood_profile.to_numpy(copy=True), obs=pd.DataFrame(index=nhood_profile.index) + ) + + # reason for scaling see https://monkeybread.readthedocs.io/en/latest/notebooks/tutorial.html#niche-analysis + if self.scale: + sc.pp.scale(adata_neighborhood, zero_center=True) + return adata_neighborhood.X + + +@d.dedent +class UtagEmbedder(NicheEmbedder): + """Compute a UTAG-style embedding by propagating features over spatial neighbors. + + The embedding is constructed by normalizing the spatial connectivity matrix, + multiplying it by ``adata.X``, and then applying PCA to the propagated + feature matrix. + + Parameters + ---------- + spatial_connectivities_key + Key in ``adata.obsp`` containing the spatial connectivity matrix. + + Notes + ----- + This follows the general UTAG idea that each observation inherits information + from its immediate spatial neighborhood before dimensionality reduction. + """ + + def __init__( + self, + spatial_connectivities_key: str, + ): + super().__init__() + self.spatial_connectivities_key = spatial_connectivities_key + + def get_embedding(self, adata: AnnData) -> NDArrayA: + """ + Performs inner product of adjacency matrix and feature matrix, + such that each observation inherits features from its immediate neighbors as described in UTAG paper. + """ + + adjacency_matrix = adata.obsp[self.spatial_connectivities_key] + new_feature_matrix = normalize(adjacency_matrix, norm="l1", axis=1) @ adata.X + adata_utag = ad.AnnData(X=new_feature_matrix) + sc.tl.pca(adata_utag) # note: unlike with flavor 'neighborhood' dim reduction is performed here + return adata_utag.obsm["X_pca"] + + +# TODO: This function requires some work later on. Right now keeping the implementation just like how +# it was before the refactor, and in that case, when use_rep was provided, then it simply returned +# that as the embedding, so no cellcharter algorithm used in that case +@d.dedent +class CellcharterEmbedder(NicheEmbedder): + """Compute a CellCharter-style embedding from spatially aggregated features. + + The embedding can either be derived from a precomputed representation in + ``adata.obsm`` or constructed by aggregating features across multi-hop + spatial neighborhoods. + + Parameters + ---------- + distance + Number of neighborhood hops to aggregate. + aggregation + Aggregation strategy to apply to neighborhood features, such as + ``"mean"`` or ``"variance"``. + spatial_connectivities_key + Key in ``adata.obsp`` containing the spatial connectivity matrix. + n_components + Number of components to keep from the input representation when ``use_rep`` + is provided. + use_rep + Key in ``adata.obsm`` pointing to the representation to use. If ``None``, + a spatially aggregated embedding is constructed from ``adata.X``. + + Notes + ----- + When ``use_rep`` is ``None``, PCA is applied to the concatenated aggregated + feature matrix to produce the final embedding. + """ + + def __init__( + self, + distance: int | None, + aggregation: str | None, + spatial_connectivities_key: str | None, + n_components: int | None, + use_rep: str | None, + ): + super().__init__() + self.distance = distance + self.aggregation = aggregation + self.spatial_connectivities_key = spatial_connectivities_key + self.n_components = n_components + self.use_rep = use_rep + + def _setdiag(self, adjacency_matrix: sps.spmatrix, value: int) -> sps.spmatrix: + """remove self-loops""" + + if issparse(adjacency_matrix): + adjacency_matrix = adjacency_matrix.tolil() + adjacency_matrix.setdiag(value) + adjacency_matrix = adjacency_matrix.tocsr() + if value == 0: + adjacency_matrix.eliminate_zeros() + return adjacency_matrix + + def _hop( + self, + adj_hop: sps.spmatrix, + adj: sps.spmatrix, + adj_visited: sps.spmatrix = None, + ) -> tuple[sps.spmatrix, sps.spmatrix]: + """get nearest neighbor of neighbors""" + + adj_hop = adj_hop @ adj + + if adj_visited is not None: + adj_hop = adj_hop > adj_visited + adj_visited = adj_visited + adj_hop + + return adj_hop, adj_visited + + def _normalize(self, adj: sps.spmatrix) -> sps.spmatrix: + """normalize adjacency matrix such that nodes with high degree don't disproportionately affect aggregation""" + + deg = np.array(np.sum(adj, axis=1)).squeeze() + with np.errstate(divide="ignore"): + deg_inv = 1 / deg + deg_inv[deg_inv == float("inf")] = 0 + + return spdiags(deg_inv, 0, len(deg_inv), len(deg_inv)) * adj + + def _aggregate(self, adata: AnnData, normalized_adjacency_matrix: sps.spmatrix, aggregation: str = "mean") -> Any: + """aggregate count and adjacency matrix either by mean or variance""" + # TODO: add support for other aggregation methods + if aggregation == "mean": + aggregated_matrix = normalized_adjacency_matrix @ adata.X + elif aggregation == "variance": + mean_matrix = (normalized_adjacency_matrix @ adata.X).toarray() + X_to_arr = adata.X.toarray() + mean_squared_matrix = normalized_adjacency_matrix @ (X_to_arr * X_to_arr) + aggregated_matrix = mean_squared_matrix - mean_matrix * mean_matrix + else: + raise ValueError(f"Invalid aggregation method '{aggregation}'. Please choose either 'mean' or 'variance'.") + + return aggregated_matrix + + # this will hold an if block checking if use_rep is not None. If not None, then it will simply + # return that representation from adata + # Also a note for user, n_components is only used when use_rep is not None. It is the number of + # components from that representation to use as the embedding + # aggregation is only used when use_rep is None + def get_embedding(self, adata: AnnData) -> NDArrayA: + """adapted from https://github.com/CSOgroup/cellcharter/blob/main/src/cellcharter/gr/_aggr.py + and https://github.com/CSOgroup/cellcharter/blob/main/src/cellcharter/tl/_gmm.py""" + + if self.use_rep is not None: + # Use provided embedding from adata.obsm + assert_key_in_adata(adata, self.use_rep, attr="obsm") + embedding = adata.obsm[self.use_rep] + # Ensure embedding has the right number of components + if embedding.shape[1] < self.n_components: + raise ValueError( + f"Embedding has {embedding.shape[1]} components, but n_components={self.n_components}. Please provide an embedding with at least {self.n_components} components." + ) + # Use only the first n_components + embedding = embedding[:, : self.n_components] + else: + logg.warning( + "CellCharter recommends to use a dimensionality reduced embedding of the data, e.g. a scVI embedding. Since 'use_rep' is not provided, PCA will be used as proxy - performance may be suboptimal." + ) + adjacency_matrix = adata.obsp[self.spatial_connectivities_key] + layers = list(range(self.distance + 1)) + + aggregated_matrices = [] + adj_hop = self._setdiag(adjacency_matrix, 0) # Remove self-loops, set diagonal to 0 + adj_visited = self._setdiag(adjacency_matrix.copy(), 1) # Track visited neighbors + for k in layers: + if k == 0: + # get original count matrix (not aggregated) + aggregated_matrices.append(adata.X) + else: + # get count and adjacency matrix for k-hop (neighbor of neighbor of neighbor ...) and aggregate them + if k > 1: + adj_hop, adj_visited = self._hop(adj_hop, adjacency_matrix, adj_visited) + adj_hop_norm = self._normalize(adj_hop) + aggregated_matrix = self._aggregate(adata, adj_hop_norm, self.aggregation) + aggregated_matrices.append(aggregated_matrix) + + concatenated_matrix = hstack(aggregated_matrices) # Stack all matrices horizontally + arr = concatenated_matrix.toarray() # Densify + + arr_ad = ad.AnnData(X=arr) + sc.tl.pca(arr_ad) + embedding = arr_ad.obsm["X_pca"] + + return embedding + + +############ +### clusterer classes +############ + + +class NicheClusterer: + """Base class for clustering embeddings into niche assignments. + + Subclasses must implement :meth:`cluster`, which assigns cluster labels + and stores them in ``adata.obs``. + """ + + @abstractmethod + def cluster(self, adata: AnnData, embedding: NDArrayA) -> list[str]: + """Adds column/s in adata.obs with the clustering done. Returns the names of the columns just added.""" + + +@d.dedent +class LeidenClusterer(NicheClusterer): + """Cluster embeddings using the Leiden algorithm. + + Parameters + ---------- + n_neighbors + Number of neighbors used to construct the kNN graph. + resolutions + Resolution parameter(s) for Leiden clustering. Can be a single + float value or list of floats. + base_colname + Base name for columns added to ``adata.obs``. Resolution is + appended to this to unique identify columns for each resolution. + + Notes + ----- + A separate clustering is computed for each resolution, producing multiple + niche annotation columns. + """ + + def __init__( + self, + n_neighbors: int, + resolutions: float | list[float], + base_colname: str = "niche_leiden", + ): + super().__init__() + self.n_neighbors = n_neighbors + self.resolutions = resolutions if isinstance(resolutions, list) else [resolutions] + self.base_colname = base_colname + + def cluster(self, adata: AnnData, embedding: NDArrayA) -> list: + # first create an adata object using the embedding provided + adata_embedding = ad.AnnData(X=embedding, obs=pd.DataFrame(index=adata.obs.index)) + + # required for leiden clustering (note: no dim reduction performed in original implementation) + sc.pp.neighbors(adata_embedding, n_neighbors=self.n_neighbors, use_rep="X") + + # For each resolution, apply leiden on neighborhood profile. Each cluster label equals to a niche label + niche_keys = [] + for res in self.resolutions: + 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, + ) + + adata.obs[niche_key] = list( + adata_embedding.obs[niche_key] + ) # since constrain all embedders to return embedding with numrows==numcells and in same order, this should be fine + + return niche_keys + + +@d.dedent +class GMMClusterer(NicheClusterer): + """Cluster embeddings with a Gaussian mixture model. + + Parameters + ---------- + n_components + Number of mixture components. + random_state + Random seed used by the Gaussian mixture model. + base_colname + Name of the output column added to ``adata.obs``. + + Notes + ----- + Cluster assignments are stored as categorical niche labels in ``adata.obs``. + """ + + def __init__( + self, + n_components: int, + random_state: int, + base_colname: str = "niche_gmm", + ): + super().__init__() + self.n_components = n_components + self.random_state = random_state + self.base_colname = base_colname + + def cluster(self, adata: AnnData, embedding: NDArrayA) -> list: + """Returns niche labels generated by GMM clustering. + Compared to cellcharter this approach is simplified by using sklearn's GaussianMixture model without stability analysis. + """ + # cluster concatenated matrix with GMM, each cluster label equals to a niche label + gmm = GaussianMixture( + n_components=self.n_components, + random_state=self.random_state, + init_params="random_from_data", + ) + gmm.fit(embedding) + niches = gmm.predict(embedding) + + if self.base_colname in adata.obs.columns: + logg.info(f"Overwriting existing column '{self.base_colname}'") + + adata.obs[self.base_colname] = pd.Categorical(niches) + return [self.base_colname] + + +############ +### postprocessor classes +############ + + +class NichePostprocessor: + """Base class for postprocessing niche assignments. + + Postprocessors operate on clustering results stored in ``adata.obs`` and + typically generate new columns derived from existing niche columns. + """ + + def __init__(self, suffix: str): + self.suffix = suffix + + @abstractmethod + def postprocess(self, adata: AnnData, result_columns: list[str]) -> list[str]: + """Logic to postprocess adata and return the names of columns added.""" + # should append add self.suffix to the columns added + + +@d.dedent +class MinNicheSizePostprocessor(NichePostprocessor): + """Filter niches below a minimum size threshold. + + Parameters + ---------- + min_niche_size + Minimum number of observations required for a niche. + suffix + Suffix appended to result column names. + + Notes + ----- + Niche labels with fewer than ``min_niche_size`` observations are replaced + with ``"not_a_niche"``. + """ + + def __init__(self, min_niche_size: int, suffix: str = "_size_filter"): + super().__init__(suffix=suffix) + self.min_niche_size = min_niche_size + + def postprocess(self, adata: AnnData, result_columns: list[str]) -> list[str]: + new_result_columns = [] + # filter niches with n_cells < min_niche_size + for result_column in result_columns: + # copy into new column + new_result_column = result_column + self.suffix + new_result_columns.append(new_result_column) + adata.obs[new_result_column] = list(adata.obs[result_column]) + + counts_by_niche = adata.obs[new_result_column].value_counts() + to_filter = counts_by_niche[counts_by_niche < self.min_niche_size].index + + if new_result_column in adata.obs.columns: + logg.info(f"Overwriting existing column '{new_result_column}'") + + adata.obs[new_result_column] = adata.obs[new_result_column].apply( + lambda x, to_filter=to_filter: "not_a_niche" if x in to_filter else x + ) + adata.obs[new_result_column] = adata.obs.index.map(adata.obs[new_result_column]).fillna("not_a_niche") + + return new_result_columns + + +@d.dedent +class MaskPostprocessor(NichePostprocessor): + """Mask selected observations from niche assignments. + + Parameters + ---------- + mask + Boolean mask or index specifying observations to exclude. + suffix + Suffix appended to result column names. + + Notes + ----- + Observations included in ``mask`` are assigned the label ``"not_a_niche"``. + + Examples + ----- + Mask can look like the following. Here, the index values would correspond to adata.obs.index. + The entries that are False are the ones ignored. + mask = Series([False, False, True], index = ["a", "b", "c"]) + """ + + def __init__(self, mask: pd.Series, suffix: str = "_mask"): + super().__init__(suffix=suffix) + self.mask = mask + + def postprocess(self, adata: AnnData, result_columns: list[str]) -> list[str]: + new_result_columns = [] + # mask obs to exclude cells for which no niche shall be assigned + for result_column in result_columns: + # copy into new column + new_result_column = result_column + self.suffix + new_result_columns.append(new_result_column) + adata.obs[new_result_column] = list(adata.obs[result_column]) + + if new_result_column in adata.obs.columns: + logg.info(f"Overwriting existing column '{new_result_column}'") + + to_filter = self.mask[self.mask.index.isin(adata.obs.index)] + adata.obs.loc[~to_filter, new_result_column] = "not_a_niche" + + return new_result_columns + + +@d.dedent +class RenamePostprocessor(NichePostprocessor): + """Rename niche labels by adding a prefix. + + Parameters + ---------- + prefix_for_niches + Prefix added to each niche label. + suffix + Suffix appended to result column names. + + Notes + ----- + This is useful when combining results across subsets (e.g. libraries) + to ensure unique niche identifiers. + """ + + def __init__(self, prefix_for_niches: str, suffix: str = "_renamed"): + super().__init__(suffix=suffix) + self.prefix_for_niches = prefix_for_niches + + def postprocess(self, adata: AnnData, result_columns: list[str]) -> list[str]: + new_result_columns = [] + for result_column in result_columns: + # copy into new column + new_result_column = result_column + self.suffix + new_result_columns.append(new_result_column) + + if new_result_column in adata.obs.columns: + logg.info(f"Overwriting existing column '{new_result_column}'") + + adata.obs[new_result_column] = self.prefix_for_niches + adata.obs[result_column].astype(str) + + return new_result_columns diff --git a/tests/conftest.py b/tests/conftest.py index 93f32958f..cd26a2ef5 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -117,6 +117,63 @@ def dummy_adata() -> AnnData: return adata +@pytest.fixture() +def dummy_adata2() -> AnnData: + r = np.random.RandomState(100) + adata = AnnData(r.rand(10, 100), obs={"celltype": r.choice(["foo", "bar", "baz"], size=10)}) + adata.obs.index = ["a", "b", "c", "d", "e", "f", "g", "h", "i", "j"] + + # Spatial layout of the 10 data points in the grid: + # + # Y ▲ + # │ + # 5 │ a b c Points: (2,5), (3,5), (5,5) + # │ + # 4 │ d e Points: (1,4), (4,4) + # │ + # 3 │ f g Points: (3,3), (5,3) + # │ + # 2 │ h Points: (4,2) + # │ + # 1 │ i j Points: (1,1), (5,1) + # │ + # └─────────────────────────► + # 1 2 3 4 5 X + # + # + # celltypes- + # a bar + # b baz + # c foo + # d foo + # e foo + # f baz + # g baz + # h baz + # i bar + # j baz + + adata.obsm["spatial"] = np.array( + [ + [2, 5], + [3, 5], + [5, 5], + [1, 4], + [4, 4], + [3, 3], + [5, 3], + [4, 2], + [1, 1], + [5, 1], + ] + ) + + # using a radius of 1.5 will lead to, say e being a neighbor of b,c,f,g but not h. So all side-adjacent and diagonal-adjacent + # locations will be neighbors + sq.gr.spatial_neighbors_radius(adata, radius=1.5) + return adata + + @pytest.fixture() def adata_intmat() -> AnnData: graph = csr_matrix( diff --git a/tests/graph/test_niche.py b/tests/graph/test_niche.py index 5ba93ed50..f86e191f0 100644 --- a/tests/graph/test_niche.py +++ b/tests/graph/test_niche.py @@ -1,16 +1,283 @@ from __future__ import annotations from anndata import AnnData -from pandas.testing import assert_frame_equal -from scipy.sparse import issparse +from pandas import Categorical, DataFrame, Series +from scanpy.pp import neighbors +from scipy.sparse import csr_matrix +from spatialdata import SpatialData +from spatialdata.models import TableModel from squidpy.gr import calculate_niche, spatial_neighbors_knn -from squidpy.gr._niche import _calculate_neighborhood_profile, _utag +from squidpy.gr._niche import NhoodProfileEmbedder, UtagEmbedder SPATIAL_CONNECTIVITIES_KEY = "spatial_connectivities" N_NEIGHBORS = 20 GROUPS = "celltype_mapped_refined" +# test if calculate_niche() gives appropriate output for dummy_adata2 for the different flavors + + +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() + + +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() + + +def test_niche_calc_cellcharter_dummy_adata(dummy_adata2: AnnData): + "Check whether niche calculation using cellcharter approach works as intended for dummy_adata2." + + # since cellcharter throws an error if the object's expression matrix is not sparse, first ensure that is the case + dummy_adata2.X = csr_matrix(dummy_adata2.X) + + calculate_niche(dummy_adata2, flavor="cellcharter", distance=2, aggregation="mean", random_state=0) + + assert "cellcharter_niche" in dummy_adata2.obs.columns + + expected_niches = Series( + Categorical([8, 4, 0, 7, 2, 9, 5, 6, 1, 3], 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() + + +def test_niche_calc_spatialleiden_dummy_adata(dummy_adata2: AnnData): + "Check whether niche calculation using spatialleiden approach works as intended for dummy_adata2." + + # need the latent_connectivities_key, meaning have to run the graph construction + neighbors(dummy_adata2, n_neighbors=3, use_rep="X") + + calculate_niche( + dummy_adata2, + flavor="spatialleiden", + latent_connectivities_key="connectivities", + spatial_connectivities_key="spatial_connectivities", + resolutions=1.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() + + +# more special test cases + + +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", + ] + + calculate_niche( + dummy_adata2, flavor="neighborhood", groups="celltype", n_neighbors=3, resolutions=1.5, library_key="batch" + ) + + assert "nhood_niche_res=1.5_renamed" 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_renamed"]).all() + + +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." + + # 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", + ] + + calculate_niche( + dummy_adata2, + flavor="spatialleiden", + latent_connectivities_key="connectivities", + spatial_connectivities_key="spatial_connectivities", + resolutions=1.0, + library_key="batch", + ) + + assert "spatialleiden_res=1.0_renamed" 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_renamed"]).all() + + +def test_niche_calc_nhood_multipostprocessor_dummy_adata(dummy_adata2: AnnData): + "Check whether niche calculation using neighborhood profile approach works as intended for dummy_adata2, when using both, mask and min_niche_size postprocessors" + mask = Series( + [False, False, True, True, True, True, True, True, True, True], + index=["a", "b", "c", "d", "e", "f", "g", "h", "i", "j"], + ) + calculate_niche( + dummy_adata2, + flavor="neighborhood", + groups="celltype", + n_neighbors=3, + resolutions=1.0, + mask=mask, + min_niche_size=3, + ) + assert "nhood_niche_res=1.0_mask_size_filter" 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_mask_size_filter"]).all() + + +def test_niche_calc_nhood_dummy_sdata(dummy_adata2: AnnData): + "Check whether niche calculation works as intended for the spatialdata version of dummy_adata2." + + # make adata into sdata object + adata_for_sdata = TableModel.parse(dummy_adata2) + sdata = SpatialData( + # images={"hne": img_for_sdata}, + # shapes={"spots": shapes_for_sdata}, + tables={"adata": adata_for_sdata}, + ) + + 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() + + +# older tests + + +def test_calculate_neighborhood_profile(dummy_adata2: AnnData): + "calculate_neighborhood_profile function needs to be tested, as it is at the base of the functionality of neighborhood flavor" + + # since now this has become a method of an embedder class, will initialize the class and call the method for that class + embedder1 = NhoodProfileEmbedder( + groups="celltype", + spatial_connectivities_key="spatial_connectivities", + scale=True, # doesnt make a difference in this test function + distance=1, # doesnt make a difference in this test function + abs_nhood=True, + n_hop_weights=None, # doesnt make a difference in this test function + ) + + embedder2 = NhoodProfileEmbedder( + groups="celltype", + spatial_connectivities_key="spatial_connectivities", + scale=True, # doesnt make a difference in this test function + distance=1, # doesnt make a difference in this test function + abs_nhood=False, + n_hop_weights=None, # doesnt make a difference in this test function + ) + + matrix = dummy_adata2.obsp["spatial_connectivities"].tocoo() + nhood_profile = embedder1._calculate_neighborhood_profile(dummy_adata2, matrix) + relative_nhood_profile = embedder2._calculate_neighborhood_profile(dummy_adata2, matrix) + + # nhood_profile and relative_nhood_profile should have same entries as this manually determined version for dummy_adata2 + expected_nhood_profile = DataFrame( + { + 0: [0, 1, 0, 1, 0, 0, 0, 0, 0, 0], + 1: [1, 0, 0, 0, 3, 1, 1, 3, 0, 1], + 2: [1, 1, 1, 0, 1, 1, 1, 0, 0, 0], + }, + index=list("abcdefghij"), + dtype="float", + ) + total_neighs = expected_nhood_profile.sum(axis=1) + expected_relative_nhood_profile = expected_nhood_profile.div(total_neighs, axis=0) + expected_relative_nhood_profile = expected_relative_nhood_profile.fillna(0.0) + + # compare + assert (nhood_profile.values == expected_nhood_profile.values).all() + assert (relative_nhood_profile.values == expected_relative_nhood_profile.values).all() + def test_niche_calc_nhood(adata_seqfish: AnnData): """Check whether niche calculation using neighborhood profile approach works as intended.""" @@ -23,7 +290,7 @@ def test_niche_calc_nhood(adata_seqfish: AnnData): resolutions=[0.1], min_niche_size=100, ) - niches = adata_seqfish.obs["nhood_niche_res=0.1"] + niches = adata_seqfish.obs["nhood_niche_res=0.1_size_filter"] # assert no nans, more niche labels than non-niche labels, and at least 100 obs per niche assert niches.isna().sum() == 0 @@ -35,9 +302,28 @@ def test_niche_calc_nhood(adata_seqfish: AnnData): # get obs x neighbor matrix from sparse matrix matrix = adata_seqfish.obsp[SPATIAL_CONNECTIVITIES_KEY].tocoo() + # since now this has become a method of an embedder class, will initialize the class and call the method for that class + embedder1 = NhoodProfileEmbedder( + groups=GROUPS, + spatial_connectivities_key="spatial_connectivities", + scale=True, # doesnt make a difference in this test function + distance=1, # doesnt make a difference in this test function + abs_nhood=True, + n_hop_weights=None, # doesnt make a difference in this test function + ) + + embedder2 = NhoodProfileEmbedder( + groups=GROUPS, + spatial_connectivities_key="spatial_connectivities", + scale=True, # doesnt make a difference in this test function + distance=1, # doesnt make a difference in this test function + abs_nhood=False, + n_hop_weights=None, # doesnt make a difference in this test function + ) + # get obs x category matrix where each column is the absolute/relative frequency of a category in the neighborhood - rel_nhood_profile = _calculate_neighborhood_profile(adata_seqfish, groups=GROUPS, matrix=matrix, abs_nhood=False) - abs_nhood_profile = _calculate_neighborhood_profile(adata_seqfish, groups=GROUPS, matrix=matrix, abs_nhood=True) + abs_nhood_profile = embedder1._calculate_neighborhood_profile(adata_seqfish, matrix) + rel_nhood_profile = embedder2._calculate_neighborhood_profile(adata_seqfish, matrix) # assert shape obs x groups assert rel_nhood_profile.shape == ( adata_seqfish.n_obs, @@ -62,26 +348,14 @@ def test_niche_calc_utag(adata_seqfish: AnnData): assert niches.isna().sum() == 0 assert niches.nunique() > niches_low_res.nunique() - # assert shape obs x var and sparsity in new feature matrix - new_feature_matrix = _utag( - adata_seqfish, - normalize_adj=True, - spatial_connectivity_key=SPATIAL_CONNECTIVITIES_KEY, - ) - assert new_feature_matrix.shape == adata_seqfish.X.shape - assert issparse(new_feature_matrix) + embedder = UtagEmbedder(spatial_connectivities_key=SPATIAL_CONNECTIVITIES_KEY) + embedding1 = embedder.get_embedding(adata_seqfish) # this is pca of the utag feature matrix + + # assert shape[0] of matrices same + assert embedding1.shape[0] == adata_seqfish.X.shape[0] spatial_neighbors_knn(adata_seqfish, n_neighs=40) - new_feature_matrix_more_neighs = _utag( - adata_seqfish, - normalize_adj=True, - spatial_connectivity_key=SPATIAL_CONNECTIVITIES_KEY, - ) + embedding2 = embedder.get_embedding(adata_seqfish) # this is pca of the utag feature matrix - # matrix products should differ when using different amount of neighbors - try: - assert_frame_equal(new_feature_matrix, new_feature_matrix_more_neighs) - except AssertionError: - pass - else: - raise AssertionError + # matrix products (and hence pca) should differ when using different amount of neighbors + assert (embedding1 != embedding2).all()