diff --git a/docs/api/scanpy_gpu.md b/docs/api/scanpy_gpu.md index e3d8b9bed..f680fd079 100644 --- a/docs/api/scanpy_gpu.md +++ b/docs/api/scanpy_gpu.md @@ -2,6 +2,23 @@ These functions offer accelerated near drop-in replacements for common tools provided by [`scanpy`](https://scanpy.readthedocs.io/en/stable/api/index.html) {cite}`Wolf2018`. +## Scanpy backend + +With Scanpy versions that support computational backends, RAPIDS-singlecell is +available as the `rapids-singlecell` backend with the aliases `cuda`, `rapids`, +and `rapids_singlecell`. + +```python +import scanpy as sc + +sc.settings.backend = "cuda" +``` + +The backend exposes RAPIDS-singlecell's `pp` and `tl` functions, plus {func}`rapids_singlecell.get.aggregate`, for Scanpy's backend dispatcher. +The data must already be on the GPU: move it with {func}`~rapids_singlecell.get.anndata_to_GPU` before calling the accelerated functions (see {doc}`/usage_principles`). Scanpy's `copy` argument, and `subset`/`inplace` of {func}`~rapids_singlecell.pp.highly_variable_genes`, behave as in Scanpy; other arguments that RAPIDS-singlecell does not support are dropped with a warning. + +Backend calls use RAPIDS-singlecell's native defaults and return types. In particular, {func}`~rapids_singlecell.pp.normalize_total` and {func}`~rapids_singlecell.pp.normalize_pearson_residuals` return the normalized GPU matrix directly with `inplace=False`, rather than Scanpy's result dictionary. {func}`~rapids_singlecell.pp.calculate_qc_metrics` defaults to `inplace=False` through the backend, matching Scanpy. + ## Preprocessing `pp` Filtering of highly-variable genes, batch-effect correction, per-cell normalization. diff --git a/docs/api/squidpy_gpu.md b/docs/api/squidpy_gpu.md index c11e5922c..eb82bf58b 100644 --- a/docs/api/squidpy_gpu.md +++ b/docs/api/squidpy_gpu.md @@ -3,6 +3,24 @@ {mod}`squidpy.gr` is a tool for the analysis of spatial molecular data {cite}`Palla2022`. {mod}`rapids_singlecell.gr` accelerates some of these functions. +## Squidpy backend + +With Squidpy versions that support computational backends, RAPIDS-singlecell is +available as the `rapids-singlecell` backend with the aliases `cuda`, `rapids`, +and `rapids_singlecell`. + +```python +import squidpy as sq + +sq.settings.backend = "cuda" +``` + +The backend exposes RAPIDS-singlecell's {mod}`rapids_singlecell.gr` functions for Squidpy's backend dispatcher, including the flavor-specific {func}`~rapids_singlecell.gr.calculate_niche_neighborhood`, {func}`~rapids_singlecell.gr.calculate_niche_utag` and {func}`~rapids_singlecell.gr.calculate_niche_cellcharter`. + +Like their Squidpy counterparts, {func}`~rapids_singlecell.gr.spatial_autocorr`, {func}`~rapids_singlecell.gr.co_occurrence`, {func}`~rapids_singlecell.gr.ligrec` and the `calculate_niche_*` functions accept a {class}`~spatialdata.SpatialData` together with `table_key`. `flavor="spatialleiden"` of {func}`squidpy.gr.calculate_niche` has no GPU implementation and raises `NotImplementedError` on this backend; `calculate_niche_spatialleiden` is not exposed, so Squidpy runs it on the CPU. + +Niche calls use RAPIDS-singlecell's native defaults and validation. The deprecated {func}`~rapids_singlecell.gr.calculate_niche` supplies `n_neighbors=15` and `resolutions=(0.5,)` when omitted and emits its native deprecation warning; prefer the flavor-specific functions. The backend translates Squidpy's `data`, `copy`, and `rng` arguments to the native names. + ```{eval-rst} .. module:: rapids_singlecell.gr .. currentmodule:: rapids_singlecell diff --git a/docs/conf.py b/docs/conf.py index 8f4a4bf2f..e532ba2c0 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -81,6 +81,7 @@ "pylibraft", "dask", "cuvs", + "spatialdata", ] default_role = "literal" napoleon_google_docstring = False @@ -128,6 +129,7 @@ "statsmodels": ("https://www.statsmodels.org/stable/", None), "omnipath": ("https://omnipath.readthedocs.io/en/latest/", None), "dask": ("https://docs.dask.org/en/stable/", None), + "spatialdata": ("https://spatialdata.scverse.org/en/stable/", None), } # List of patterns, relative to source directory, that match files and diff --git a/docs/release-notes/0.18.0.md b/docs/release-notes/0.18.0.md new file mode 100644 index 000000000..953ffee8f --- /dev/null +++ b/docs/release-notes/0.18.0.md @@ -0,0 +1,15 @@ +### 0.18.0 {small}`the-future` + +```{rubric} Bug fixes +``` +* Fix {func}`~rapids_singlecell.gr.ligrec` raising ``KeyError`` when the cluster column has integer categories {pr}`628` {smaller}`S Dicks` + +```{rubric} Features +``` +* Register RAPIDS-singlecell as a computational backend for Scanpy and Squidpy through the ``scanpy.backends`` and ``squidpy.backends`` entry points. Select it with ``sc.settings.backend = "cuda"`` / ``sq.settings.backend = "cuda"`` (aliases ``rapids`` and ``rapids_singlecell``) or per call with ``backend="cuda"``. All {mod}`~rapids_singlecell.pp`, {mod}`~rapids_singlecell.tl` and {mod}`~rapids_singlecell.gr` functions plus {func}`~rapids_singlecell.get.aggregate` are exposed, including the flavor-specific niche functions. Backend calls use RAPIDS-singlecell's native defaults and return types {pr}`628` {smaller}`S Dicks` +* {func}`~rapids_singlecell.gr.spatial_autocorr`, {func}`~rapids_singlecell.gr.co_occurrence`, {func}`~rapids_singlecell.gr.ligrec` and the ``calculate_niche_*`` functions accept a {class}`~spatialdata.SpatialData` together with ``table_key``, following {mod}`squidpy` {pr}`628` {smaller}`S Dicks` + +```{rubric} Misc +``` +* Add native ``copy`` support to {func}`~rapids_singlecell.pp.filter_cells`, {func}`~rapids_singlecell.pp.filter_genes`, {func}`~rapids_singlecell.pp.regress_out`, {func}`~rapids_singlecell.tl.diffmap`, {func}`~rapids_singlecell.tl.draw_graph` and {func}`~rapids_singlecell.tl.rank_genes_groups`, and ``subset``/``inplace`` support to {func}`~rapids_singlecell.pp.highly_variable_genes` {pr}`628` {smaller}`S Dicks` +* Scanpy compat: {func}`~rapids_singlecell.pp.calculate_qc_metrics` accepts ``inplace``; with ``inplace=False`` it returns the per-cell and per-gene metrics as two DataFrames instead of writing them to ``adata`` {pr}`628` {smaller}`S Dicks` diff --git a/docs/release-notes/index.md b/docs/release-notes/index.md index 7b6dacfb8..aeaf69a5a 100644 --- a/docs/release-notes/index.md +++ b/docs/release-notes/index.md @@ -2,6 +2,10 @@ # Release notes +## Version 0.18.0 +```{include} /release-notes/0.18.0.md +``` + ## Version 0.17.0 ```{include} /release-notes/0.17.0.md ``` diff --git a/pyproject.toml b/pyproject.toml index 9b39dee3b..f91eedecc 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -84,6 +84,12 @@ dev = [ "pre-commit", ] +[project.entry-points."scanpy.backends"] +rapids-singlecell = "rapids_singlecell._backends.scanpy" + +[project.entry-points."squidpy.backends"] +rapids-singlecell = "rapids_singlecell._backends.squidpy" + [project.urls] Documentation = "https://rapids-singlecell.readthedocs.io" Source = "https://github.com/scverse/rapids_singlecell" diff --git a/src/rapids_singlecell/_backends/__init__.py b/src/rapids_singlecell/_backends/__init__.py new file mode 100644 index 000000000..9d48db4f9 --- /dev/null +++ b/src/rapids_singlecell/_backends/__init__.py @@ -0,0 +1 @@ +from __future__ import annotations diff --git a/src/rapids_singlecell/_backends/scanpy.py b/src/rapids_singlecell/_backends/scanpy.py new file mode 100644 index 000000000..662278d0e --- /dev/null +++ b/src/rapids_singlecell/_backends/scanpy.py @@ -0,0 +1,85 @@ +"""Scanpy backend exports with RAPIDS-singlecell's native semantics.""" + +from __future__ import annotations + +from functools import partial + +from rapids_singlecell.get import aggregate +from rapids_singlecell.preprocessing import ( + bbknn, + filter_cells, + filter_genes, + filter_highly_variable, + flag_gene_family, + harmony_integrate, + highly_variable_genes, + log1p, + neighbors, + normalize_pearson_residuals, + normalize_total, + pca, + regress_out, + scale, + scrublet, + scrublet_simulate_doublets, + sqrt, +) +from rapids_singlecell.preprocessing import ( + calculate_qc_metrics as _calculate_qc_metrics, +) +from rapids_singlecell.tools import ( + diffmap, + draw_graph, + embedding_density, + ingest, + kmeans, + leiden, + louvain, + rank_genes_groups, + rank_genes_groups_logreg, + score_genes, + score_genes_cell_cycle, + tsne, + umap, +) + +name = "rapids-singlecell" +aliases = ["cuda", "rapids", "rapids_singlecell"] + +calculate_qc_metrics = partial(_calculate_qc_metrics, inplace=False) + + +__all__ = [ + "aggregate", + "bbknn", + "calculate_qc_metrics", + "diffmap", + "draw_graph", + "embedding_density", + "filter_cells", + "filter_genes", + "filter_highly_variable", + "flag_gene_family", + "harmony_integrate", + "highly_variable_genes", + "ingest", + "kmeans", + "leiden", + "log1p", + "louvain", + "neighbors", + "normalize_pearson_residuals", + "normalize_total", + "pca", + "rank_genes_groups", + "rank_genes_groups_logreg", + "regress_out", + "scale", + "score_genes", + "score_genes_cell_cycle", + "scrublet", + "scrublet_simulate_doublets", + "sqrt", + "tsne", + "umap", +] diff --git a/src/rapids_singlecell/_backends/squidpy.py b/src/rapids_singlecell/_backends/squidpy.py new file mode 100644 index 000000000..e22194806 --- /dev/null +++ b/src/rapids_singlecell/_backends/squidpy.py @@ -0,0 +1,87 @@ +"""Squidpy backend using native defaults and small niche argument bridges.""" + +from __future__ import annotations + +import inspect +import warnings +from functools import wraps + +from rapids_singlecell._utils._random import _seed_from_rng +from rapids_singlecell.squidpy_gpu import calculate_niche as _calculate_niche +from rapids_singlecell.squidpy_gpu import calculate_niche_cellcharter as _cellcharter +from rapids_singlecell.squidpy_gpu import calculate_niche_neighborhood as _neighborhood +from rapids_singlecell.squidpy_gpu import calculate_niche_utag as _utag +from rapids_singlecell.squidpy_gpu import co_occurrence, spatial_autocorr +from rapids_singlecell.squidpy_gpu import ligrec as _ligrec + +name = "rapids-singlecell" +aliases = ["cuda", "rapids", "rapids_singlecell"] + +SQUIDPY_LIGREC_CPU_ONLY = { + "seed": None, + "n_jobs": None, + "show_progress_bar": True, + "numba_parallel": None, +} + + +@wraps(_ligrec) +def ligrec(*args, **kwargs): + ignored = [] + for key, default in SQUIDPY_LIGREC_CPU_ONLY.items(): + if key in kwargs and kwargs.pop(key) != default: + ignored.append(key) + if ignored: + warnings.warn( + f"Parameters {', '.join(ignored)} have no effect on the " + "rapids-singlecell backend.", + UserWarning, + stacklevel=2, + ) + return _ligrec(*args, **kwargs) + + +def _niche_adapter(func, *, flavor_api: bool = False): + """Expose native niche signatures with Squidpy's argument names.""" + names = {"adata": "data"} + defaults = {} + if flavor_api: + names.update(inplace="copy", random_state="rng") + defaults.update(copy=False, rng=None) + + @wraps(func) + def adapted(data, **kwargs): + if kwargs.get("flavor") == "spatialleiden": + raise NotImplementedError("Use `backend='cpu'` for flavor='spatialleiden'.") + if flavor_api: + if "copy" in kwargs: + kwargs["inplace"] = not kwargs.pop("copy") + if (rng := kwargs.pop("rng", None)) is not None: + kwargs["random_state"] = _seed_from_rng(rng, allow_none=False) + return func(data, **kwargs) + + signature = inspect.signature(func) + parameters = [] + for parameter in signature.parameters.values(): + name = names.get(parameter.name, parameter.name) + parameters.append( + parameter.replace(name=name, default=defaults.get(name, parameter.default)) + ) + adapted.__signature__ = signature.replace(parameters=parameters) + return adapted + + +calculate_niche = _niche_adapter(_calculate_niche) +calculate_niche_neighborhood = _niche_adapter(_neighborhood, flavor_api=True) +calculate_niche_utag = _niche_adapter(_utag, flavor_api=True) +calculate_niche_cellcharter = _niche_adapter(_cellcharter, flavor_api=True) + +__all__ = [ + "calculate_niche", + "calculate_niche_cellcharter", + "calculate_niche_neighborhood", + "calculate_niche_utag", + "co_occurrence", + "ligrec", + "spatial_autocorr", +] diff --git a/src/rapids_singlecell/preprocessing/_hvg/__init__.py b/src/rapids_singlecell/preprocessing/_hvg/__init__.py index 24c0b326d..c77145402 100644 --- a/src/rapids_singlecell/preprocessing/_hvg/__init__.py +++ b/src/rapids_singlecell/preprocessing/_hvg/__init__.py @@ -3,8 +3,11 @@ from typing import TYPE_CHECKING, Literal import numpy as np +from anndata import AnnData +from cupyx.scipy.sparse import issparse from rapids_singlecell._settings import Default, resolve_default +from rapids_singlecell.get import _get_obs_rep, _set_obs_rep from rapids_singlecell.preprocessing._utils import _sanitize_column from ._cutoffs import _Cutoffs @@ -17,7 +20,7 @@ from ._seurat_v3 import _highly_variable_genes_seurat_v3 if TYPE_CHECKING: - from anndata import AnnData + import pandas as pd flavors = Literal[ "seurat", @@ -47,7 +50,9 @@ def highly_variable_genes( chunksize: int = 1000, n_samples: int = 10000, batch_key: str | None = None, -) -> None: + subset: bool = False, + inplace: bool = True, +) -> pd.DataFrame | None: """\ Annotate highly variable genes :cite:p:`Satija2015,Zheng2017,Stuart2019,Lause2021,Andrews2019`. @@ -117,10 +122,15 @@ def highly_variable_genes( of enrichment of zeros for each gene (only for `flavor='poisson_gene_selection'`). batch_key If specified, highly-variable genes are selected within each batch separately and merged. + subset + Restrict the result to highly variable genes. + inplace + Write metrics to `adata.var`. If `False`, return the metrics as a DataFrame. Returns ------- - updates `adata.var` with the following fields: + Returns the metrics if `inplace=False`; otherwise updates `adata.var` + with the following fields and subsets `adata` if requested: `highly_variable` : bool boolean indicator of highly-variable genes @@ -148,6 +158,17 @@ def highly_variable_genes( """ flavor = resolve_default(flavor) + if not inplace: + adata = AnnData( + X=adata.X, + obs=adata.obs.copy(), + var=adata.var[[]].copy(), + layers=adata.layers, + ) + X = _get_obs_rep(adata, layer=layer) + if issparse(X) and not X.has_canonical_format: + _set_obs_rep(adata, X.copy(), layer=layer) + if batch_key is not None: _sanitize_column(adata, batch_key) @@ -221,3 +242,9 @@ def highly_variable_genes( adata.var["highly_variable_intersection"] = df[ "highly_variable_intersection" ] + + if not inplace: + return adata.var.loc[adata.var["highly_variable"]] if subset else adata.var + if subset: + adata._inplace_subset_var(adata.var["highly_variable"]) + return None diff --git a/src/rapids_singlecell/preprocessing/_qc.py b/src/rapids_singlecell/preprocessing/_qc.py index 249ce3a04..d6889385f 100644 --- a/src/rapids_singlecell/preprocessing/_qc.py +++ b/src/rapids_singlecell/preprocessing/_qc.py @@ -3,6 +3,7 @@ from typing import TYPE_CHECKING import cupy as cp +import pandas as pd from cupyx.scipy import sparse from rapids_singlecell._compat import DaskArray @@ -24,7 +25,8 @@ def calculate_qc_metrics( qc_vars: str | list = None, log1p: bool = True, layer: str = None, -) -> None: + inplace: bool = True, +) -> tuple[pd.DataFrame, pd.DataFrame] | None: """\ Calculates basic qc Parameters :cite:p:`McCarthy2017`. @@ -46,9 +48,13 @@ def calculate_qc_metrics( Set to `False` to skip computing `log1p` transformed annotations. layer If provided, use :attr:`~anndata.AnnData.layers` for expression values instead of :attr:`~anndata.AnnData.X`. + inplace + Whether to place the calculated metrics in :attr:`~anndata.AnnData.obs` and :attr:`~anndata.AnnData.var`. + If `False`, return them as two :class:`~pandas.DataFrame` instead, like :func:`scanpy.pp.calculate_qc_metrics`. Returns ------- + If `inplace = False`, returns `(obs_metrics, var_metrics)`. Otherwise adds the following columns in :attr:`~anndata.AnnData.obs` : `total_{var_type}_by_{expr_type}` E.g. 'total_genes_by_counts'. Number of genes with positive counts in a cell. @@ -78,24 +84,26 @@ def calculate_qc_metrics( sums_cells, sums_genes, genes_per_cell, cells_per_gene = _basic_qc(X) # .var - adata.var[f"n_cells_by_{expr_type}"] = cp.asnumpy(cells_per_gene) - adata.var[f"total_{expr_type}"] = cp.asnumpy(sums_genes) + var_metrics = pd.DataFrame(index=adata.var_names) + var_metrics[f"n_cells_by_{expr_type}"] = cp.asnumpy(cells_per_gene) + var_metrics[f"total_{expr_type}"] = cp.asnumpy(sums_genes) mean_array = sums_genes / adata.n_obs - adata.var[f"mean_{expr_type}"] = cp.asnumpy(mean_array) - adata.var[f"pct_dropout_by_{expr_type}"] = cp.asnumpy( + var_metrics[f"mean_{expr_type}"] = cp.asnumpy(mean_array) + var_metrics[f"pct_dropout_by_{expr_type}"] = cp.asnumpy( (1 - cells_per_gene / adata.n_obs) * 100 ) if log1p: - adata.var[f"log1p_total_{expr_type}"] = cp.asnumpy(cp.log1p(sums_genes)) - adata.var[f"log1p_mean_{expr_type}"] = cp.asnumpy(cp.log1p(mean_array)) + var_metrics[f"log1p_total_{expr_type}"] = cp.asnumpy(cp.log1p(sums_genes)) + var_metrics[f"log1p_mean_{expr_type}"] = cp.asnumpy(cp.log1p(mean_array)) # .obs - adata.obs[f"n_{var_type}_by_{expr_type}"] = cp.asnumpy(genes_per_cell) - adata.obs[f"total_{expr_type}"] = cp.asnumpy(sums_cells) + obs_metrics = pd.DataFrame(index=adata.obs_names) + obs_metrics[f"n_{var_type}_by_{expr_type}"] = cp.asnumpy(genes_per_cell) + obs_metrics[f"total_{expr_type}"] = cp.asnumpy(sums_cells) if log1p: - adata.obs[f"log1p_n_{var_type}_by_{expr_type}"] = cp.asnumpy( + obs_metrics[f"log1p_n_{var_type}_by_{expr_type}"] = cp.asnumpy( cp.log1p(genes_per_cell) ) - adata.obs[f"log1p_total_{expr_type}"] = cp.asnumpy(cp.log1p(sums_cells)) + obs_metrics[f"log1p_total_{expr_type}"] = cp.asnumpy(cp.log1p(sums_cells)) if qc_vars: if isinstance(qc_vars, str): @@ -104,15 +112,21 @@ def calculate_qc_metrics( mask = cp.array(adata.var[qc_var], dtype=cp.bool_) sums_cells_sub = _geneset_qc(X, mask) - adata.obs[f"total_{expr_type}_{qc_var}"] = cp.asnumpy(sums_cells_sub) - adata.obs[f"pct_{expr_type}_{qc_var}"] = cp.asnumpy( + obs_metrics[f"total_{expr_type}_{qc_var}"] = cp.asnumpy(sums_cells_sub) + obs_metrics[f"pct_{expr_type}_{qc_var}"] = cp.asnumpy( sums_cells_sub / sums_cells * 100 ) if log1p: - adata.obs[f"log1p_total_{expr_type}_{qc_var}"] = cp.asnumpy( + obs_metrics[f"log1p_total_{expr_type}_{qc_var}"] = cp.asnumpy( cp.log1p(sums_cells_sub) ) + if not inplace: + return obs_metrics, var_metrics + adata.obs[obs_metrics.columns] = obs_metrics + adata.var[var_metrics.columns] = var_metrics + return None + def _basic_qc( X: ArrayTypesDask, diff --git a/src/rapids_singlecell/preprocessing/_regress_out.py b/src/rapids_singlecell/preprocessing/_regress_out.py index 8baf60cf0..19beb302c 100644 --- a/src/rapids_singlecell/preprocessing/_regress_out.py +++ b/src/rapids_singlecell/preprocessing/_regress_out.py @@ -26,7 +26,8 @@ def regress_out( inplace: bool = True, batchsize: int | Literal["all"] | None = None, verbose: bool = False, -) -> Union[cp.ndarray, None]: # noqa: UP007 + copy: bool = False, +) -> Union[cp.ndarray, AnnData, None]: # noqa: UP007 """ Use linear regression to adjust for the effects of unwanted noise and variation. @@ -54,6 +55,8 @@ def regress_out( verbose Print debugging information + copy + Return a corrected AnnData copy when `inplace=True`. Returns ------- @@ -64,6 +67,7 @@ def regress_out( raise ValueError("batchsize must be `int`, `None` or `'all'`") if isinstance(adata, AnnData): + adata = adata.copy() if copy and inplace else adata view_to_actual(adata) X = _get_obs_rep(adata, layer=layer) @@ -99,6 +103,7 @@ def regress_out( if inplace: _set_obs_rep(adata, X, layer=layer) + return adata if copy else None else: return X diff --git a/src/rapids_singlecell/preprocessing/_simple.py b/src/rapids_singlecell/preprocessing/_simple.py index 650e92bda..2bff54a7a 100644 --- a/src/rapids_singlecell/preprocessing/_simple.py +++ b/src/rapids_singlecell/preprocessing/_simple.py @@ -63,7 +63,8 @@ def filter_genes( max_cells: int | None = None, inplace: bool = True, verbose: bool = True, -) -> tuple[np.ndarray, np.ndarray] | None: + copy: bool = False, +) -> tuple[np.ndarray, np.ndarray] | AnnData | None: """\ Filter genes based on number of cells or counts. @@ -91,6 +92,8 @@ def filter_genes( Perform computation inplace or return result. verbose Print number of discarded genes + copy + Return a filtered AnnData copy when `inplace=True`. Returns ------- @@ -154,9 +157,11 @@ def filter_genes( print(msg) if isinstance(data, AnnData) and inplace: + data = data.copy() if copy else data col = "n_counts" if (min_cells is None and max_cells is None) else "n_cells" data.var[col] = number_per_gene.get() data._inplace_subset_var(gene_subset.get()) + return data if copy else None else: return gene_subset.get(), number_per_gene.get() @@ -170,7 +175,8 @@ def filter_cells( max_genes: int | None = None, inplace: bool = True, verbose: bool = True, -) -> tuple[np.ndarray, np.ndarray] | None: + copy: bool = False, +) -> tuple[np.ndarray, np.ndarray] | AnnData | None: """\ Filter cell outliers based on counts and numbers of genes expressed. @@ -198,6 +204,8 @@ def filter_cells( Perform computation inplace or return result. verbose Print number of discarded cells + copy + Return a filtered AnnData copy when `inplace=True`. Returns ------- @@ -294,9 +302,11 @@ def filter_cells( print(msg) if isinstance(data, AnnData) and inplace: + data = data.copy() if copy else data col = "n_counts" if (min_genes is None and max_genes is None) else "n_genes" data.obs[col] = number_per_cell.get() data._inplace_subset_obs(cell_subset.get()) + return data if copy else None else: return cell_subset.get(), number_per_cell.get() diff --git a/src/rapids_singlecell/squidpy_gpu/_autocorr.py b/src/rapids_singlecell/squidpy_gpu/_autocorr.py index 5010079d8..6accc342b 100644 --- a/src/rapids_singlecell/squidpy_gpu/_autocorr.py +++ b/src/rapids_singlecell/squidpy_gpu/_autocorr.py @@ -16,12 +16,13 @@ from ._gearysc import _gearys_C_cupy from ._moransi import _morans_I_cupy -from ._utils import _p_value_calc +from ._utils import _extract_adata_if_sdata, _p_value_calc if TYPE_CHECKING: from collections.abc import Sequence from anndata import AnnData + from spatialdata import SpatialData def _to_cupy(vals, *, use_sparse: bool, dtype): @@ -49,7 +50,7 @@ def _to_cupy(vals, *, use_sparse: bool, dtype): def spatial_autocorr( - adata: AnnData, + adata: AnnData | SpatialData, *, connectivity_key: str = "spatial_connectivities", genes: str | Sequence[str] | None = None, @@ -64,6 +65,7 @@ def spatial_autocorr( dtype: np.dtype | None = None, multi_gpu: bool | list[int] | str | None = None, copy: bool = False, + table_key: str | None = None, ) -> pd.DataFrame | None: """ Calculate spatial autocorrelation for genes in an AnnData object. @@ -79,7 +81,7 @@ def spatial_autocorr( Parameters ---------- adata - Annotated data matrix. + Annotated data matrix or :class:`~spatialdata.SpatialData`. connectivity_key Key of the connectivity matrix in `adata.obsp`, by default "spatial_connectivities". genes @@ -112,12 +114,16 @@ def spatial_autocorr( - str: Comma-separated GPU IDs (e.g., "0,2") copy If True, return the results as a DataFrame instead of storing them in `adata.uns`, by default False. + table_key + Key in :attr:`spatialdata.SpatialData.tables` selecting the table to use. + Required if ``adata`` is a :class:`~spatialdata.SpatialData`. Returns ------- DataFrame containing the autocorrelation scores, p-values, and corrected p-values for each gene. \ If `copy` is False, the results are stored in `adata.uns` and None is returned. """ + adata = _extract_adata_if_sdata(adata, table_key=table_key) if genes is None: if "highly_variable" in adata.var: genes = adata[:, adata.var["highly_variable"]].var_names.values diff --git a/src/rapids_singlecell/squidpy_gpu/_co_oc.py b/src/rapids_singlecell/squidpy_gpu/_co_oc.py index 1e21e5d82..e0e33e76b 100644 --- a/src/rapids_singlecell/squidpy_gpu/_co_oc.py +++ b/src/rapids_singlecell/squidpy_gpu/_co_oc.py @@ -15,20 +15,26 @@ parse_device_ids, ) -from ._utils import _assert_categorical_obs, _assert_spatial_basis +from ._utils import ( + _assert_categorical_obs, + _assert_spatial_basis, + _extract_adata_if_sdata, +) if TYPE_CHECKING: from anndata import AnnData + from spatialdata import SpatialData def co_occurrence( - adata: AnnData, + adata: AnnData | SpatialData, cluster_key: str, *, spatial_key: str = "spatial", interval: int | np.ndarray | cp.ndarray = 50, multi_gpu: bool | list[int] | str | None = None, copy: bool = False, + table_key: str | None = None, ) -> tuple[np.ndarray, np.ndarray] | None: """ Compute co-occurrence probability of clusters. @@ -36,7 +42,7 @@ def co_occurrence( Parameters ---------- adata - Annotated data object. + Annotated data object or :class:`~spatialdata.SpatialData`. cluster_key Key for the cluster labels. spatial_key @@ -53,6 +59,9 @@ def co_occurrence( - str: Comma-separated GPU IDs (e.g., "0,2") copy If ``True``, return the co-occurrence probability and the distance thresholds intervals. + table_key + Key in :attr:`spatialdata.SpatialData.tables` selecting the table to use. + Required if ``adata`` is a :class:`~spatialdata.SpatialData`. Returns ------- @@ -66,6 +75,7 @@ def co_occurrence( computed at ``interval``. """ + adata = _extract_adata_if_sdata(adata, table_key=table_key) _assert_categorical_obs(adata, key=cluster_key) _assert_spatial_basis(adata, key=spatial_key) spatial = cp.array(adata.obsm[spatial_key]).astype(np.float32) diff --git a/src/rapids_singlecell/squidpy_gpu/_ligrec.py b/src/rapids_singlecell/squidpy_gpu/_ligrec.py index f1b730b81..3c8a8434d 100644 --- a/src/rapids_singlecell/squidpy_gpu/_ligrec.py +++ b/src/rapids_singlecell/squidpy_gpu/_ligrec.py @@ -3,6 +3,7 @@ from collections.abc import Iterable, Mapping, Sequence from itertools import product from typing import ( + TYPE_CHECKING, Literal, ) @@ -14,7 +15,14 @@ from cupyx.scipy.sparse import issparse as cpissparse from scipy.sparse import csc_matrix, issparse -from ._utils import _assert_categorical_obs, _create_sparse_df +from ._utils import ( + _assert_categorical_obs, + _create_sparse_df, + _extract_adata_if_sdata, +) + +if TYPE_CHECKING: + from spatialdata import SpatialData SOURCE = "source" TARGET = "target" @@ -118,7 +126,7 @@ def _check_tuple_needles(needles, haystack, *, msg: str, reraise: bool = True): def ligrec( - adata: AnnData, + adata: AnnData | SpatialData, cluster_key: str, *, clusters: list | None = None, @@ -136,6 +144,7 @@ def ligrec( interactions_params: dict = {}, transmitter_params: dict = {"categories": "ligand"}, receiver_params: dict = {"categories": "receptor"}, + table_key: str | None = None, ) -> pd.DataFrame | None: """\ Perform the permutation test as described in [Efremova et al., 2020]. @@ -143,7 +152,7 @@ def ligrec( Parameters ---------- adata - Annotated data object. + Annotated data object or :class:`~spatialdata.SpatialData`. cluster_key Key in :attr:`~anndata.AnnData.obs` where clustering is stored. @@ -216,6 +225,10 @@ def ligrec( Keyword arguments for :func:`omnipath.interactions.import_intercell_network()` \ defining the receiver side of intercellular connections. + table_key + Key in :attr:`spatialdata.SpatialData.tables` selecting the table to use. \ + Required if ``adata`` is a :class:`~spatialdata.SpatialData`. + Returns ------- If `copy = True`, returns a dict with following keys: @@ -233,6 +246,7 @@ def ligrec( interacting components was 0 or it didn't pass the threshold percentage of \ cells being expressed within a given cluster. """ + adata = _extract_adata_if_sdata(adata, table_key=table_key) # Get and Check interactions if interactions is None: interactions = _get_interactions( @@ -416,7 +430,7 @@ def find_min_gene_in_complex(_complex: str | None) -> str | None: mat = adata.raw[filter_obs, filter_var].X else: mat = adata[filter_obs, filter_var].X - cluster_obs = adata.obs.loc[filter_obs, cluster_key] + cluster_obs = filtered_data["clusters"][filter_obs] cluster_obs = cluster_obs.cat.remove_unused_categories() cat = cluster_obs.cat diff --git a/src/rapids_singlecell/squidpy_gpu/_niche.py b/src/rapids_singlecell/squidpy_gpu/_niche.py index a451acee4..6e9458d5f 100644 --- a/src/rapids_singlecell/squidpy_gpu/_niche.py +++ b/src/rapids_singlecell/squidpy_gpu/_niche.py @@ -15,9 +15,13 @@ import rapids_singlecell as rsc from rapids_singlecell._keys import _embedding_keys +from ._utils import _extract_adata_if_sdata + if TYPE_CHECKING: from collections.abc import Callable, Sequence + from spatialdata import SpatialData + __all__ = [ "calculate_niche", @@ -47,7 +51,7 @@ def _normalize_resolutions( def calculate_niche_neighborhood( - adata: AnnData, + adata: AnnData | SpatialData, *, groups: str, resolutions: float | Sequence[float], @@ -61,6 +65,7 @@ def calculate_niche_neighborhood( mask: pd.Series | None = None, library_key: str | None = None, inplace: bool = True, + table_key: str | None = None, ) -> AnnData | None: """\ Compute spatial niches from cell-type neighborhood profiles on the GPU. @@ -101,7 +106,11 @@ def calculate_niche_neighborhood( per sample and labels are prefixed with ``lib=_``. inplace Write the niche columns to ``adata``. If ``False``, return a modified copy. + table_key + Key in :attr:`spatialdata.SpatialData.tables` selecting the table to use. + Required if ``adata`` is a :class:`~spatialdata.SpatialData`. """ + adata = _extract_adata_if_sdata(adata, table_key=table_key) resolutions = _normalize_resolutions(resolutions) _check_key(adata, spatial_connectivities_key) if groups is None: @@ -138,7 +147,7 @@ def calculate_niche_neighborhood( def calculate_niche_utag( - adata: AnnData, + adata: AnnData | SpatialData, *, resolutions: float | Sequence[float], n_neighbors: int = 15, @@ -147,6 +156,7 @@ def calculate_niche_utag( mask: pd.Series | None = None, library_key: str | None = None, inplace: bool = True, + table_key: str | None = None, ) -> AnnData | None: """\ Compute spatial niches from UTAG-smoothed expression on the GPU. @@ -177,7 +187,11 @@ def calculate_niche_utag( per sample and labels are prefixed with ``lib=_``. inplace Write the niche columns to ``adata``. If ``False``, return a modified copy. + table_key + Key in :attr:`spatialdata.SpatialData.tables` selecting the table to use. + Required if ``adata`` is a :class:`~spatialdata.SpatialData`. """ + adata = _extract_adata_if_sdata(adata, table_key=table_key) resolutions = _normalize_resolutions(resolutions) _check_key(adata, spatial_connectivities_key) if n_neighbors < 1: @@ -200,7 +214,7 @@ def calculate_niche_utag( def calculate_niche_cellcharter( - adata: AnnData, + adata: AnnData | SpatialData, *, distance: int = 3, aggregation: Literal["mean", "variance"] = "mean", @@ -212,6 +226,7 @@ def calculate_niche_cellcharter( mask: pd.Series | None = None, library_key: str | None = None, inplace: bool = True, + table_key: str | None = None, ) -> AnnData | None: """\ Compute spatial niches with the CellCharter approach on the GPU. @@ -251,7 +266,11 @@ def calculate_niche_cellcharter( per sample and labels are prefixed with ``lib=_``. inplace Write the niche columns to ``adata``. If ``False``, return a modified copy. + table_key + Key in :attr:`spatialdata.SpatialData.tables` selecting the table to use. + Required if ``adata`` is a :class:`~spatialdata.SpatialData`. """ + adata = _extract_adata_if_sdata(adata, table_key=table_key) if use_rep is None: _check_key(adata, spatial_connectivities_key) if distance < 0: @@ -302,7 +321,7 @@ def calculate_niche_cellcharter( ) ) def calculate_niche( - adata: AnnData, + adata: AnnData | SpatialData, *, flavor: Literal["neighborhood", "utag", "cellcharter"], groups: str | None = None, @@ -322,6 +341,7 @@ def calculate_niche( random_state: int = 42, inplace: bool = True, copy: bool | None = None, + table_key: str | None = None, **kwargs, ) -> AnnData | None: """\ @@ -392,9 +412,12 @@ def calculate_niche( Write the niche columns to ``adata``. If ``False``, return a modified copy. copy Deprecated alias for ``inplace``; ``copy=True`` is ``inplace=False``. + table_key + Table to use when ``adata`` is a :class:`~spatialdata.SpatialData`. kwargs Accepts the removed ``gmm_init`` argument, which is ignored with a warning. """ + adata = _extract_adata_if_sdata(adata, table_key=table_key) if flavor not in FLAVORS: raise ValueError( f"Unknown flavor '{flavor}'. Use 'neighborhood', 'utag', or 'cellcharter'." diff --git a/src/rapids_singlecell/squidpy_gpu/_utils.py b/src/rapids_singlecell/squidpy_gpu/_utils.py index 17a3a73c1..7c575afc5 100644 --- a/src/rapids_singlecell/squidpy_gpu/_utils.py +++ b/src/rapids_singlecell/squidpy_gpu/_utils.py @@ -1,5 +1,6 @@ from __future__ import annotations +import sys from typing import ( TYPE_CHECKING, Any, @@ -12,6 +13,32 @@ from scipy import stats from scipy.sparse import issparse, spmatrix +if TYPE_CHECKING: + from anndata import AnnData + from spatialdata import SpatialData + + +def _extract_adata_if_sdata( + adata: AnnData | SpatialData, *, table_key: str | None = None +) -> AnnData: + """\ + Resolve a :class:`~spatialdata.SpatialData` to the table ``table_key``. + + An :class:`~anndata.AnnData` is returned unchanged. ``spatialdata`` is only + looked up in :data:`sys.modules`, so it is never imported here. + """ + spatialdata = sys.modules.get("spatialdata") + if spatialdata is None or not isinstance(adata, spatialdata.SpatialData): + return adata + if table_key is None: + raise TypeError("missing required keyword-only argument: 'table_key'") + if table_key not in adata.tables: + raise ValueError( + f"Table {table_key!r} not found in SpatialData. " + f"Available tables: {list(adata.tables)}" + ) + return adata.tables[table_key] + def _check_precision_issues(score: cp.ndarray, dtype: np.dtype) -> None: """Check for numerical issues in autocorrelation scores. diff --git a/src/rapids_singlecell/tools/_diffmap.py b/src/rapids_singlecell/tools/_diffmap.py index 2614bae74..ce6ac850d 100644 --- a/src/rapids_singlecell/tools/_diffmap.py +++ b/src/rapids_singlecell/tools/_diffmap.py @@ -161,7 +161,8 @@ def diffmap( sort: Literal["decrease", "increase"] = "decrease", density_normalize: bool = True, rng: SeedLike | RNGLike | None = None, -) -> None: + copy: bool = False, +) -> AnnData | None: """ Diffusion Maps :cite:p:`Coifman2005,Haghverdi2015`. @@ -190,6 +191,8 @@ def diffmap( rng Random seed or :class:`~numpy.random.Generator` for reproducibility. The superseded `random_state` argument is still accepted. + copy + Return an annotated copy instead of updating `adata`. Returns ------- @@ -202,6 +205,7 @@ def diffmap( Array of size (number of eigen vectors). Eigenvalues of transition matrix. """ + adata = adata.copy() if copy else adata rng = np.random.default_rng(rng) connectivities = _load_connectivities(adata, neighbors_key) transitions_sym, _Z = _compute_transitions( @@ -214,3 +218,4 @@ def diffmap( # the scanpy 1 key holds the bare evals, `key_added` holds a dict adata.uns[keys.uns] = evals.get() if key_added is None else {"evals": evals.get()} adata.obsm[keys.obsm] = evecs.get() + return adata if copy else None diff --git a/src/rapids_singlecell/tools/_draw_graph.py b/src/rapids_singlecell/tools/_draw_graph.py index abd99ed12..1abe5e02c 100644 --- a/src/rapids_singlecell/tools/_draw_graph.py +++ b/src/rapids_singlecell/tools/_draw_graph.py @@ -33,7 +33,8 @@ def draw_graph( max_iter: int = 500, rng: SeedLike | RNGLike | None = None, key_added: str | Default | None = Default(("draw_graph", "key_added")), -) -> None: + copy: bool = False, +) -> AnnData | None: """ Force-directed graph drawing :cite:p:`Fruchterman1991,Jacomy2014`. @@ -64,6 +65,8 @@ def draw_graph( The superseded `random_state` argument is still accepted. key_added Template controlling where coordinates and parameters are stored. + copy + Return an annotated copy instead of updating `adata`. Returns ------- @@ -72,6 +75,7 @@ def draw_graph( X_draw_graph_layout_fa : `adata.obsm` Coordinates of graph layout. """ + adata = adata.copy() if copy else adata rng = np.random.default_rng(rng) meta_random_state = {"random_state": rng.arg} if isinstance(rng, _LegacyRng) else {} @@ -132,3 +136,4 @@ def draw_graph( keys = _embedding_keys("draw_graph", key_added, layout=layout) adata.uns[keys.uns] = {"params": {"layout": layout, **meta_random_state}} adata.obsm[keys.obsm] = positions.get() # Format output + return adata if copy else None diff --git a/src/rapids_singlecell/tools/_rank_genes_groups/__init__.py b/src/rapids_singlecell/tools/_rank_genes_groups/__init__.py index 2f0095186..751a534ef 100644 --- a/src/rapids_singlecell/tools/_rank_genes_groups/__init__.py +++ b/src/rapids_singlecell/tools/_rank_genes_groups/__init__.py @@ -75,8 +75,9 @@ def rank_genes_groups( n_bins: int | None = None, bin_range: Literal["log1p", "auto"] | None = None, skip_empty_groups: bool = False, + copy: bool = False, **kwds, -) -> None: +) -> AnnData | None: """ Rank genes for characterizing groups using GPU acceleration. @@ -197,6 +198,8 @@ def rank_genes_groups( Skip selected groups with fewer than two observations after filtering. This is useful for perturbation workflows where a per-cell-type slice keeps categories that are empty or singleton in that slice. + copy + Return an annotated copy instead of updating `adata`. **kwds Additional arguments passed to the method. For `'logreg'`, these are passed to :class:`cuml.linear_model.LogisticRegression`. @@ -228,6 +231,7 @@ def rank_genes_groups( `adata.uns['rank_genes_groups' | key_added]['pts_rest']` Fraction of cells expressing genes in rest. Only if `pts=True` and `reference='rest'`. """ + adata = adata.copy() if copy else adata method = resolve_default(method) mean_in_log_space = resolve_default(mean_in_log_space) @@ -360,7 +364,7 @@ def rank_genes_groups( test_obj.pts_rest.T, index=test_obj.var_names, columns=groups_names ) - return None + return adata if copy else None @deprecated(Deprecation("0.14.1", "Use `rank_genes_groups(method='logreg')` instead.")) diff --git a/tests/test_autocorr.py b/tests/test_autocorr.py index 4fb504d9b..0d691e8e0 100644 --- a/tests/test_autocorr.py +++ b/tests/test_autocorr.py @@ -48,6 +48,23 @@ def test_autocorr_consistency(mode): assert not np.array_equal(idx_df, idx_adata) +def test_autocorr_spatialdata(): + """Check that a SpatialData table gives the same result as the AnnData.""" + spatialdata = pytest.importorskip("spatialdata") + file = Path(__file__).parent / Path("_data/dummy.h5ad") + dummy_adata = read_h5ad(file) + sdata = spatialdata.SpatialData(tables={"table": dummy_adata.copy()}) + + expected = spatial_autocorr(dummy_adata, mode="moran", copy=True) + actual = spatial_autocorr(sdata, mode="moran", copy=True, table_key="table") + + np.testing.assert_allclose(actual.to_numpy(), expected.to_numpy()) + spatial_autocorr(sdata, mode="moran", table_key="table") + assert MORAN_I in sdata.tables["table"].uns + with pytest.raises(TypeError, match="table_key"): + spatial_autocorr(sdata, mode="moran") + + @pytest.mark.parametrize("mode", ["moran", "geary"]) @pytest.mark.parametrize("dtype", [np.float32, np.float64]) def test_autocorr_sparse(mode, dtype): diff --git a/tests/test_backend_conformance.py b/tests/test_backend_conformance.py new file mode 100644 index 000000000..bfb11c09e --- /dev/null +++ b/tests/test_backend_conformance.py @@ -0,0 +1,29 @@ +"""Run scanpy's and squidpy's backend conformance suites against the RSC backend.""" + +from __future__ import annotations + +import pytest + +SCANPY_CONFORMANCE_BLOCKED = ( + "scanpy's conformance suite feeds CPU AnnData, which rsc.pp/tl reject, and its " + "6-cell dataset crashes cuML's spectral UMAP initialisation" +) + + +@pytest.mark.parametrize( + "module", + [ + pytest.param( + "scanpy.testing", marks=pytest.mark.skip(reason=SCANPY_CONFORMANCE_BLOCKED) + ), + "squidpy.testing.backend_conformance", + ], +) +def test_conformance(module): + testing = pytest.importorskip(module) + if not hasattr(testing, "validate_backend"): + pytest.skip(f"{module} has no backend conformance suite") + + results = testing.validate_backend("rapids-singlecell") + for name, status in results.items(): + assert status == "PASSED", f"{name}: {status}" diff --git a/tests/test_co_oc.py b/tests/test_co_oc.py index 6ef89fcb2..b3340e3c5 100644 --- a/tests/test_co_oc.py +++ b/tests/test_co_oc.py @@ -41,6 +41,24 @@ def test_co_occurrence(adata: AnnData): assert arr.shape[1] == arr.shape[0] == adata.obs["leiden"].unique().shape[0] +def test_co_occurrence_spatialdata(adata: AnnData): + """Check that a SpatialData table gives the same result as the AnnData.""" + spatialdata = pytest.importorskip("spatialdata") + sdata = spatialdata.SpatialData(tables={"table": adata.copy()}) + + expected = co_occurrence(adata, cluster_key="leiden", copy=True) + actual = co_occurrence(sdata, cluster_key="leiden", copy=True, table_key="table") + + for actual_value, expected_value in zip(actual, expected, strict=True): + np.testing.assert_allclose(actual_value, expected_value) + co_occurrence(sdata, cluster_key="leiden", table_key="table") + assert "leiden_co_occurrence" in sdata.tables["table"].uns + with pytest.raises(TypeError, match="table_key"): + co_occurrence(sdata, cluster_key="leiden") + with pytest.raises(ValueError, match="not found"): + co_occurrence(sdata, cluster_key="leiden", table_key="missing") + + def test_co_occurrence_reproducibility(adata: AnnData): """Check co_occurrence reproducibility results.""" arr_1, interval_1 = co_occurrence(adata, cluster_key="leiden", copy=True) diff --git a/tests/test_ligrec.py b/tests/test_ligrec.py index ebacbae57..987e3d73c 100644 --- a/tests/test_ligrec.py +++ b/tests/test_ligrec.py @@ -234,7 +234,12 @@ def test_pvals_in_correct_range( assert np.nanmax(r["pvalues"].values) <= 1.0, np.nanmax(r["pvalues"].values) assert np.nanmin(r["pvalues"].values) >= 0, np.nanmin(r["pvalues"].values) - def test_result_correct_index(self, adata: AnnData, interactions: Interactions_t): + @pytest.mark.parametrize("integer_categories", [False, True]) + def test_result_correct_index( + self, adata: AnnData, interactions: Interactions_t, integer_categories: bool + ): + if integer_categories: + adata.obs[_CK] = adata.obs[_CK].cat.codes.astype("category") r = ligrec(adata, _CK, interactions=interactions, n_perms=5, copy=True) np.testing.assert_array_equal(r["means"].index, r["pvalues"].index) @@ -243,6 +248,9 @@ def test_result_correct_index(self, adata: AnnData, interactions: Interactions_t np.testing.assert_array_equal(r["means"].columns, r["pvalues"].columns) assert not np.array_equal(r["means"].columns, r["metadata"].columns) assert not np.array_equal(r["pvalues"].columns, r["metadata"].columns) + assert set(r["means"].columns.get_level_values("cluster_1")) == set( + map(str, adata.obs[_CK].cat.categories) + ) def test_result_is_sparse(self, adata: AnnData, interactions: Interactions_t): interactions = pd.DataFrame(interactions, columns=["source", "target"]) diff --git a/tests/test_niche.py b/tests/test_niche.py index 7d79af08e..a9350425f 100644 --- a/tests/test_niche.py +++ b/tests/test_niche.py @@ -155,6 +155,22 @@ def test_copy_returns_new_object(adata): assert list(adata.obs.columns) == before +def test_spatialdata_table(adata): + """A SpatialData table gives the same niches as the AnnData and is updated in place.""" + spatialdata = pytest.importorskip("spatialdata") + sdata = spatialdata.SpatialData(tables={"table": adata.copy()}) + + calculate_niche_utag(sdata, resolutions=0.5, n_neighbors=10, table_key="table") + calculate_niche_utag(adata, resolutions=0.5, n_neighbors=10) + + np.testing.assert_array_equal( + sdata.tables["table"].obs["utag_niche_res=0.5"].astype(str), + adata.obs["utag_niche_res=0.5"].astype(str), + ) + with pytest.raises(TypeError, match="table_key"): + calculate_niche_utag(sdata, resolutions=0.5, n_neighbors=10) + + def test_multiple_resolutions(adata): calculate_niche( adata, diff --git a/tests/test_qc_metrics.py b/tests/test_qc_metrics.py index dba85cb0d..f87971b2b 100644 --- a/tests/test_qc_metrics.py +++ b/tests/test_qc_metrics.py @@ -55,6 +55,15 @@ def test_qc_metrics(dtype): assert np.allclose(cudata.obs[col], old_obs[col]) for col in cudata.var: assert np.allclose(cudata.var[col], old_var[col]) + # inplace=False returns the metrics and leaves adata untouched + plain = AnnData(X=cudata.X.copy(), var=cudata.var[["mito", "negative"]].copy()) + obs_metrics, var_metrics = rsc.pp.calculate_qc_metrics( + plain, qc_vars=["mito", "negative"], inplace=False + ) + assert list(plain.obs.columns) == [] + assert list(plain.var.columns) == ["mito", "negative"] + pd.testing.assert_frame_equal(obs_metrics, cudata.obs[obs_metrics.columns]) + pd.testing.assert_frame_equal(var_metrics, cudata.var[var_metrics.columns]) # with log1p=False cudata = AnnData( X=sparse_gpu.csr_matrix( diff --git a/tests/test_scanpy_backend.py b/tests/test_scanpy_backend.py new file mode 100644 index 000000000..573d4c718 --- /dev/null +++ b/tests/test_scanpy_backend.py @@ -0,0 +1,195 @@ +from __future__ import annotations + +import inspect +import tomllib +from copy import copy +from functools import wraps +from pathlib import Path + +import cupy as cp +import numpy as np +import pytest +import scanpy as sc +from anndata import AnnData + +import rapids_singlecell as rsc +from rapids_singlecell._backends import scanpy as scanpy_backend + + +def _public_functions(module): + return { + name + for name, value in vars(module).items() + if not name.startswith("_") and inspect.isfunction(value) + } + + +@pytest.mark.parametrize("module", [rsc.pp, rsc.tl]) +def test_scanpy_backend_exports_public_scanpy_api(module): + assert _public_functions(module) <= set(scanpy_backend.__all__) + + +def test_scanpy_backend_exports_aggregate(): + assert scanpy_backend.aggregate is rsc.get.aggregate + assert "aggregate" in scanpy_backend.__all__ + + +@pytest.mark.parametrize("name", ["log1p", "pca", "scale"]) +def test_public_pp_data_first_signature(name): + params = inspect.signature(getattr(rsc.pp, name)).parameters + first_param = next(iter(params)) + + assert first_param == "data" + assert "adata" not in params + + +@pytest.mark.parametrize( + "name", + [ + "filter_cells", + "filter_genes", + "regress_out", + "diffmap", + "draw_graph", + "rank_genes_groups", + ], +) +def test_scanpy_backend_adds_copy(name): + params = inspect.signature(getattr(scanpy_backend, name)).parameters + + assert params["copy"].default is False + assert set( + inspect.signature( + getattr(rsc.pp, name, None) or getattr(rsc.tl, name) + ).parameters + ) <= set(params) + + +def test_scanpy_backend_copy_returns_filtered_copy(): + adata = AnnData(cp.asarray(np.arange(1, 13, dtype=np.float32).reshape(4, 3))) + + # the dispatcher passes the data object by keyword + result = scanpy_backend.filter_cells(data=adata, min_counts=10, copy=True) + + assert result is not adata + assert (result.n_obs, adata.n_obs) == (3, 4) + assert scanpy_backend.filter_cells(adata, min_counts=10) is None + assert adata.n_obs == 3 + + +def test_scanpy_backend_hvg_subset_and_inplace(): + rng = np.random.default_rng(0) + adata = AnnData(cp.asarray(rng.poisson(1.0, (200, 50)).astype(np.float32))) + + metrics = scanpy_backend.highly_variable_genes(adata, n_top_genes=10, inplace=False) + + assert "highly_variable" not in adata.var + assert metrics["highly_variable"].sum() == 10 + assert ( + len( + scanpy_backend.highly_variable_genes( + adata, n_top_genes=10, inplace=False, subset=True + ) + ) + == 10 + ) + scanpy_backend.highly_variable_genes(adata, n_top_genes=10, subset=True) + assert adata.n_vars == 10 + + +def test_scanpy_backend_entrypoint_is_declared(): + pyproject = tomllib.loads( + (Path(__file__).resolve().parents[1] / "pyproject.toml").read_text() + ) + + assert pyproject["project"]["entry-points"]["scanpy.backends"] == { + "rapids-singlecell": "rapids_singlecell._backends.scanpy" + } + + +def test_scanpy_backend_dispatch_smoke(monkeypatch): + scanpy_backends = pytest.importorskip("scanpy._backends") + + registry = scanpy_backends.dispatcher._registry + dispatch_impl = scanpy_backends.dispatcher._dispatch_impl + old_backend = scanpy_backends.settings.backend + old_state = { + "_backends": copy(registry._backends), + "_alias_map": copy(registry._alias_map), + "_load_errors": copy(registry._load_errors), + "_registration_errors": copy(registry._registration_errors), + "_warned_untrusted": copy(registry._warned_untrusted), + "_discovered": registry._discovered, + "_sig_cache": copy(dispatch_impl._sig_cache), + } + + @wraps(scanpy_backend.normalize_total) + def fake_normalize_total(adata: AnnData, **kwargs) -> None: + adata.X *= kwargs["target_sum"] + adata.uns["scanpy_backend_called"] = "normalize_total" + + @wraps(scanpy_backend.scale) + def fake_scale(data: AnnData, **kwargs) -> None: + data.X *= 2 + data.uns["scanpy_scale_backend_called"] = kwargs + + @wraps(scanpy_backend.pca) + def fake_pca(data: AnnData, n_comps: int | None = None, **kwargs) -> None: + data.uns["scanpy_pca_backend_called"] = { + "n_comps": n_comps, + **kwargs, + } + + monkeypatch.setattr(scanpy_backend, "normalize_total", fake_normalize_total) + monkeypatch.setattr(scanpy_backend, "scale", fake_scale) + monkeypatch.setattr(scanpy_backend, "pca", fake_pca) + + try: + scanpy_backends.settings._backend_var.set("cpu") + registry._backends.clear() + registry._alias_map.clear() + registry._load_errors.clear() + registry._registration_errors.clear() + registry._warned_untrusted.clear() + registry._discovered = True + registry._register_backend( + scanpy_backend, + entrypoint_name="rapids-singlecell", + distribution_name="rapids-singlecell", + object_ref="rapids_singlecell._backends.scanpy", + ) + dispatch_impl._sig_cache.clear() + dispatch_impl._update_signatures() + + adata = AnnData(np.ones((2, 2), dtype=np.float32)) + sc.pp.normalize_total(adata, target_sum=3, backend="cuda") + sc.pp.scale(adata, max_value=5, backend="cuda") + sc.pp.pca(adata, n_comps=1, backend="cuda") + + np.testing.assert_allclose(adata.X, 6) + assert adata.uns["scanpy_backend_called"] == "normalize_total" + assert adata.uns["scanpy_scale_backend_called"]["max_value"] == 5 + assert adata.uns["scanpy_pca_backend_called"]["n_comps"] == 1 + + gpu = AnnData(cp.ones((2, 2), dtype=cp.float32)) + obs_metrics, var_metrics = sc.pp.calculate_qc_metrics(gpu, backend="cuda") + assert "total_counts" in obs_metrics and "total_counts" in var_metrics + assert "total_counts" not in gpu.obs + sc.pp.calculate_qc_metrics(gpu, inplace=True, backend="cuda") + assert "total_counts" in gpu.obs + finally: + scanpy_backends.settings._backend_var.set(old_backend) + registry._backends.clear() + registry._backends.update(old_state["_backends"]) + registry._alias_map.clear() + registry._alias_map.update(old_state["_alias_map"]) + registry._load_errors.clear() + registry._load_errors.update(old_state["_load_errors"]) + registry._registration_errors.clear() + registry._registration_errors.update(old_state["_registration_errors"]) + registry._warned_untrusted.clear() + registry._warned_untrusted.update(old_state["_warned_untrusted"]) + registry._discovered = old_state["_discovered"] + dispatch_impl._sig_cache.clear() + dispatch_impl._sig_cache.update(old_state["_sig_cache"]) + dispatch_impl._update_signatures() diff --git a/tests/test_squidpy_backend.py b/tests/test_squidpy_backend.py new file mode 100644 index 000000000..3c1898bd5 --- /dev/null +++ b/tests/test_squidpy_backend.py @@ -0,0 +1,237 @@ +from __future__ import annotations + +import inspect +import tomllib +from copy import copy +from functools import wraps +from pathlib import Path + +import numpy as np +import pytest +from anndata import AnnData + +import rapids_singlecell as rsc +from rapids_singlecell._backends import squidpy as squidpy_backend + + +def test_squidpy_backend_identity_and_exports(): + assert squidpy_backend.name == "rapids-singlecell" + assert squidpy_backend.aliases == ["cuda", "rapids", "rapids_singlecell"] + assert squidpy_backend.__all__ == [ + "calculate_niche", + "calculate_niche_cellcharter", + "calculate_niche_neighborhood", + "calculate_niche_utag", + "co_occurrence", + "ligrec", + "spatial_autocorr", + ] + + +def test_squidpy_backend_exports_public_gr_api(): + public = { + name + for name, value in vars(rsc.gr).items() + if not name.startswith("_") and inspect.isfunction(value) + } + + assert public <= set(squidpy_backend.__all__) + + +def test_squidpy_backend_entrypoint_is_declared(): + pyproject = tomllib.loads( + (Path(__file__).resolve().parents[1] / "pyproject.toml").read_text() + ) + + assert pyproject["project"]["entry-points"]["squidpy.backends"] == { + "rapids-singlecell": "rapids_singlecell._backends.squidpy" + } + + +def test_calculate_niche_routes_to_flavor_function(monkeypatch): + adata = AnnData(np.ones((2, 2), dtype=np.float32)) + captured = {} + + @wraps(squidpy_backend._utag) + def fake_utag(data, **kwargs): + captured["data"] = data + captured.update(kwargs) + return None if kwargs["inplace"] else data + + monkeypatch.setattr(squidpy_backend, "_utag", fake_utag) + + with pytest.warns(UserWarning, match="groups are not used for flavor 'utag'"): + result = squidpy_backend.calculate_niche( + adata, + flavor="utag", + groups="cluster", + n_neighbors=10, + resolutions=0.5, + layer_ratio=2.0, + inplace=False, + ) + + assert result is adata + assert captured.pop("data") is adata + assert captured == { + "spatial_connectivities_key": "spatial_connectivities", + "inplace": False, + "table_key": None, + "n_neighbors": 10, + "resolutions": 0.5, + } + + +def test_calculate_niche_cellcharter_maps_squidpy_conventions(monkeypatch): + adata = AnnData(np.ones((2, 2), dtype=np.float32)) + captured = {} + + @wraps(squidpy_backend._cellcharter) + def fake_cellcharter(data, **kwargs): + captured["data"] = data + captured.update(kwargs) + return data + + monkeypatch.setattr(squidpy_backend, "_cellcharter", fake_cellcharter) + + result = squidpy_backend.calculate_niche_cellcharter( + adata, distance=2, rng=np.random.default_rng(0), copy=True + ) + + assert result is adata + assert captured.pop("data") is adata + assert isinstance(captured.pop("random_state"), int) + assert captured == { + "inplace": False, + "distance": 2, + "aggregation": "mean", + "spatial_connectivities_key": "spatial_connectivities", + "n_components": 10, + "use_rep": None, + "min_niche_size": None, + "mask": None, + "library_key": None, + "table_key": None, + } + + +def test_calculate_niche_rejects_spatialleiden(): + adata = AnnData(np.ones((2, 2), dtype=np.float32)) + + with pytest.raises(NotImplementedError, match="spatialleiden"): + squidpy_backend.calculate_niche(adata, flavor="spatialleiden") + + +def test_calculate_niche_forwards_spatialdata_table(monkeypatch): + spatialdata = pytest.importorskip("spatialdata") + adata = AnnData(np.ones((2, 2), dtype=np.float32)) + sdata = spatialdata.SpatialData(tables={"table": adata}) + captured = {} + + @wraps(squidpy_backend._cellcharter) + def fake_cellcharter(data, **kwargs): + captured["data"] = data + captured.update(kwargs) + + monkeypatch.setattr(squidpy_backend, "_cellcharter", fake_cellcharter) + + squidpy_backend.calculate_niche( + sdata, flavor="cellcharter", distance=2, aggregation="mean", table_key="table" + ) + + assert captured["data"] is sdata + assert captured["table_key"] == "table" + + +def test_calculate_niche_requires_squidpy_arguments(): + adata = AnnData(np.ones((2, 2), dtype=np.float32)) + + with pytest.raises(ValueError, match="'resolutions' is required for flavor 'utag'"): + squidpy_backend.calculate_niche(adata, flavor="utag", n_neighbors=10) + + +def test_ligrec_drops_cpu_only_options(monkeypatch): + adata = AnnData(np.ones((2, 2), dtype=np.float32)) + captured = {} + + def fake_ligrec(data, cluster_key, **kwargs): + captured.update(kwargs) + + monkeypatch.setattr(squidpy_backend, "_ligrec", fake_ligrec) + + with pytest.warns(UserWarning, match="seed, n_jobs have no effect"): + squidpy_backend.ligrec(adata, "cluster", seed=0, n_jobs=4, n_perms=5) + + assert captured == {"n_perms": 5} + assert inspect.signature(squidpy_backend.ligrec) == inspect.signature(rsc.gr.ligrec) + + +def test_squidpy_backend_dispatch_smoke(monkeypatch): + squidpy_backends = pytest.importorskip("squidpy._backends") + sq = pytest.importorskip("squidpy") + + registry = squidpy_backends.dispatcher._registry + dispatch_impl = squidpy_backends.dispatcher._dispatch_impl + old_backend = squidpy_backends.settings.backend + old_state = { + "_backends": copy(registry._backends), + "_alias_map": copy(registry._alias_map), + "_load_errors": copy(registry._load_errors), + "_registration_errors": copy(registry._registration_errors), + "_warned_untrusted": copy(registry._warned_untrusted), + "_discovered": registry._discovered, + "_sig_cache": copy(dispatch_impl._sig_cache), + } + + @wraps(squidpy_backend._utag) + def fake_utag(data, **kwargs): + data.uns["squidpy_backend_called"] = kwargs + return None if kwargs["inplace"] else data + + monkeypatch.setattr(squidpy_backend, "_utag", fake_utag) + + try: + squidpy_backends.settings._backend_var.set("cpu") + registry._backends.clear() + registry._alias_map.clear() + registry._load_errors.clear() + registry._registration_errors.clear() + registry._warned_untrusted.clear() + registry._discovered = True + registry._register_backend( + squidpy_backend, + entrypoint_name="rapids-singlecell", + distribution_name="rapids-singlecell", + object_ref="rapids_singlecell._backends.squidpy", + ) + dispatch_impl._sig_cache.clear() + dispatch_impl._update_signatures() + + adata = AnnData(np.ones((2, 2), dtype=np.float32)) + result = sq.gr.calculate_niche( + adata, + flavor="utag", + n_neighbors=10, + resolutions=0.5, + backend="cuda", + inplace=False, + ) + + assert result is adata + assert adata.uns["squidpy_backend_called"]["inplace"] is False + finally: + squidpy_backends.settings._backend_var.set(old_backend) + registry._backends.clear() + registry._backends.update(old_state["_backends"]) + registry._alias_map.clear() + registry._alias_map.update(old_state["_alias_map"]) + registry._load_errors.clear() + registry._load_errors.update(old_state["_load_errors"]) + registry._registration_errors.clear() + registry._registration_errors.update(old_state["_registration_errors"]) + registry._warned_untrusted.clear() + registry._warned_untrusted.update(old_state["_warned_untrusted"]) + registry._discovered = old_state["_discovered"] + dispatch_impl._sig_cache.clear() + dispatch_impl._sig_cache.update(old_state["_sig_cache"]) + dispatch_impl._update_signatures()