diff --git a/docs/api.md b/docs/api.md index e86551a..0d35936 100644 --- a/docs/api.md +++ b/docs/api.md @@ -74,7 +74,7 @@ the infercnvpy equivalents operate on CNV data. ### Spatial -Tools for spatial data. They read a cells x cells neighbour graph from `adata.obsp`, +Tools for spatial data. They read a cells x cells neighbor graph from `adata.obsp`, e.g. `obsp["spatial_connectivities"]` from {func}`squidpy.gr.spatial_neighbors_knn`. ```{eval-rst} diff --git a/docs/installation.md b/docs/installation.md index b71987d..64bd3a2 100644 --- a/docs/installation.md +++ b/docs/installation.md @@ -57,7 +57,7 @@ uv add git+https://github.com/scverse/infercnvpy.git - `infercnvpy[dask]` installs dask and the lazy readers of anndata (`anndata[dask,lazy]`) for {doc}`large and lazy datasets `. - `infercnvpy[gtf]` installs [gtfparse][] for {func}`infercnvpy.io.genomic_position_from_gtf`. -- `infercnvpy[spatial]` installs [squidpy][] for the spatial neighbour graph that {func}`infercnvpy.tl.smooth_cnv`, +- `infercnvpy[spatial]` installs [squidpy][] for the spatial neighbor graph that {func}`infercnvpy.tl.smooth_cnv`, {func}`infercnvpy.tl.spatial_regions` and {func}`infercnvpy.pl.spatial_cnv` read from `adata.obsp`. The functions themselves work with any graph in `adata.obsp`, so squidpy is only needed to build it. - `infercnvpy[copykat]` installs [rpy2][] for {func}`infercnvpy.tl.copykat`. diff --git a/docs/large_data.md b/docs/large_data.md index 61e85ee..f286983 100644 --- a/docs/large_data.md +++ b/docs/large_data.md @@ -9,7 +9,9 @@ as a dense array. The expression matrix can be - a [dask](https://www.dask.org/) array, also with sparse `scipy` blocks, e.g. a matrix loaded lazily with {func}`anndata.experimental.read_elem_lazy` or {func}`anndata.experimental.read_lazy`, or a table from [spatialdata-io](https://spatialdata.scverse.org/projects/io/), -- the matrix of a backed AnnData (`anndata.read_h5ad(..., backed="r")`). +- the matrix of a backed AnnData (`anndata.read_h5ad(..., backed="r")`), stored as CSR or dense. Reading a block + of rows of a CSC matrix reads all of it, so a backed CSC matrix is read into memory once, as CSR, with a warning. + Convert it before writing the file (`adata.X = adata.X.tocsr()`) to read it block by block. Install the optional dependency with `pip install "infercnvpy[dask]"` to use dask arrays. @@ -41,12 +43,12 @@ With `lazy=True`, the result for a dask array stays a dask array of CSR blocks, {meth}`dask.array.Array.map_blocks`. The noise band is computed first, with one dask computation. `lazy=True` has to be passed explicitly, also when a {class}`distributed.Client` is active. In this mode the dask scheduler decides how many blocks run at the same time, each in a single thread, and `n_jobs` is ignored. infercnvpy never starts numba -threads inside a dask worker thread, because numba's default threading layer on macOS aborts the process when it is -used from several threads at once. +threads inside a dask task, with any scheduler: the scheduler's threads or processes already use the cores, and +numba's default threading layer on macOS aborts the process when it is used from several threads at once. :::{note} -infercnv expects log-normalised expression (see {ref}`input-data`). For a dask array, {func}`scanpy.pp.normalize_total` -and {func}`scanpy.pp.log1p` stay lazy. A backed AnnData cannot be normalised in place by scanpy, so normalise the data +infercnv expects log-normalized expression (see {ref}`input-data`). For a dask array, {func}`scanpy.pp.normalize_total` +and {func}`scanpy.pp.log1p` stay lazy. A backed AnnData cannot be normalized in place by scanpy, so normalize the data before you write the file. ::: @@ -54,7 +56,7 @@ before you write the file. Spatial technologies often produce tables with a lazy dask `X` of CSR blocks that are indexed by Ensembl ID. Next to the genes, such a table can contain control probes or codewords without a genomic position, which you remove before -normalisation. +normalization. ```python import anndata as ad diff --git a/docs/notebooks/reproduce_infercnv.ipynb b/docs/notebooks/reproduce_infercnv.ipynb index 98b4fcb..3b581e0 100644 --- a/docs/notebooks/reproduce_infercnv.ipynb +++ b/docs/notebooks/reproduce_infercnv.ipynb @@ -202,7 +202,7 @@ "id": "artificial-animation", "metadata": {}, "source": [ - "It contains malignant cells from four tumours and two types of non-malignant cells." + "It contains malignant cells from four tumors and two types of non-malignant cells." ] }, { diff --git a/docs/notebooks/spatial_xenium.ipynb b/docs/notebooks/spatial_xenium.ipynb index e5c6f88..53a48e7 100644 --- a/docs/notebooks/spatial_xenium.ipynb +++ b/docs/notebooks/spatial_xenium.ipynb @@ -9,22 +9,22 @@ "\n", "Imaging-based spatial transcriptomics (Xenium, CosMx, MERSCOPE) measures a few thousand genes in hundreds of\n", "thousands of cells, with a few hundred counts per cell. inferCNV-style profiles of single cells are noisy at that\n", - "depth, but the tissue coordinates add information: neighbouring tumour cells usually share their copy number state.\n", + "depth, but the tissue coordinates add information: neighboring tumor cells usually share their copy number state.\n", "This tutorial shows the spatial layer of infercnvpy:\n", "\n", - "- {func}`infercnvpy.tl.smooth_cnv` averages each cell's CNV profile with its spatial neighbours,\n", - "- {func}`infercnvpy.tl.spatial_regions` makes per-cell labels (here tumour/normal calls) spatially coherent: a cell\n", - " takes the label of its neighbourhood when most of its neighbours disagree with it (more than 75% at the default\n", + "- {func}`infercnvpy.tl.smooth_cnv` averages each cell's CNV profile with its spatial neighbors,\n", + "- {func}`infercnvpy.tl.spatial_regions` makes per-cell labels (here tumor/normal calls) spatially coherent: a cell\n", + " takes the label of its neighborhood when most of its neighbors disagree with it (more than 75% at the default\n", " `beta`),\n", "- {func}`infercnvpy.pl.spatial_cnv` plots the mean CNV of a chromosome or arm on the tissue.\n", "\n", - "All three read a spatial neighbour graph from `adata.obsp`, here from {func}`squidpy.gr.spatial_neighbors_knn`; any\n", + "All three read a spatial neighbor graph from `adata.obsp`, here from {func}`squidpy.gr.spatial_neighbors_knn`; any\n", "cells x cells graph works. To run this tutorial, install squidpy, gtfparse and the lazy readers of anndata with the optional extras, and s3fs to\n", "read the data from S3: `pip install \"infercnvpy[spatial,gtf,dask]\" s3fs`.\n", "\n", - "The ideas are not new. InSituCNV {cite}`Jensen2025` smooths infercnvpy profiles over kNN neighbours on Xenium and CosMx,\n", - "SlideCNA {cite}`Zhang2025` pools spatial neighbours on Slide-seq, and STARCH {cite}`Elyanow2021` and CalicoST\n", - "{cite}`Ma2024` make neighbouring cells prefer the same clone label. In a benchmark of CNV inference on spatial data,\n", + "The ideas are not new. InSituCNV {cite}`Jensen2025` smooths infercnvpy profiles over kNN neighbors on Xenium and CosMx,\n", + "SlideCNA {cite}`Zhang2025` pools spatial neighbors on Slide-seq, and STARCH {cite}`Elyanow2021` and CalicoST\n", + "{cite}`Ma2024` make neighboring cells prefer the same clone label. In a benchmark of CNV inference on spatial data,\n", "{cite:t}`Han2026` found that aggregating cells on imaging platforms improved CNV recovery compared to native\n", "single-cell resolution.\n", "\n", @@ -122,8 +122,8 @@ "id": "2f7e8f19", "metadata": {}, "source": [ - "We drop cells with fewer than 50 counts and log-normalise to 10,000 counts per cell. The CNV values are log fold\n", - "changes, so the normalisation target sets their scale: with the default target (the median of the remaining cells),\n", + "We drop cells with fewer than 50 counts and log-normalize to 10,000 counts per cell. The CNV values are log fold\n", + "changes, so the normalization target sets their scale: with the default target (the median of the remaining cells),\n", "most values are `log1p` of 0 or 1 count, which compresses the CNV values." ] }, @@ -183,7 +183,7 @@ "expression data and annotate the clusters with marker gene scores. The Prime 5K panel leaves out many classic,\n", "highly expressed markers, e.g. PTPRC, KRT5 and COL1A1, which are missing from the filtered lists printed below. We\n", "use T, B and myeloid cell markers for immune cells, endothelial markers, fibroblast markers and squamous epithelial\n", - "markers, as far as they are in the panel. Immune cells, endothelial cells and fibroblasts are the reference; the epithelial cells contain the tumour." + "markers, as far as they are in the panel. Immune cells, endothelial cells and fibroblasts are the reference; the epithelial cells contain the tumor." ] }, { @@ -506,8 +506,8 @@ "so we use running windows of 51 genes, every 5 genes, instead of the defaults\n", "(100 and 10). We turn off the noise filter (`dynamic_threshold=None`), because it sets small values to 0 and\n", "{func}`infercnvpy.tl.smooth_cnv` would then average those zeros. Then we build a spatial graph of the 10 nearest\n", - "neighbours of each cell with squidpy and average each cell's CNV profile with its neighbours. By default, the cell\n", - "itself counts as one more neighbour." + "neighbors of each cell with squidpy and average each cell's CNV profile with its neighbors. By default, the cell\n", + "itself counts as one more neighbor." ] }, { @@ -674,11 +674,11 @@ "id": "e088f7f9", "metadata": {}, "source": [ - "## Tumour calls\n", + "## Tumor calls\n", "\n", - "We score each cell by the mean absolute value of its smoothed CNV profile and call it tumour when the score is above\n", + "We score each cell by the mean absolute value of its smoothed CNV profile and call it tumor when the score is above\n", "the 95th percentile of the scores of the reference cells, so the score and the threshold are both per cell. By\n", - "construction, 5% of the reference cells are called tumour, and normal cells of types outside the reference can exceed\n", + "construction, 5% of the reference cells are called tumor, and normal cells of types outside the reference can exceed\n", "the threshold as well." ] }, @@ -790,8 +790,8 @@ "id": "8e291530", "metadata": {}, "source": [ - "Most epithelial cells are called tumour. The fibroblasts, immune and endothelial cells among the tumour calls (share\n", - "printed above) are the 5% of reference cells above the threshold. The tumour calls show the loss of 3p, the gain of\n", + "Most epithelial cells are called tumor. The fibroblasts, immune and endothelial cells among the tumor calls (share\n", + "printed above) are the 5% of reference cells above the threshold. The tumor calls show the loss of 3p, the gain of\n", "3q and losses on chromosomes 4, 5 and 9." ] }, @@ -835,8 +835,8 @@ "source": [ "## Spatially coherent regions\n", "\n", - "{func}`infercnvpy.tl.spatial_regions` makes labels spatially coherent: a cell takes the label of its neighbourhood\n", - "when enough neighbours disagree with it (for two labels and the default `beta=2`, more than 75% of them).\n", + "{func}`infercnvpy.tl.spatial_regions` makes labels spatially coherent: a cell takes the label of its neighborhood\n", + "when enough neighbors disagree with it (for two labels and the default `beta=2`, more than 75% of them).\n", "The calls from the smoothed profiles are already spatially coherent, and few cells change." ] }, @@ -930,9 +930,9 @@ "source": [ "The prior matters more for noisy labels. Calls of single cells from the unsmoothed CNV profiles, with the same rule,\n", "are scattered over the section (left in the figure below). With `beta=4` (a cell changes its label when more than\n", - "62.5% of its neighbours disagree), `tl.spatial_regions` removes about half of the tumour calls (8,788 to 4,710 in the\n", - "table below), including 89% of the calls of non-epithelial cells (middle). It doesn't recover tumour cells that the\n", - "unsmoothed profiles miss, so the regularised calls from the smoothed profiles (right) include more than twice as many\n", + "62.5% of its neighbors disagree), `tl.spatial_regions` removes about half of the tumor calls (8,788 to 4,710 in the\n", + "table below), including 89% of the calls of non-epithelial cells (middle). It doesn't recover tumor cells that the\n", + "unsmoothed profiles miss, so the regularized calls from the smoothed profiles (right) include more than twice as many\n", "epithelial cells." ] }, @@ -1064,13 +1064,13 @@ "source": [ "## Notes\n", "\n", - "- With 10 neighbours and the default weight of the cell itself, a cell's own profile is 1/11 of its smoothed profile, so calls from\n", - " smoothed profiles are calls of neighbourhoods. Smoothing mixes the profiles of tumour and normal cells at the\n", - " tumour boundary and of immune cells that infiltrate the tumour, over about `n_hops` cells. Use the unsmoothed\n", + "- With 10 neighbors and the default weight of the cell itself, a cell's own profile is 1/11 of its smoothed profile, so calls from\n", + " smoothed profiles are calls of neighborhoods. Smoothing mixes the profiles of tumor and normal cells at the\n", + " tumor boundary and of immune cells that infiltrate the tumor, over about `n_hops` cells. Use the unsmoothed\n", " `X_cnv` for per-cell questions in mixed regions.\n", - "- The spatial graph defines the neighbourhood. A kNN graph adapts to the cell density; a radius graph\n", + "- The spatial graph defines the neighborhood. A kNN graph adapts to the cell density; a radius graph\n", " ({func}`squidpy.gr.spatial_neighbors_radius`) keeps the physical scale fixed. `n_hops=2` in `tl.smooth_cnv`\n", - " widens the neighbourhood further.\n", + " widens the neighborhood further.\n", "- `tl.spatial_regions` changes labels only as long as a single change makes them more coherent, so the result stays\n", " close to the input labels and is not a new clustering. Larger `beta` gives smoother regions." ] diff --git a/src/infercnvpy/_numba.py b/src/infercnvpy/_numba.py index c267102..bd8e0a9 100644 --- a/src/infercnvpy/_numba.py +++ b/src/infercnvpy/_numba.py @@ -1,10 +1,12 @@ """numba helpers: a parallel kernel with a serial twin, and `n_jobs` as the number of numba threads. -A ``parallel=True`` kernel only runs on the main thread. Called from any other thread (a dask worker, a thread pool), -the serial twin runs instead: numba's default ``workqueue`` threading layer (the only one in the macOS wheels) aborts -the process on concurrent use, and a caller that brings its own threads should not get numba threads on top. +A ``parallel=True`` kernel only runs on the main thread and outside dask tasks (`dask_task`). Anywhere else (a dask +worker, a thread pool, the main thread of a dask worker process) the serial twin runs instead: numba's default +``workqueue`` threading layer (the only one in the macOS wheels) aborts the process on concurrent use, and a caller +that brings its own threads or processes should not get numba threads on top. -All kernels use numba's disk cache, with data files named after their cache key (see `_KeyedCacheFile`). +All kernels use numba's disk cache, with data files named after their cache key (see `_KeyedCacheFile`), and numpy's +error model, so that both twins divide by zero alike. """ from __future__ import annotations @@ -12,7 +14,7 @@ import hashlib import threading from contextlib import contextmanager -from functools import update_wrapper +from functools import update_wrapper, wraps from types import FunctionType from typing import TYPE_CHECKING, Any @@ -25,7 +27,7 @@ else: from numba import prange -__all__ = ["n_kernel_chunks", "njit", "njit_parallel", "numba_threads", "prange"] +__all__ = ["dask_task", "n_kernel_chunks", "njit", "njit_parallel", "numba_threads", "prange", "runs_serially"] class _KeyedCacheFile(IndexDataCacheFile): @@ -58,18 +60,18 @@ def __init__(self, py_func): def njit(fn: Callable[..., Any], /, **options: Any) -> Any: - """`numba.njit` with ``nogil=True`` and the disk cache of `_KeyedCacheFile`.""" - dispatcher = numba.njit(fn, nogil=True, **options) - dispatcher._cache = _KeyedFunctionCache(dispatcher.py_func) # type: ignore[attr-defined] # what `cache=True` sets + """`numba.njit` with ``nogil=True``, numpy's error model and the disk cache of `_KeyedCacheFile`.""" + dispatcher = numba.njit(fn, nogil=True, error_model="numpy", **options) + if hasattr(dispatcher, "py_func"): # with NUMBA_DISABLE_JIT=1, `njit` returns `fn` itself + dispatcher._cache = _KeyedFunctionCache(dispatcher.py_func) # type: ignore[attr-defined] # what `cache=True` sets return dispatcher class njit_parallel: """Compile a function with `njit` twice, ``parallel=True`` and serial. - Calling it runs the parallel version on the main thread and the serial one on any other thread, or always the - serial one with ``serial=True`` (for code that runs inside dask tasks). The two versions are distinct function - objects with distinct names, so they don't share a cache entry. + Calling it runs the serial version where `runs_serially` says so, otherwise the parallel one. The two versions are + distinct function objects with distinct names, so they don't share a cache entry. """ def __init__(self, fn: Callable[..., Any], /): @@ -78,12 +80,31 @@ def __init__(self, fn: Callable[..., Any], /): self.serial = njit(_renamed(fn, "serial")) update_wrapper(self, fn) - def __call__(self, *args, serial: bool = False): - return (self.serial if serial or not _on_main_thread() else self.parallel)(*args) + def __call__(self, *args): + return (self.serial if runs_serially() else self.parallel)(*args) -def _on_main_thread() -> bool: - return threading.current_thread() is threading.main_thread() +_scope = threading.local() + + +def runs_serially() -> bool: + """Whether kernels called here run serially: off the main thread, or inside a dask task (see `dask_task`).""" + return threading.current_thread() is not threading.main_thread() or getattr(_scope, "in_dask_task", False) + + +def dask_task[F: Callable[..., Any]](fn: F) -> F: + """Mark `fn` as a dask task: the kernels it calls run serially, whatever thread or process the scheduler uses.""" + + @wraps(fn) + def task(*args, **kwargs): + outer = getattr(_scope, "in_dask_task", False) + _scope.in_dask_task = True + try: + return fn(*args, **kwargs) + finally: + _scope.in_dask_task = outer + + return task # type: ignore[return-value] def _renamed(fn: FunctionType, suffix: str) -> FunctionType: @@ -93,9 +114,9 @@ def _renamed(fn: FunctionType, suffix: str) -> FunctionType: return new -def n_kernel_chunks(n_rows: int, *, serial: bool = False, per_thread: int = 4) -> int: +def n_kernel_chunks(n_rows: int, *, per_thread: int = 4) -> int: """Number of row chunks for a `prange` kernel: `per_thread` per numba thread, or 1 if it runs serially.""" - if serial or not _on_main_thread(): + if runs_serially(): return 1 return max(1, min(n_rows, per_thread * numba.get_num_threads())) diff --git a/src/infercnvpy/_util.py b/src/infercnvpy/_util.py index 4f73a66..fcbce9f 100644 --- a/src/infercnvpy/_util.py +++ b/src/infercnvpy/_util.py @@ -1,8 +1,10 @@ import inspect +import re import sys import warnings from collections.abc import Callable from functools import wraps +from textwrap import dedent, indent from typing import Any import numpy as np @@ -16,8 +18,10 @@ def _as_matrix(x, what: str = "the expression matrix", *, lazy: bool = False): """Narrow an AnnData array to an in-memory numpy array or scipy sparse matrix. With `lazy`, also accept a dask array or a backed (on-disk) matrix, see `_is_dask` and `_is_backed`. Without it, - such a matrix raises an error that tells the user how to load it. + such a matrix raises an error that tells the user how to load it. A DataFrame (as obsm allows) gives its values. """ + if isinstance(x, pd.DataFrame): + return x.to_numpy() if isinstance(x, np.ndarray) or sp.issparse(x) or (lazy and (_is_dask(x) or _is_backed(x))): return x if _is_dask(x) or _is_backed(x): @@ -86,6 +90,21 @@ def writable(v): return {k: writable(v) for k, v in params.items()} +def _splice_doc[F: Callable[..., Any]](**blocks: str) -> Callable[[F], F]: + """Replace each line `{name}` of the docstring by the text `blocks[name]`, at the indentation of that line.""" + + def decorate(fn: F) -> F: + doc = fn.__doc__ or "" + for name, block in blocks.items(): + match = re.search(rf"^( *)\{{{name}\}}$", doc, flags=re.MULTILINE) + assert match is not None, f"{fn.__name__} has no line {{{name}}}" + doc = doc.replace(match.group(0), indent(dedent(block).strip("\n"), match.group(1))) + fn.__doc__ = doc + return fn + + return decorate + + def _deprecated_positional[F: Callable[..., Any]](*names: str, until: str = "0.8.0") -> Callable[[F], F]: """Accept the keyword-only arguments `names` positionally, in this order, with a FutureWarning. diff --git a/src/infercnvpy/datasets/__init__.py b/src/infercnvpy/datasets/__init__.py index 183579a..d215298 100644 --- a/src/infercnvpy/datasets/__init__.py +++ b/src/infercnvpy/datasets/__init__.py @@ -5,6 +5,7 @@ from textwrap import dedent, indent import pooch +import scanpy as sc from anndata import AnnData, read_h5ad from scverse_misc.datasets import DatasetEntry, DownloadCB, fetch, parse_registry, register_loader @@ -16,13 +17,17 @@ # scverse-misc keys loaders by type name across all packages in the process, so ours uses the package name. @register_loader("infercnvpy") def _load_h5ad(entry: DatasetEntry, target: Path, download: DownloadCB, /, **kwargs: object) -> AnnData: - # Download into the cache dir itself rather than `target` (`/`), - # so files cached by earlier infercnvpy versions are reused. + # Download into the cache dir itself rather than `target` (`/`): `_fetch` passes scanpy's + # `datasetdir` as the cache dir when infercnvpy 0.6 downloaded the file there, so that file is reused. return read_h5ad(download(entry.file(suffix=".h5ad"), dest=target.parent)) def _fetch(name: str) -> AnnData: - cache_dir = os.environ.get(DATASET_ENV_VAR) or pooch.os_cache("infercnvpy") + cache_dir: str | Path | None = os.environ.get(DATASET_ENV_VAR) + if not cache_dir: + (file,) = _DATASETS[name].files + old = Path(sc.settings.datasetdir) # where infercnvpy 0.6 downloaded maynard2020_3k + cache_dir = old if (old / file.name).is_file() else pooch.os_cache("infercnvpy") return fetch(_DATASETS[name], cache_dir, base_url=_BASE_URL) @@ -33,7 +38,8 @@ def _fetch(name: str) -> AnnData: By default, the dataset will be downloaded into your operating system's default cache directory (See :func:`pooch.os_cache` for more details). If it has already been - downloaded, it will be retrieved from the cache. + downloaded, it will be retrieved from the cache, or from :attr:`scanpy.settings.datasetdir` + if an earlier version of infercnvpy downloaded it there. You can override the default cache dir by setting the `{DATASET_ENV_VAR}` environment variable to a path of your preference. diff --git a/src/infercnvpy/io/_genepos.py b/src/infercnvpy/io/_genepos.py index c9758e6..695d79f 100644 --- a/src/infercnvpy/io/_genepos.py +++ b/src/infercnvpy/io/_genepos.py @@ -74,21 +74,7 @@ def genomic_position_from_biomart( remedy=f"check that the gene ids in adata are {biomart_gene_id} values (see `adata_gene_id`)", ) - tmp_var = _as_df(adata.var).copy() - orig_index_name = tmp_var.index.name - TMP_INDEX_NAME = "adata_var_index" - tmp_var.index.name = TMP_INDEX_NAME - tmp_var.reset_index(inplace=True) - var_annotated = tmp_var.merge( - biomart_annot, - how="left", - left_on=TMP_INDEX_NAME if adata_gene_id is None else adata_gene_id, - right_on=biomart_gene_id, - validate="one_to_one", - ) - var_annotated.set_index(TMP_INDEX_NAME, inplace=True) - var_annotated.index.name = orig_index_name - + var_annotated = _merge_var(adata, biomart_annot, adata_gene_id, biomart_gene_id) if inplace: adata.var = var_annotated return None @@ -133,6 +119,8 @@ def genomic_position_from_gtf( 'genomic_position_from_gtf needs the optional dependency gtfparse. Install it with `pip install "infercnvpy[gtf]"`.' ) from None gtf = gtfparse.read_gtf(gtf_file, usecols=["seqname", "feature", "start", "end", "gene_id", "gene_name"]) + if hasattr(gtf, "to_pandas"): # gtfparse 2 returns a polars DataFrame + gtf = gtf.to_pandas() gtf = ( gtf.loc[ gtf["feature"] == "gene", @@ -157,26 +145,29 @@ def genomic_position_from_gtf( remedy="check `gtf_gene_id` and `adata_gene_id`, and that the GTF file matches the genome build of the data", ) + var_annotated = _add_chr_prefix(_merge_var(adata, gtf, adata_gene_id, gtf_gene_id)) + if inplace: + adata.var = var_annotated + return None + return var_annotated + + +def _merge_var(adata: AnnData, annot: pd.DataFrame, adata_gene_id: str | None, annot_gene_id: str) -> pd.DataFrame: + """`adata.var` with the columns of `annot` whose `annot_gene_id` matches `adata_gene_id` (default: var_names).""" tmp_var = _as_df(adata.var).copy() orig_index_name = tmp_var.index.name TMP_INDEX_NAME = "adata_var_index" tmp_var.index.name = TMP_INDEX_NAME tmp_var.reset_index(inplace=True) var_annotated = tmp_var.merge( - gtf, + annot, how="left", left_on=TMP_INDEX_NAME if adata_gene_id is None else adata_gene_id, - right_on=gtf_gene_id, + right_on=annot_gene_id, validate="one_to_one", ) var_annotated.set_index(TMP_INDEX_NAME, inplace=True) var_annotated.index.name = orig_index_name - - var_annotated = _add_chr_prefix(var_annotated) - - if inplace: - adata.var = var_annotated - return None return var_annotated diff --git a/src/infercnvpy/io/_scevan.py b/src/infercnvpy/io/_scevan.py index fe615fb..7b1b068 100644 --- a/src/infercnvpy/io/_scevan.py +++ b/src/infercnvpy/io/_scevan.py @@ -96,7 +96,7 @@ def read_scevan( Sets the following fields: `adata.obs["{key_added}_class"]`, `adata.obs["{key_added}_confident_normal"]` - Only with `scevan_res_table`: tumour/normal class of each cell and whether it is a + Only with `scevan_res_table`: tumor/normal class of each cell and whether it is a confident normal cell. `adata.obs["{key_added}_subclone"]` Only with `scevan_res_table` from a run with subclone calling: the subclone of each cell. diff --git a/src/infercnvpy/pl/_chromosome_heatmap.py b/src/infercnvpy/pl/_chromosome_heatmap.py index 045cd62..ab0e880 100644 --- a/src/infercnvpy/pl/_chromosome_heatmap.py +++ b/src/infercnvpy/pl/_chromosome_heatmap.py @@ -1,7 +1,7 @@ from __future__ import annotations import warnings -from dataclasses import dataclass +from dataclasses import dataclass, fields from pathlib import Path from typing import TYPE_CHECKING, Any, Literal @@ -21,7 +21,8 @@ from infercnvpy._logging import SKIP_OWN_FRAMES from infercnvpy._numba import n_kernel_chunks, numba_threads -from infercnvpy._util import _as_matrix, _is_backed, _is_dask, _obs_column +from infercnvpy._util import _as_matrix, _is_backed, _is_dask, _obs_column, _splice_doc +from infercnvpy.pl._chromosome_heatmap_legacy import legacy_heatmap, legacy_summary from infercnvpy.pl._heatmap_kernels import binned_sums, column_ranges, gram, pixel_rows, sampled_rows if TYPE_CHECKING: @@ -32,7 +33,7 @@ from matplotlib.figure import Figure # Sizes in inches, as `scanpy.pl.heatmap` lays them out at figsize=(16, 10). -_GROUP_W = 0.155 # colour blocks of the groups +_GROUP_W = 0.155 # color blocks of the groups _CBAR_W = 0.155 _DENDRO_W = 0.775 _CHR_H = 0.115 # chromosome brackets @@ -40,16 +41,59 @@ _PAD_H = 0.05 _CBAR_MAX_H = 4.0 -#: Colour bar label: infercnv, smooth_cnv and SCEVAN all report CNV relative to the reference cells. +#: Color bar label: infercnv, smooth_cnv and SCEVAN all report CNV relative to the reference cells. CBAR_LABEL = "CNV (relative to reference)" -#: :func:`scanpy.pl.heatmap` arguments that `mode="fast"` implements; any other one switches to `mode="legacy"`. -FAST_KWARGS = frozenset( - {"swap_axes", "standard_scale", "var_group_rotation", "dendrogram", "vmin", "vmax", "vcenter", "norm"} - | {"show_gene_labels", "log"} -) - +#: Docstring blocks shared by `chromosome_heatmap` and `chromosome_heatmap_summary`. +_HEATMAP_DOC = { + "data": """ + adata + annotated data matrix + groupby + group the cells by a categorical variable from adata.obs. It usually makes + sense to either group by unsupervised clustering obtained from + :func:`infercnvpy.tl.leiden` (the default) or a cell-type label. + use_rep + Key under which the result from :func:`infercnvpy.tl.infercnv` are stored. With `mode="fast"`, + `adata.obsm[f"X_{use_rep}"]` can also be a dask array or a backed matrix, which is read block by block. + cmap + colormap to use + figsize + (width, height) tuple in inches. Ignored if `ax` is given. + """, + "fast_args": """ + `mode="fast"` implements the :func:`scanpy.pl.heatmap` arguments `swap_axes`, `standard_scale`, + `var_group_rotation`, `dendrogram`, `vmin`, `vmax`, `vcenter`, `norm`, `show_gene_labels` (ignored) and + `log`; with any other argument (or a numeric, not categorical, `groupby` column), it warns and uses + `mode="legacy"`, or raises an error if the CNV matrix is not in memory. + """, + "figure": """ + ax + Draw the heatmap into this axes (`mode="fast"` only). The group colors, chromosome labels, dendrogram and + color bar go into axes split off from it. + n_jobs + Number of threads for reading the CNV matrix (`mode="fast"`). `None` uses all cores. + show + Whether to show the figure. Defaults to :attr:`scanpy.settings.autoshow`. + save + If `True` or a `str`, save the figure to :attr:`scanpy.settings.figdir`. A string is appended to the default + filename `heatmap`; it sets the file type if it ends in `.pdf`, `.png` or `.svg`, otherwise + :attr:`scanpy.settings.file_format_figs` does. + return_fig + Return the :class:`~matplotlib.figure.Figure`. + """, + "returns": """ + Returns + ------- + The figure if `return_fig` is true. Otherwise, if `show` is false, a dictionary of axes with the keys + `heatmap_ax`, `groupby_ax` (group colors), `gene_groups_ax` (chromosome labels), `colorbar_ax` (not with + `mode="legacy"`) and, with `dendrogram`, `dendrogram_ax`. + """, +} + + +@_splice_doc(**_HEATMAP_DOC) def chromosome_heatmap( adata: AnnData, *, @@ -79,64 +123,33 @@ def chromosome_heatmap( Parameters ---------- - adata - annotated data matrix - groupby - group the cells by a categorical variable from adata.obs. It usually makes - sense to either group by unsupervised clustering obtained from - :func:`infercnvpy.tl.leiden` (the default) or a cell-type label. - use_rep - Key under which the result from :func:`infercnvpy.tl.infercnv` are stored. With `mode="fast"`, - `adata.obsm[f"X_{use_rep}"]` can also be a dask array or a backed matrix, which is read block by block. - cmap - colormap to use - figsize - (width, height) tuple in inches. Ignored if `ax` is given. + {data} mode `"fast"` draws the heatmap with infercnvpy at the resolution of the figure. `"legacy"` draws every cell with :func:`scanpy.pl.heatmap`, which needs the CNV matrix as a dense table and several times its memory. - `mode="fast"` implements the :func:`scanpy.pl.heatmap` arguments `swap_axes`, `standard_scale`, - `var_group_rotation`, `dendrogram`, `vmin`, `vmax`, `vcenter`, `norm`, `show_gene_labels` (ignored) and - `log`; with any other argument (or a numeric, not categorical, `groupby` column), it warns and uses - `mode="legacy"`, or raises an error if the CNV matrix is not in memory. + {fast_args} aggregate What a pixel row shows when there are more cells than pixel rows (`mode="fast"` only). `"sample"` shows the cell in the middle of the cells that fall into the row, as the full-resolution plot downscaled with - nearest-neighbour interpolation, with the colour scale from the values of all cells. `"mean"` shows their - mean, with the colour scale from the range of the drawn (averaged) values unless `vmin` and `vmax` are + nearest-neighbor interpolation, with the color scale from the values of all cells. `"mean"` shows their + mean, with the color scale from the range of the drawn (averaged) values unless `vmin` and `vmax` are given. - ax - Draw the heatmap into this axes (`mode="fast"` only). The group colours, chromosome labels, dendrogram and - colour bar go into axes split off from it. - n_jobs - Number of threads for reading the CNV matrix (`mode="fast"`). `None` uses all cores. - show - Whether to show the figure. Defaults to :attr:`scanpy.settings.autoshow`. - save - If `True` or a `str`, save the figure to :attr:`scanpy.settings.figdir`. A string is appended to the default - filename `heatmap`; it sets the file type if it ends in `.pdf`, `.png` or `.svg`, otherwise - :attr:`scanpy.settings.file_format_figs` does. - return_fig - Return the :class:`~matplotlib.figure.Figure`. + {figure} **kwargs - Arguments of :func:`scanpy.pl.heatmap`, see `mode`. By default, the colour scale is a - :class:`~matplotlib.colors.TwoSlopeNorm` centred at `vcenter=0`, with limits from the data. Pass the same + Arguments of :func:`scanpy.pl.heatmap`, see `mode`. By default, the color scale is a + :class:`~matplotlib.colors.TwoSlopeNorm` centered at `vcenter=0`, with limits from the data. Pass the same `vmin` and `vmax` to compare figures across modes, representations or datasets. `dendrogram=True` orders the groups as :func:`scanpy.tl.dendrogram` with its defaults clusters them (as `mode="legacy"`; with more windows than `scanpy.settings.N_PCS`, it clusters the groups' means of a PCA of all cells) and draws the dendrogram next to the heatmap. A dendrogram stored in `adata.uns[f"dendrogram_{groupby}"]` is used if it has the same groups; otherwise the computed one is stored there, as `mode="legacy"` does. - Returns - ------- - The figure if `return_fig` is true. Otherwise, if `show` is false, a dictionary of axes with the keys - `heatmap_ax`, `groupby_ax` (group colours), `gene_groups_ax` (chromosome labels), `colorbar_ax` (not with - `mode="legacy"`) and, with `dendrogram`, `dendrogram_ax`. + {returns} """ if aggregate not in ("sample", "mean"): raise ValueError(f"`aggregate` must be 'sample' or 'mean', got {aggregate!r}.") if _use_legacy("chromosome_heatmap", adata, groupby, use_rep, mode, kwargs): - return _legacy(_legacy_heatmap, adata, groupby, use_rep, cmap, figsize, ax, show, save, return_fig, kwargs) + return _legacy(legacy_heatmap, adata, groupby, use_rep, cmap, figsize, ax, show, save, return_fig, kwargs) opts = _Options.pop(kwargs) own = ax is None fig, ax = _figure(ax, figsize) @@ -167,6 +180,7 @@ def chromosome_heatmap( return _finish(fig, axes, show=show, save=save, return_fig=return_fig, own=own) +@_splice_doc(**_HEATMAP_DOC) def chromosome_heatmap_summary( adata: AnnData, *, @@ -192,50 +206,20 @@ def chromosome_heatmap_summary( Parameters ---------- - adata - annotated data matrix - groupby - group the cells by a categorical variable from adata.obs. It usually makes - sense to either group by unsupervised clustering obtained from - :func:`infercnvpy.tl.leiden` (the default) or a cell-type label. - use_rep - Key under which the result from :func:`infercnvpy.tl.infercnv` are stored. With `mode="fast"`, - `adata.obsm[f"X_{use_rep}"]` can also be a dask array or a backed matrix, which is read block by block. - cmap - colormap to use - figsize - (width, height) tuple in inches. Ignored if `ax` is given. + {data} mode `"fast"` computes the group means in one pass over the CNV matrix and draws them with infercnvpy. - `"legacy"` uses :func:`scanpy.pl.heatmap`. `mode="fast"` implements the :func:`scanpy.pl.heatmap` arguments - `swap_axes`, `standard_scale`, `var_group_rotation`, `dendrogram`, `vmin`, `vmax`, `vcenter`, `norm`, - `show_gene_labels` (ignored) and `log`; with any other argument (or a numeric, not categorical, `groupby` - column), it warns and uses `mode="legacy"`, or raises an error if the CNV matrix is not in memory. - ax - Draw the heatmap into this axes (`mode="fast"` only). The group colours, chromosome labels, dendrogram and - colour bar go into axes split off from it. - n_jobs - Number of threads for reading the CNV matrix (`mode="fast"`). `None` uses all cores. - show - Whether to show the figure. Defaults to :attr:`scanpy.settings.autoshow`. - save - If `True` or a `str`, save the figure to :attr:`scanpy.settings.figdir`. A string is appended to the default - filename `heatmap`; it sets the file type if it ends in `.pdf`, `.png` or `.svg`, otherwise - :attr:`scanpy.settings.file_format_figs` does. - return_fig - Return the :class:`~matplotlib.figure.Figure`. + `"legacy"` uses :func:`scanpy.pl.heatmap`. + {fast_args} + {figure} **kwargs - Arguments of :func:`scanpy.pl.heatmap`, see `mode`. By default, the colour scale is a - :class:`~matplotlib.colors.TwoSlopeNorm` centred at `vcenter=0` from the smallest to the largest group mean. - - Returns - ------- - The figure if `return_fig` is true. Otherwise, if `show` is false, a dictionary of axes with the keys - `heatmap_ax`, `groupby_ax` (group colours), `gene_groups_ax` (chromosome labels), `colorbar_ax` (not with - `mode="legacy"`) and, with `dendrogram`, `dendrogram_ax`. + Arguments of :func:`scanpy.pl.heatmap`, see `mode`. By default, the color scale is a + :class:`~matplotlib.colors.TwoSlopeNorm` centered at `vcenter=0` from the smallest to the largest group mean. + + {returns} """ if _use_legacy("chromosome_heatmap_summary", adata, groupby, use_rep, mode, kwargs): - return _legacy(_legacy_summary, adata, groupby, use_rep, cmap, figsize, ax, show, save, return_fig, kwargs) + return _legacy(legacy_summary, adata, groupby, use_rep, cmap, figsize, ax, show, save, return_fig, kwargs) opts = _Options.pop(kwargs) own = ax is None fig, ax = _figure(ax, figsize) @@ -314,6 +298,10 @@ def color_norm(self, lo: float, hi: float) -> Normalize: ) +#: :func:`scanpy.pl.heatmap` arguments that `mode="fast"` implements; any other one switches to `mode="legacy"`. +FAST_KWARGS = frozenset(f.name for f in fields(_Options)) | {"show_gene_labels"} + + # --- mode selection ------------------------------------------------------------------------------------------------- @@ -377,18 +365,6 @@ def _legacy( return _finish(fig, axes, show=show, save=save, return_fig=return_fig) -def _legacy_heatmap(adata: AnnData, **kwargs) -> dict[str, Axes]: - from infercnvpy.pl._chromosome_heatmap_legacy import legacy_heatmap - - return legacy_heatmap(adata, **kwargs) - - -def _legacy_summary(adata: AnnData, **kwargs) -> dict[str, Axes]: - from infercnvpy.pl._chromosome_heatmap_legacy import legacy_summary - - return legacy_summary(adata, **kwargs) - - # --- data ----------------------------------------------------------------------------------------------------------- @@ -415,7 +391,7 @@ def _categories(adata: AnnData, groupby: str) -> pd.Categorical: def _group_colors(adata: AnnData, groupby: str, n: int) -> list: - """The colours scanpy stores for `groupby`, or the default palette scanpy would pick.""" + """The colors scanpy stores for `groupby`, or the default palette scanpy would pick.""" colors = adata.uns.get(f"{groupby}_colors") if colors is not None and len(colors) >= n: return list(colors[:n]) @@ -442,7 +418,7 @@ def _group_means(X: Any, codes: np.ndarray, n_groups: int) -> np.ndarray: def _order_by_dendrogram( adata: AnnData, groupby: str, cat: pd.Categorical, palette: list, compute: Callable[[], dict[str, Any]] ) -> tuple[pd.Categorical, list, dict[str, Any] | None]: - """Reorder the categories (and their colours) as the leaves of the dendrogram of the non-empty groups. + """Reorder the categories (and their colors) as the leaves of the dendrogram of the non-empty groups. The dendrogram is `adata.uns[f"dendrogram_{groupby}"]` from :func:`scanpy.tl.dendrogram` if that has exactly these groups, otherwise `compute()`, which is stored there as `mode="legacy"` does. Unused categories go last. @@ -517,7 +493,7 @@ def _draw( ) -> dict[str, Axes]: """Draw `values` (pixel rows x windows, `rows` rows per category of `cat`) with its annotations. - Cells run down the heatmap, or along it with `swap_axes`. `palette` has a colour per category, `dendro` the + Cells run down the heatmap, or along it with `swap_axes`. `palette` has a color per category, `dendro` the dendrogram of the non-empty groups in their order. """ keep = np.flatnonzero(rows) @@ -679,7 +655,7 @@ def _extent(proto: Text, text: str, renderer: RendererBase, *, along_x: bool) -> def _keep_labels(pos: np.ndarray, extent: np.ndarray, priority: np.ndarray) -> np.ndarray: - """Which labels (centred at `pos`, `extent` long) to draw so that none overlap: greedily, highest priority first. + """Which labels (centered at `pos`, `extent` long) to draw so that none overlap: greedily, highest priority first. Ties go to the earlier label. Labels that do not overlap any other are always drawn. """ diff --git a/src/infercnvpy/pl/_chromosome_heatmap_legacy.py b/src/infercnvpy/pl/_chromosome_heatmap_legacy.py index 9257713..2f2eda7 100644 --- a/src/infercnvpy/pl/_chromosome_heatmap_legacy.py +++ b/src/infercnvpy/pl/_chromosome_heatmap_legacy.py @@ -29,44 +29,23 @@ def legacy_heatmap( """Every cell as a row of :func:`scanpy.pl.heatmap`, in the original order within each group.""" X = _as_matrix(adata.obsm[f"X_{use_rep}"], f"`adata.obsm['X_{use_rep}']`") obs = _as_df(adata.obs) - # an actual AnnData (not a view) that shares `uns`, so the colours scanpy picks are stored in `adata.uns` + # an actual AnnData (not a view) that shares `uns`, so the colors scanpy picks are stored in `adata.uns` order = _stable_group_order(AnnData(X=X, obs=obs, uns=adata.uns), groupby, kwargs.get("dendrogram", False)) tmp_adata = AnnData(X=X[order], obs=obs.iloc[order], uns=adata.uns) - - # re-sort, as saving & loading anndata destroys the order - chr_pos_dict = dict(sorted(adata.uns[use_rep]["chr_pos"].items(), key=lambda x: x[1])) - chr_pos = list(chr_pos_dict.values()) - - # center color map at 0 - tmp_data = X if isinstance(X, np.ndarray) else X.data - vmin = kwargs.pop("vmin", None) - vmax = kwargs.pop("vmax", None) - if vmin is None: - vmin = np.nanmin(tmp_data) - if vmax is None: - vmax = np.nanmax(tmp_data) - kwargs["norm"] = TwoSlopeNorm(0, vmin=vmin, vmax=vmax) - - # add chromosome annotations - var_group_positions = list(zip(chr_pos, chr_pos[1:] + [tmp_adata.shape[1]], strict=False)) - - return_ax_dic = sc.pl.heatmap( + values = X if isinstance(X, np.ndarray) else X.data + return _scanpy_heatmap( tmp_adata, - var_names=tmp_adata.var.index.values, + adata, + use_rep, + np.nanmin(values), + np.nanmax(values), + 0, groupby=groupby, - figsize=figsize, cmap=cmap, - show_gene_labels=False, - var_group_positions=var_group_positions, - var_group_labels=list(chr_pos_dict.keys()), - show=False, - save=False, + figsize=figsize, **kwargs, ) - return_ax_dic["heatmap_ax"].vlines(chr_pos[1:], lw=0.6, ymin=0, ymax=tmp_adata.shape[0]) - return return_ax_dic - def _stable_group_order(adata: AnnData, groupby: str, dendrogram: bool | str) -> np.ndarray: """The cells sorted by group with a stable sort, in the group order `scanpy.pl.heatmap` will use. @@ -118,21 +97,42 @@ def _get_group_mean(group): X=np.vstack([np.repeat(_get_group_mean(group), 10, axis=0) for group in groups]), obs=tmp_obs, uns=adata.uns ) - chr_pos_dict = dict(sorted(adata.uns[use_rep]["chr_pos"].items(), key=lambda x: x[1])) - chr_pos = list(chr_pos_dict.values()) + return _scanpy_heatmap( + tmp_adata, + adata, + use_rep, + np.min(tmp_adata.X), + np.max(tmp_adata.X), + -1, + groupby=groupby, + cmap=cmap, + figsize=figsize, + **kwargs, + ) - # center color map at 0 - vmin = kwargs.pop("vmin", None) - vmax = kwargs.pop("vmax", None) - if vmin is None: - vmin = np.min(tmp_adata.X) - if vmax is None: - vmax = np.max(tmp_adata.X) - kwargs["norm"] = TwoSlopeNorm(0, vmin=vmin, vmax=vmax) - # add chromosome annotations - var_group_positions = list(zip(chr_pos, chr_pos[1:] + [tmp_adata.shape[1]], strict=False)) +def _scanpy_heatmap( + tmp_adata: AnnData, + adata: AnnData, + use_rep: str, + lo: float, + hi: float, + ymin: float, + *, + groupby: str, + cmap: str | Colormap, + figsize: tuple[float, float], + **kwargs, +) -> dict[str, Axes]: + """Draw `tmp_adata` with :func:`scanpy.pl.heatmap`, with chromosome annotations. + The color scale is centered at 0, from `vmin` to `vmax` (default `lo` and `hi`). + """ + # re-sort, as saving & loading anndata destroys the order + chr_pos_dict = dict(sorted(adata.uns[use_rep]["chr_pos"].items(), key=lambda x: x[1])) + chr_pos = list(chr_pos_dict.values()) + vmin, vmax = kwargs.pop("vmin", None), kwargs.pop("vmax", None) + kwargs["norm"] = TwoSlopeNorm(0, vmin=lo if vmin is None else vmin, vmax=hi if vmax is None else vmax) return_ax_dic = sc.pl.heatmap( tmp_adata, var_names=tmp_adata.var.index.values, @@ -140,12 +140,11 @@ def _get_group_mean(group): figsize=figsize, cmap=cmap, show_gene_labels=False, - var_group_positions=var_group_positions, + var_group_positions=list(zip(chr_pos, chr_pos[1:] + [tmp_adata.shape[1]], strict=False)), var_group_labels=list(chr_pos_dict.keys()), show=False, save=False, **kwargs, ) - - return_ax_dic["heatmap_ax"].vlines(chr_pos[1:], lw=0.6, ymin=-1, ymax=tmp_adata.shape[0]) + return_ax_dic["heatmap_ax"].vlines(chr_pos[1:], lw=0.6, ymin=ymin, ymax=tmp_adata.shape[0]) return return_ax_dic diff --git a/src/infercnvpy/pl/_heatmap_kernels.py b/src/infercnvpy/pl/_heatmap_kernels.py index fc0c3ce..8518d56 100644 --- a/src/infercnvpy/pl/_heatmap_kernels.py +++ b/src/infercnvpy/pl/_heatmap_kernels.py @@ -2,12 +2,12 @@ The cells are sorted by group and every group is split into contiguous runs of cells, one per pixel row (`pixel_rows`). A pixel row then shows either one cell of its run (`sampled_rows`: the cell in the middle, which is -what matplotlib's nearest-neighbour downscaling shows when every cell is drawn) or the mean of the run (`binned_sums`: +what matplotlib's nearest-neighbor downscaling shows when every cell is drawn) or the mean of the run (`binned_sums`: `_bin_sums_csr` / `_bin_sums_dense` add up the rows of each run straight from the matrix, CSR without densifying). A figure is ~1000 pixels tall, so the result is a small dense matrix whatever the number of cells. -The matrix is read in row blocks (`infercnvpy.tl._infercnv._Rows`): an in-memory matrix as one block with the parallel -kernel, a dask array block by block and a backed matrix `CHUNKSIZE` rows at a time, with the serial kernel. +The matrix is read in row blocks (`infercnvpy.tl._infercnv._Rows`): an in-memory matrix as one block, a dask array +block by block and a backed matrix `CHUNKSIZE` rows at a time. """ from __future__ import annotations @@ -78,14 +78,6 @@ def _blocks(X: Any) -> _Rows: return _Rows(X, max(X.shape[0], 1) if isinstance(X, np.ndarray) or sp.issparse(X) else CHUNKSIZE) -def _each_block(source: _Rows, fn) -> None: - """Call `fn(X, row0, n_rows, start)` for every block of `source`, in order (backed and dask blocks read ahead).""" - blocks = list(range(source.n_blocks)) - for i, (Xb, row0) in zip(blocks, source.read_ahead(blocks), strict=True): - s, e = int(source.bounds[i]), int(source.bounds[i + 1]) - fn(Xb, row0, e - s, s) - - def binned_sums(X: Any, order: np.ndarray, bounds: np.ndarray) -> np.ndarray: """Sum of the rows of `X` in each run `order[bounds[b]:bounds[b + 1]]` (float64, runs x columns). @@ -95,53 +87,42 @@ def binned_sums(X: Any, order: np.ndarray, bounds: np.ndarray) -> np.ndarray: sums = np.zeros((n_bins, X.shape[1]), dtype=np.float64) row_bin = np.full(X.shape[0], -1, dtype=np.int64) row_bin[order] = np.repeat(np.arange(n_bins), np.diff(bounds)) - source = _blocks(X) - serial = source.kind != "in memory" - - def add(Xb, row0, n, start): - local = np.flatnonzero(row_bin[start : start + n] >= 0) - b = row_bin[start + local] + for Xb, row0, s, e in _blocks(X).read_ahead(): + local = np.flatnonzero(row_bin[s:e] >= 0) + b = row_bin[s + local] by_bin = np.argsort(b, kind="stable") bins, first = np.unique(b[by_bin], return_index=True) o = (row0 + local[by_bin]).astype(np.int64) bnd = np.append(first, len(o)).astype(np.int64) if sp.issparse(Xb): - _bin_sums_csr(Xb.indptr, Xb.indices, Xb.data, o, bnd, bins, sums, serial=serial) + _bin_sums_csr(Xb.indptr, Xb.indices, Xb.data, o, bnd, bins, sums) else: - _bin_sums_dense(Xb, o, bnd, bins, sums, serial=serial) - - _each_block(source, add) + _bin_sums_dense(Xb, o, bnd, bins, sums) return sums def column_ranges(X: Any) -> tuple[np.ndarray, np.ndarray]: """Minimum and maximum of each column of `X` over all rows (with the zeros of a sparse matrix), NaNs ignored.""" lo, hi = np.full(X.shape[1], np.inf), np.full(X.shape[1], -np.inf) - - def update(Xb, row0, n, start): - rows = Xb[row0 : row0 + n] + for Xb, row0, s, e in _blocks(X).read_ahead(): + rows = Xb[row0 : row0 + e - s] if sp.issparse(rows): b_lo, b_hi = rows.nanmin(axis=0).toarray().ravel(), rows.nanmax(axis=0).toarray().ravel() else: b_lo, b_hi = np.fmin.reduce(rows, axis=0), np.fmax.reduce(rows, axis=0) np.fmin(lo, b_lo, out=lo) np.fmax(hi, b_hi, out=hi) - - _each_block(_blocks(X), update) return lo, hi def gram(X: Any) -> tuple[np.ndarray, np.ndarray]: """`X.T @ X` and the column sums of `X`, in float64, `CHUNKSIZE` rows at a time.""" G, sums = np.zeros((X.shape[1], X.shape[1])), np.zeros(X.shape[1]) - - def add(Xb, row0, n, start): - for s in range(row0, row0 + n, CHUNKSIZE): - rows = Xb[s : min(s + CHUNKSIZE, row0 + n)].astype(np.float64) - G[:] += (rows.T @ rows).toarray() if sp.issparse(rows) else rows.T @ rows - sums[:] += np.asarray(rows.sum(axis=0)).ravel() - - _each_block(_blocks(X), add) + for Xb, row0, s, e in _blocks(X).read_ahead(): + for r in range(row0, row0 + e - s, CHUNKSIZE): + rows = Xb[r : min(r + CHUNKSIZE, row0 + e - s)].astype(np.float64) + G += (rows.T @ rows).toarray() if sp.issparse(rows) else rows.T @ rows + sums += np.asarray(rows.sum(axis=0)).ravel() return G, sums @@ -153,18 +134,14 @@ def sampled_rows(X: Any, picks: np.ndarray) -> tuple[np.ndarray, float, float]: out = np.empty((len(picks), X.shape[1]), dtype=np.float64) sorted_picks = np.argsort(picks, kind="stable") ranges = [] - - def take(Xb, row0, n, start): - values = Xb.data if sp.issparse(Xb) else Xb[row0 : row0 + n] + for Xb, row0, s, e in _blocks(X).read_ahead(): + values = Xb.data if sp.issparse(Xb) else Xb[row0 : row0 + e - s] if values.size: ranges.append((np.nanmin(values), np.nanmax(values))) - lo, hi = np.searchsorted(picks[sorted_picks], [start, start + n]) + lo, hi = np.searchsorted(picks[sorted_picks], [s, e]) dest = sorted_picks[lo:hi] - rows = Xb[row0 + picks[dest] - start] + rows = Xb[row0 + picks[dest] - s] out[dest] = rows.toarray() if sp.issparse(rows) else rows - - source = _blocks(X) - _each_block(source, take) if not ranges: return out, np.nan, np.nan return out, float(min(r[0] for r in ranges)), float(max(r[1] for r in ranges)) diff --git a/src/infercnvpy/pl/_spatial.py b/src/infercnvpy/pl/_spatial.py index 4aeb66a..426b560 100644 --- a/src/infercnvpy/pl/_spatial.py +++ b/src/infercnvpy/pl/_spatial.py @@ -50,8 +50,8 @@ def spatial_cnv( """Plot the mean CNV of genomic regions, or columns of `adata.obs`, on the tissue. Thin wrapper around :func:`scanpy.pl.embedding` with `basis="spatial"`. Each region - becomes a panel coloured by the mean of `adata.obsm["X_{use_rep}"]` over its windows, - with a diverging colour map centred at 0 (pass `vcenter=None` to turn that off). + becomes a panel colored by the mean of `adata.obsm["X_{use_rep}"]` over its windows, + with a diverging color map centered at 0 (pass `vcenter=None` to turn that off). As :func:`scanpy.pl.spatial` without an image, the panels have an equal aspect and the y axis points down. `adata.obsm["X_{use_rep}"]` must be in memory. @@ -72,7 +72,7 @@ def spatial_cnv( basis Key of the coordinates in `adata.obsm` (`"spatial"` refers to `obsm["spatial"]`). cmap - Colour map of the regions (and other continuous panels). + Color map of the regions (and other continuous panels). **kwargs Passed to :func:`scanpy.pl.embedding`, e.g. `size`, `ncols`, `vmin`/`vmax`, `show`, `save` (a string is appended to the file name `basis`) or `return_fig`. diff --git a/src/infercnvpy/tl/_infercnv.py b/src/infercnvpy/tl/_infercnv.py index 79b4ec3..e17e34c 100644 --- a/src/infercnvpy/tl/_infercnv.py +++ b/src/infercnvpy/tl/_infercnv.py @@ -14,7 +14,7 @@ from numpy.lib.stride_tricks import sliding_window_view from infercnvpy._logging import SKIP_OWN_FRAMES, log, log_result -from infercnvpy._numba import n_kernel_chunks, numba_threads +from infercnvpy._numba import dask_task, n_kernel_chunks, numba_threads from infercnvpy._util import _as_matrix, _is_backed, _is_dask, _obs_column, _params_for_uns from infercnvpy.tl._infercnv_kernels import ( _cnv_csr, @@ -22,8 +22,6 @@ _count_kept, _fill_kept, _gene_values, - _group_expm1_sums_csr, - _group_expm1_sums_dense, _group_sums_csr, _group_sums_dense, _n_unsorted_rows, @@ -144,14 +142,14 @@ def infercnv( are NaN. This needs a dense float32 cells x genes matrix in memory. min_mean Leave out genes whose mean expression over all cells is below this value, before the - windows are formed. The mean is taken of `expm1` of the input, i.e. of the normalised - counts behind the log-normalised values (see :ref:`input-data`), not of the log values. + windows are formed. The mean is taken of `expm1` of the input, i.e. of the normalized + counts behind the log-normalized values (see :ref:`input-data`), not of the log values. This follows the `cutoff` of R inferCNV, which removes genes by their mean read count and suggests `1` for Smart-seq2 and `0.1` for 10x data. R applies its cutoff to the raw - counts before normalisation; the means here are on the scale of + counts before normalization; the means here are on the scale of :func:`scanpy.pp.normalize_total`, which are on the scale of the raw counts when the data - is normalised to the median total count (scanpy's default `target_sum=None`, as R - inferCNV normalises). Default: keep all genes. + is normalized to the median total count (scanpy's default `target_sum=None`, as R + inferCNV normalizes). Default: keep all genes. Returns ------- @@ -369,7 +367,6 @@ def __init__(self, var: pd.DataFrame, var_mask: np.ndarray, window_size: int, st self.n_windows = n_windows g0s, gs, As, Bs, Ks, w0s = (np.array(x, dtype=np.int64) for x in zip(*geometry, strict=True)) self.geometry = (g0s, gs, As, Bs, Ks, w0s, step) - self.step = step first, n_in, offset, length, ab = (np.concatenate(x).astype(np.int64) for x in zip(*scatter, strict=True)) #: largest number of windows that contain one gene self.span = int(n_in.max()) @@ -419,9 +416,9 @@ def build(cls, layout: _Layout, reference: np.ndarray, lfc_clip: float, dynamic_ def _windows( - plan: _Plan, X, rows: np.ndarray, band=None, *, out: bool = True, serial: bool = False + plan: _Plan, X, rows: np.ndarray, band=None, *, out: bool = True ) -> tuple[np.ndarray, np.ndarray, np.ndarray]: - """Median-centred windows (float32) of the rows `rows` of `X`, their (sum, sum of squares), and their counts. + """Median-centered windows (float32) of the rows `rows` of `X`, their (sum, sum of squares), and their counts. The counts are the numbers of values outside the noise `band`, all 0 without a band. With `out=False`, the windows are not kept (the first array has no rows). @@ -432,15 +429,15 @@ def _windows( rowstats = np.empty((n_rows, 2), dtype=np.float64) counts = np.empty(n_rows, dtype=np.int64) thr, centre = (np.inf, 0.0) if band is None else band - n_chunks = n_kernel_chunks(n_rows, serial=serial) + n_chunks = n_kernel_chunks(n_rows) shared = ( rows, lay.col_map, plan.lo, plan.hi, plan.base, plan.cap, *lay.geometry, *lay.scatter, lay.window_weight, plan.base_sums, plan.max_nnz, thr, centre, windows, rowstats, counts, n_chunks, ) # fmt: skip if sp.issparse(X): - _cnv_csr(X.data, X.indices, X.indptr, *shared, serial=serial) + _cnv_csr(X.data, X.indices, X.indptr, *shared) else: - _cnv_dense(X, *shared, serial=serial) + _cnv_dense(X, *shared) if np.isnan(rowstats[:, 0]).any(): raise ValueError("infercnv needs finite values, but the expression matrix or the reference contains NaN.") return windows, rowstats, counts @@ -491,7 +488,7 @@ def __init__(self, n_rows: int, n_cols: int): self.indices = np.empty(0, dtype=idx_dtype) self.n = 0 - def add(self, windows: np.ndarray, counts: np.ndarray, band: tuple[float, float], serial: bool = False): + def add(self, windows: np.ndarray, counts: np.ndarray, band: tuple[float, float]): """Append the rows `windows`, with `counts` values outside the noise `band` per row.""" r0, r1 = self.n, self.n + len(windows) indptr = self.indptr[r0 : r1 + 1] @@ -502,7 +499,7 @@ def add(self, windows: np.ndarray, counts: np.ndarray, band: tuple[float, float] size = max(nnz, int(nnz / r1 * self.shape[0] * 1.05), self.data.size * 5 // 4) self.data.resize(size, refcheck=False) self.indices.resize(size, refcheck=False) - _fill_kept(windows, *band, indptr, self.data, self.indices, serial=serial) + _fill_kept(windows, *band, indptr, self.data, self.indices) self.n = r1 def finish(self) -> sp.csr_matrix: @@ -514,17 +511,19 @@ def finish(self) -> sp.csr_matrix: return sp.csr_matrix((self.data, indices, self.indptr.astype(idx_dtype)), shape=self.shape, copy=False) -def _to_csr(windows: np.ndarray, counts: np.ndarray, band: tuple[float, float], serial: bool = False) -> sp.csr_matrix: +def _to_csr(windows: np.ndarray, counts: np.ndarray, band: tuple[float, float]) -> sp.csr_matrix: """The entries of `windows` outside the noise band (`counts` per row) as a CSR matrix.""" builder = _CsrBuilder(*windows.shape) - builder.add(windows, counts, band, serial=serial) + builder.add(windows, counts, band) return builder.finish() -def _gene_values_of(plan: _Plan, out: np.ndarray, band: tuple[float, float], serial: bool = False) -> np.ndarray: - gv = np.full((out.shape[0], plan.layout.n_vars), np.nan, dtype=np.float32) - n_chunks = n_kernel_chunks(out.shape[0], serial=serial) - _gene_values(out, *plan.layout.gene_geometry, *band, gv, n_chunks, serial=serial) +def _gene_values_of(plan: _Plan, out: np.ndarray, band: tuple[float, float], gv: np.ndarray | None = None): + """Gene values of the centered windows `out` into `gv` (NaN-filled; default: a new array), see `_gene_values`.""" + if gv is None: + gv = np.full((out.shape[0], plan.layout.n_vars), np.nan, dtype=np.float32) + n_chunks = n_kernel_chunks(out.shape[0]) + _gene_values(out, *plan.layout.gene_geometry, *band, gv, n_chunks) return gv @@ -552,14 +551,14 @@ def _run(plan: _Plan, source: "_Rows", band_rows: np.ndarray | None, gene_values band = _noise_band(plan, rowstats if band_rows is None else rowstats[band_rows]) _count_kept(windows, *band, counts) if gv is not None: - gv[builder.n : builder.n + len(rows)] = _gene_values_of(plan, windows, band) + _gene_values_of(plan, windows, band, gv[builder.n : builder.n + len(rows)]) builder.add(windows, counts, band) assert band is not None return builder.finish(), gv, band def _run_lazy(plan: _Plan, source: "_Rows", band_rows: np.ndarray | None, gene_values: bool): - """The CNV matrix (and gene values) as lazy dask arrays; the serial kernels run inside dask tasks. + """The CNV matrix (and gene values) as lazy dask arrays; the kernels run serially inside the dask tasks. The noise band is computed first, by one dask computation over the blocks that hold `band_rows` (or all blocks). """ @@ -589,20 +588,23 @@ def _run_lazy(plan: _Plan, source: "_Rows", band_rows: np.ndarray | None, gene_v return res, gv, band +@dask_task def _dask_stats_task(block, plan: _Plan, local: np.ndarray | None) -> np.ndarray: X = _kernel_ready(block) - return _windows(plan, X, np.arange(X.shape[0]) if local is None else local, out=False, serial=True)[1] + return _windows(plan, X, np.arange(X.shape[0]) if local is None else local, out=False)[1] +@dask_task def _lazy_cnv_block(plan: _Plan, band, block) -> sp.csr_matrix: X = _kernel_ready(block) - windows, _, counts = _windows(plan, X, np.arange(X.shape[0]), band, serial=True) - return _to_csr(windows, counts, band, serial=True) + windows, _, counts = _windows(plan, X, np.arange(X.shape[0]), band) + return _to_csr(windows, counts, band) +@dask_task def _lazy_gene_values_block(plan: _Plan, band, block) -> np.ndarray: X = _kernel_ready(block) - return _gene_values_of(plan, _windows(plan, X, np.arange(X.shape[0]), serial=True)[0], band, serial=True) + return _gene_values_of(plan, _windows(plan, X, np.arange(X.shape[0]))[0], band) class _Rows: @@ -610,7 +612,7 @@ class _Rows: In memory and backed, blocks have `chunksize` rows; a dask array keeps its own blocks. `rows_of` reads the blocks ahead in threads (dask: with the synchronous scheduler) and yields them on this thread, so the kernels can use all - numba threads. + numba threads. A backed CSC matrix is read into memory as CSR first: each row block of it would read all of it. """ def __init__(self, X, chunksize: int): @@ -621,6 +623,15 @@ def __init__(self, X, chunksize: int): X = X.rechunk({1: -1}) sizes = list(X.chunks[0]) else: + if _is_backed(X) and getattr(X, "format", None) == "csc": + warnings.warn( + "The backed matrix is stored as CSC, so reading a block of its rows reads all of it. infercnvpy " + "reads it into memory once, as CSR. To read it block by block, store it as CSR, e.g. with " + "`adata.X = adata.X.tocsr()` before `adata.write_h5ad`.", + UserWarning, + skip_file_prefixes=SKIP_OWN_FRAMES, + ) + X = X.to_memory() self.kind = "backed" if _is_backed(X) else "in memory" if self.kind == "in memory": X = _kernel_ready(X) @@ -649,12 +660,15 @@ def rows_of(self, rows: np.ndarray | None) -> Iterator[tuple[Any, np.ndarray]]: `X` is an in-memory, kernel-ready matrix that holds the block. """ blocks = list(self.blocks_of(rows)) - for (i, local), (X, row0) in zip(blocks, self.read_ahead([i for i, _ in blocks]), strict=True): - s, e = int(self.bounds[i]), int(self.bounds[i + 1]) + for (_, local), (X, row0, s, e) in zip(blocks, self.read_ahead([i for i, _ in blocks]), strict=True): yield X, row0 + (np.arange(e - s) if local is None else local) - def read_ahead(self, blocks: list[int]) -> Iterator[tuple[Any, int]]: - """`read` of each of `blocks`. Backed and dask blocks are read ahead in threads while earlier ones are used.""" + def read_ahead(self, blocks: Sequence[int] | None = None) -> Iterator[tuple[Any, int, int, int]]: + """`read` of each of `blocks` (default: all), in order. + + Backed and dask blocks are read ahead in threads while earlier ones are used. + """ + blocks = range(self.n_blocks) if blocks is None else blocks if self.kind == "in memory" or len(blocks) < 2: yield from map(self.read, blocks) return @@ -668,13 +682,13 @@ def read_ahead(self, blocks: list[int]) -> Iterator[tuple[Any, int]]: while ahead: yield ahead.popleft().result() - def read(self, i: int) -> tuple[Any, int]: - """Block `i` as an in-memory, kernel-ready matrix, and the row of that matrix where the block starts.""" + def read(self, i: int) -> tuple[Any, int, int, int]: + """`(X, row0, s, e)`: rows `s:e`, block `i`, are the rows from `row0` on of `X` (in memory, kernel-ready).""" s, e = int(self.bounds[i]), int(self.bounds[i + 1]) if self.kind == "in memory": - return self.X, s + return self.X, s, s, e block = self.X.blocks[i].compute(scheduler="synchronous") if self.kind == "dask" else self.X[s:e] - return _kernel_ready(block), 0 + return _kernel_ready(block), 0, s, e def _kernel_ready(X): @@ -781,7 +795,7 @@ def _means(source: _Rows, codes: np.ndarray, n_groups: int, *, lazy: bool, expm1 sums = np.sum(dask.compute(*tasks), axis=0) else: sums = np.zeros((n_groups, source.X.shape[1])) - for (i, local), (X, row0) in zip(blocks, source.read_ahead([i for i, _ in blocks]), strict=True): + for (i, local), (X, row0, _, _) in zip(blocks, source.read_ahead([i for i, _ in blocks]), strict=True): local, groups = _block_groups(source, i, local, codes) sums += _group_sums(X, row0 + local, groups, n_groups, expm1=expm1) counts = np.bincount(codes[codes >= 0], minlength=n_groups) @@ -796,6 +810,7 @@ def _block_groups(source: _Rows, i: int, local: np.ndarray | None, codes: np.nda return local, codes[s + local] +@dask_task def _dask_block_sums(block, local: np.ndarray, groups: np.ndarray, n_groups: int, expm1: bool = False) -> np.ndarray: return _group_sums(_kernel_ready(block), local, groups, n_groups, expm1=expm1) @@ -807,7 +822,7 @@ def _group_sums(X, rows: np.ndarray, groups: np.ndarray, n_groups: int, *, expm1 """ sums = np.zeros((n_kernel_chunks(len(rows), per_thread=1), n_groups, X.shape[1])) if sp.issparse(X): - (_group_expm1_sums_csr if expm1 else _group_sums_csr)(X.data, X.indices, X.indptr, rows, groups, sums) + _group_sums_csr(X.data, X.indices, X.indptr, rows, groups, expm1, sums) else: - (_group_expm1_sums_dense if expm1 else _group_sums_dense)(X, rows, groups, sums) + _group_sums_dense(X, rows, groups, expm1, sums) return sums.sum(axis=0) diff --git a/src/infercnvpy/tl/_infercnv_kernels.py b/src/infercnvpy/tl/_infercnv_kernels.py index c79466e..626510f 100644 --- a/src/infercnvpy/tl/_infercnv_kernels.py +++ b/src/infercnvpy/tl/_infercnv_kernels.py @@ -15,7 +15,7 @@ squares of the written values for the noise threshold and the number of values outside the noise band. Scratch buffers are allocated once per call, one row per chunk of cells, outside the `prange` loop: allocating inside -the loop body serialises the threads on numba's allocator lock. `np.median` would allocate too, hence `_median`. +the loop body serializes the threads on numba's allocator lock. `np.median` would allocate too, hence `_median`. """ from __future__ import annotations @@ -135,6 +135,25 @@ def _normalise(win, n_windows, window_weight): win[w] /= window_weight[w] +@njit +def _kept(v, thr, centre): + """Whether a centered window value is kept: non-zero and outside the noise band `centre ± thr`.""" + return v != 0 and abs(np.float64(v) - centre) >= thr + + +@njit +def _scratch(n_chunks, n_genes, gs, n_windows, span): + """Scratch rows per chunk of cells: gene values, the two prefix sums, window sums (`span` spare), median buffer.""" + max_g = gs.max() + 1 + return ( + np.empty((n_chunks, n_genes), np.float32), + np.empty((n_chunks, max_g)), + np.empty((n_chunks, max_g)), + np.empty((n_chunks, n_windows + span)), + np.empty((n_chunks, n_windows)), + ) + + @njit def _center_emit(win, buf, thr, centre, out, r, rowstats, counts): """Subtract the median from the running means `win`, write them to `out[r]` as float32. @@ -153,7 +172,7 @@ def _center_emit(win, buf, thr, centre, out, r, rowstats, counts): out[r, w] = v total += v squares += np.float64(v) * v - if v != 0 and abs(np.float64(v) - centre) >= thr: + if _kept(v, thr, centre): kept += 1 rowstats[r, 0] = total rowstats[r, 1] = squares @@ -165,7 +184,7 @@ def _cnv_csr( data, indices, indptr, rows, col_map, lo, hi, base, cap, g0s, gs, As, Bs, Ks, w0s, step, first, weights, window_weight, base_sums, max_nnz, thr, centre, out, rowstats, counts, n_chunks, ): # fmt: skip - """Centred windows of the CSR rows `rows` into `out`, see `_center_emit`. + """Centered windows of the CSR rows `rows` into `out`, see `_center_emit`. A row with at most `max_nnz` stored values starts from `base_sums`, the weighted window sums of the baseline `base`, and adds the difference to the baseline of each stored value to the windows that contain it (`_scatter`). Other @@ -174,12 +193,7 @@ def _cnv_csr( """ n_rows = rows.size n_windows = out.shape[1] - max_g = gs.max() + 1 - x = np.empty((n_chunks, base.size), np.float32) - p1 = np.empty((n_chunks, max_g)) - p2 = np.empty((n_chunks, max_g)) - win = np.empty((n_chunks, n_windows + weights.shape[1])) - buf = np.empty((n_chunks, n_windows)) + x, p1, p2, win, buf = _scratch(n_chunks, base.size, gs, n_windows, weights.shape[1]) for t in prange(n_chunks): xt = x[t] wt = win[t] @@ -213,12 +227,7 @@ def _cnv_dense( """Same as `_cnv_csr` for the rows `rows` of the dense matrix `X`; values that differ from the baseline count.""" n_rows = rows.size n_windows = out.shape[1] - max_g = gs.max() + 1 - x = np.empty((n_chunks, base.size), np.float32) - p1 = np.empty((n_chunks, max_g)) - p2 = np.empty((n_chunks, max_g)) - win = np.empty((n_chunks, n_windows + weights.shape[1])) - buf = np.empty((n_chunks, n_windows)) + x, p1, p2, win, buf = _scratch(n_chunks, base.size, gs, n_windows, weights.shape[1]) pos = np.empty((n_chunks, max_nnz + 1), np.int64) diff = np.empty((n_chunks, max_nnz + 1)) n_cols = X.shape[1] @@ -256,8 +265,7 @@ def _count_kept(out, thr, centre, counts): for r in prange(n_rows): c = 0 for w in range(n_windows): - v = out[r, w] - if v != 0 and abs(np.float64(v) - centre) >= thr: + if _kept(out[r, w], thr, centre): c += 1 counts[r] = c @@ -270,7 +278,7 @@ def _fill_kept(out, thr, centre, indptr, data, indices): k = indptr[r] for w in range(n_windows): v = out[r, w] - if v != 0 and abs(np.float64(v) - centre) >= thr: + if _kept(v, thr, centre): data[k] = v indices[k] = w k += 1 @@ -278,9 +286,9 @@ def _fill_kept(out, thr, centre, indptr, data, indices): @njit_parallel def _gene_values(out, g0s, gs, Ks, w0s, first, n_in, cols, thr, centre, gv, n_chunks): - """Per gene: mean of the kept windows that contain it, median-centred over these genes, denoised. + """Per gene: mean of the kept windows that contain it, median-centered over these genes, denoised. - `out` holds the centred windows of a block of cells. Per gene in genomic order, `first` is the first window that + `out` holds the centered windows of a block of cells. Per gene in genomic order, `first` is the first window that contains it and `n_in` how many do (genes in none are skipped). Values are written into `gv` (cells x all genes of `adata`, prefilled with NaN) at the column `cols[gene]`. A window prefix sum makes this O(1) per gene. """ @@ -314,36 +322,11 @@ def _gene_values(out, g0s, gs, Ks, w0s, first, n_in, cols, thr, centre, gv, n_ch @njit_parallel -def _group_sums_csr(data, indices, indptr, rows, groups, sums): - """`sums[t, groups[i]] += X[rows[i]]` for the CSR matrix `X`, where `t` is the chunk of `i` (one per `sums[t]`).""" - n_chunks = sums.shape[0] - n = rows.size - for t in prange(n_chunks): - st = sums[t] - for i in range(t * n // n_chunks, (t + 1) * n // n_chunks): - s = st[groups[i]] - r = rows[i] - for k in range(indptr[r], indptr[r + 1]): - s[indices[k]] += data[k] - - -@njit_parallel -def _group_sums_dense(X, rows, groups, sums): - """Same as `_group_sums_csr` for the dense matrix `X`.""" - n_chunks = sums.shape[0] - n = rows.size - for t in prange(n_chunks): - st = sums[t] - for i in range(t * n // n_chunks, (t + 1) * n // n_chunks): - s = st[groups[i]] - r = rows[i] - for j in range(X.shape[1]): - s[j] += X[r, j] +def _group_sums_csr(data, indices, indptr, rows, groups, expm1, sums): + """`sums[t, groups[i]] += X[rows[i]]` for the CSR matrix `X`, where `t` is the chunk of `i` (one per `sums[t]`). - -@njit_parallel -def _group_expm1_sums_csr(data, indices, indptr, rows, groups, sums): - """Like `_group_sums_csr`, of `expm1(X)` (the counts behind log-normalised values; implicit zeros stay 0).""" + With `expm1`, the sums are of `expm1(X)` (the counts behind log-normalized values; implicit zeros stay 0). + """ n_chunks = sums.shape[0] n = rows.size for t in prange(n_chunks): @@ -352,12 +335,13 @@ def _group_expm1_sums_csr(data, indices, indptr, rows, groups, sums): s = st[groups[i]] r = rows[i] for k in range(indptr[r], indptr[r + 1]): - s[indices[k]] += np.expm1(np.float64(data[k])) + v = np.float64(data[k]) + s[indices[k]] += np.expm1(v) if expm1 else v @njit_parallel -def _group_expm1_sums_dense(X, rows, groups, sums): - """Same as `_group_expm1_sums_csr` for the dense matrix `X`.""" +def _group_sums_dense(X, rows, groups, expm1, sums): + """Same as `_group_sums_csr` for the dense matrix `X`.""" n_chunks = sums.shape[0] n = rows.size for t in prange(n_chunks): @@ -366,7 +350,8 @@ def _group_expm1_sums_dense(X, rows, groups, sums): s = st[groups[i]] r = rows[i] for j in range(X.shape[1]): - s[j] += np.expm1(np.float64(X[r, j])) + v = np.float64(X[r, j]) + s[j] += np.expm1(v) if expm1 else v @njit_parallel diff --git a/src/infercnvpy/tl/_ith.py b/src/infercnvpy/tl/_ith.py index cb1b074..34504cd 100644 --- a/src/infercnvpy/tl/_ith.py +++ b/src/infercnvpy/tl/_ith.py @@ -1,7 +1,7 @@ """Interquartile range of the cell-cell Pearson correlations within a group, without the cells x cells matrix. `np.percentile(np.corrcoef(X), [75, 25])` takes quantiles over all n^2 entries of the correlation matrix: the n -diagonal ones and every pair twice. Here the rows are standardized once (float32, centred, unit norm), the +diagonal ones and every pair twice. Here the rows are standardized once (float32, centered, unit norm), the correlations of each block of rows with all later rows come from one BLAS product, and a parallel numba kernel adds that block's upper triangle to thread-local histograms with weight 2. The quantiles are read off the cumulative counts (linear interpolation inside a bin, then numpy's linear interpolation between ranks), which @@ -26,7 +26,7 @@ import numpy as np import scipy.sparse as sp -from infercnvpy._numba import _on_main_thread, njit_parallel, prange +from infercnvpy._numba import njit_parallel, prange, runs_serially N_BINS = 1 << 20 # bin width 2 / 2^20 ~ 2e-6 TILE = 2048 # rows per block: a block of correlations is TILE x cells float32 @@ -34,7 +34,7 @@ @njit_parallel def _standardize_dense(X, out): # pragma: no cover - """Centre each row and scale it to unit norm; a constant row becomes NaN. `out` may be `X`.""" + """Center each row and scale it to unit norm; a constant row becomes NaN. `out` may be `X`.""" n, g = X.shape for i in prange(n): m = 0.0 @@ -84,7 +84,7 @@ def _sampled_corr(Z, rows, cols, out): # pragma: no cover def _quantiles(hist: np.ndarray, qs: tuple[float, ...]) -> list[float]: - """Numpy's default (linear) quantiles of the values summarised by `hist` over [-1, 1].""" + """Numpy's default (linear) quantiles of the values summarized by `hist` over [-1, 1].""" n_bins = len(hist) cum = np.cumsum(hist) total = int(cum[-1]) @@ -152,7 +152,7 @@ def corr_iqr( C = C[~np.eye(n, dtype=bool)] q75, q25 = np.percentile(C, [75, 25]) return q75 - q25, n_dropped, False - hists = np.zeros((numba.get_num_threads() if _on_main_thread() else 1, N_BINS), dtype=np.int64) + hists = np.zeros((1 if runs_serially() else numba.get_num_threads(), N_BINS), dtype=np.int64) for r0 in range(0, n, tile): _bin_tile(Z[r0 : r0 + tile] @ Z[r0:].T, r0, hists) hist = hists.sum(axis=0) diff --git a/src/infercnvpy/tl/_scores.py b/src/infercnvpy/tl/_scores.py index dba10b2..360c642 100644 --- a/src/infercnvpy/tl/_scores.py +++ b/src/infercnvpy/tl/_scores.py @@ -18,6 +18,7 @@ _is_dask, _obs_column, _params_for_uns, + _splice_doc, ) from ._ith import corr_iqr @@ -111,10 +112,7 @@ def _abs_group_means(X, groups: pd.Series) -> dict[Any, np.number]: codes = pd.Index(uniq).get_indexer(groups) codes[pd.isna(groups.to_numpy())] = n # cells without a group (a NaN group gets no cells, as in memory) sums = np.zeros(n + 1) - source = _Rows(X, 10_000) - blocks = list(range(source.n_blocks)) - for i, (Xb, row0) in zip(blocks, source.read_ahead(blocks), strict=True): - s, e = int(source.bounds[i]), int(source.bounds[i + 1]) + for Xb, row0, s, e in _Rows(X, 10_000).read_ahead(): row_sums = np.asarray(abs(Xb[row0 : row0 + e - s]).sum(axis=1), dtype=np.float64).ravel() sums += np.bincount(codes[s:e], weights=row_sums, minlength=n + 1) counts = np.bincount(codes, minlength=n + 1)[:n] * X.shape[1] @@ -123,7 +121,66 @@ def _abs_group_means(X, groups: pd.Series) -> dict[Any, np.number]: return dict(zip(uniq, means, strict=True)) +#: Docstring blocks shared by `ithgex` and `ithcna`. +_ITH_DOC = { + "wu2021": """ + Intratumoral heterogeneity scores based on CNAs and gene expressions + The calculations of intratumoral heterogeneity scores were inspired by a + previous study and modified as follows. First, to calculate ITHCNA, we used + the relative expression value matrix generated by inferCNV and calculated the + pairwise cell–cell distances using Pearson's correlation coefficients for each + patient. ITHCNA was defined as interquartile range (IQR) of the distribution for + all malignant cell pairs' Pearson's correlation coefficients. + """, + "pairs": """ + The correlations are taken over all ordered pairs of cells in a group, including each + cell with itself (self-pairs, see `exclude_self_pairs`), as in earlier versions. + """, + "changed": """ + .. versionchanged:: 0.7.0 + The correlations are computed block by block with a numba kernel instead of as one + cells x cells matrix, so large groups no longer run out of memory. The quantiles match + numpy to about 1e-6. Cells with a constant profile (e.g. no CNV signal at all) have no + defined correlation and are left out; before, any such cell made the group's score `NaN`. + The number of cells left out is logged per group. + `groupby` is keyword-only; passing it positionally still works with a :class:`FutureWarning` + until 0.8.0. + """, + "added": """ + .. versionadded:: 0.7.0 + `n_pairs` and `exclude_self_pairs`. + """, + "params": """ + n_jobs + Number of threads. Default: all cores (numba's setting, see :func:`numba.set_num_threads`). + n_pairs + Estimate the score of large groups from this many randomly drawn cell pairs instead of + all pairs, e.g. `1_000_000`. Pairs are drawn uniformly from the full correlation matrix, + so the estimate is consistent (it converges to the exact score as `n_pairs` grows); its + error shrinks like `1 / sqrt(n_pairs)` and the cost no longer grows with the group size. + Groups with fewer pairs are computed exactly. Default: all pairs (exact). + random_state + Seed for drawing the pairs with `n_pairs`. + exclude_self_pairs + Leave out the correlation of each cell with itself. These self-pairs are always 1; + including them (the default, as in earlier versions) raises the upper quartile for small + groups, while for large groups the difference is negligible. + + Returns + ------- + If `inplace=True`, returns `None` and sets the field below. + Otherwise returns a dictionary with one value per group. Groups with a + single cell have no cell pairs to correlate and get a score of `NaN`. + Sets the following fields: + + `adata.obs[key_added]` : float + Score of the group each cell belongs to, `NaN` for single-cell groups. + """, +} + + @_deprecated_positional("groupby") +@_splice_doc(**_ITH_DOC) def ithgex( adata: AnnData, *, @@ -144,35 +201,21 @@ def ithgex( The score is defined as follows: - Intratumoral heterogeneity scores based on CNAs and gene expressions - The calculations of intratumoral heterogeneity scores were inspired by a - previous study and modified as follows. First, to calculate ITHCNA, we used - the relative expression value matrix generated by inferCNV and calculated the - pairwise cell–cell distances using Pearson's correlation coefficients for each - patient. ITHCNA was defined as interquartile range (IQR) of the distribution for - all malignant cell pairs' Pearson's correlation coefficients. **Similarly, we also - used gene expression profiles of cancer cells of each patient to construct the - distribution of the intratumoral distances. ITHGEX was assigned as the IQR of the - distribution.** + {wu2021} + **Similarly, we also used gene expression profiles of cancer cells of each patient + to construct the distribution of the intratumoral distances. ITHGEX was assigned as + the IQR of the distribution.** (from :cite:`Wu2021`) - The correlations are taken over all ordered pairs of cells in a group, including each - cell with itself (self-pairs, see `exclude_self_pairs`), as in earlier versions. Cells - with a constant profile, e.g. all-zero rows, have no defined correlation and are left out. + {pairs} + Cells with a constant profile, e.g. all-zero rows, have no defined correlation and are left out. - .. versionchanged:: 0.7.0 - The correlations are computed block by block with a numba kernel instead of as one - cells x cells matrix, so large groups no longer run out of memory. The quantiles match - numpy to about 1e-6. Cells with a constant profile (e.g. no CNV signal at all) have no - defined correlation and are left out; before, any such cell made the group's score `NaN`. - The number of cells left out is logged per group. - `groupby` is keyword-only; passing it positionally still works with a :class:`FutureWarning` - until 0.8.0. `use_raw` is a plain boolean, `False` by default; `None` still works as - `False` with a :class:`FutureWarning` until 0.8.0. - - .. versionadded:: 0.7.0 - `n_pairs` and `exclude_self_pairs`. + {changed} + `use_raw` is a plain boolean, `False` by default; `None` still works as `False` with a + :class:`FutureWarning` until 0.8.0. + + {added} Parameters ---------- @@ -189,30 +232,7 @@ def ithgex( If True, store the result in adata, otherwise return it. key_added Key under which the score will be stored in `adata.obs`. - n_jobs - Number of threads. Default: all cores (numba's setting, see :func:`numba.set_num_threads`). - n_pairs - Estimate the score of large groups from this many randomly drawn cell pairs instead of - all pairs, e.g. `1_000_000`. Pairs are drawn uniformly from the full correlation matrix, - so the estimate is consistent (it converges to the exact score as `n_pairs` grows); its - error shrinks like `1 / sqrt(n_pairs)` and the cost no longer grows with the group size. - Groups with fewer pairs are computed exactly. Default: all pairs (exact). - random_state - Seed for drawing the pairs with `n_pairs`. - exclude_self_pairs - Leave out the correlation of each cell with itself. These self-pairs are always 1; - including them (the default, as in earlier versions) raises the upper quartile for small - groups, while for large groups the difference is negligible. - - Returns - ------- - If `inplace=True`, returns `None` and sets the field below. - Otherwise returns a dictionary with one value per group. Groups with a - single cell have no cell pairs to correlate and get a score of `NaN`. - Sets the following fields: - - `adata.obs[key_added]` : float - Score of the group each cell belongs to, `NaN` for single-cell groups. + {params} """ if use_raw is None: warnings.warn( @@ -227,6 +247,7 @@ def ithgex( @_deprecated_positional("groupby") +@_splice_doc(**_ITH_DOC) def ithcna( adata: AnnData, *, @@ -246,35 +267,21 @@ def ithcna( The score is defined as follows: - Intratumoral heterogeneity scores based on CNAs and gene expressions - The calculations of intratumoral heterogeneity scores were inspired by a - previous study and modified as follows. First, to calculate ITHCNA, we used - the relative expression value matrix generated by inferCNV and calculated the - pairwise cell–cell distances using Pearson's correlation coefficients for each - patient. ITHCNA was defined as interquartile range (IQR) of the distribution for - all malignant cell pairs' Pearson's correlation coefficients. + {wu2021} (from :cite:`Wu2021`) - The correlations are taken over all ordered pairs of cells in a group, including each - cell with itself (self-pairs, see `exclude_self_pairs`), as in earlier versions. Cells - with a constant profile, e.g. the all-zero rows of a CNV matrix where every value lies + {pairs} + Cells with a constant profile, e.g. the all-zero rows of a CNV matrix where every value lies in the noise band of :func:`infercnvpy.tl.infercnv`, have no defined correlation and are left out. - .. versionchanged:: 0.7.0 - The correlations are computed block by block with a numba kernel instead of as one - cells x cells matrix, so large groups no longer run out of memory. The quantiles match - numpy to about 1e-6. Cells with a constant profile (e.g. no CNV signal at all) have no - defined correlation and are left out; before, any such cell made the group's score `NaN`. - The number of cells left out is logged per group. - `groupby` is keyword-only; passing it positionally still works with a :class:`FutureWarning` - until 0.8.0. `use_rep` is the suffix of the key in `adata.obsm`, `"cnv"` for - `adata.obsm["X_cnv"]`, as in the other functions; the full key (`"X_cnv"`) still works with a + {changed} + `use_rep` is the suffix of the key in `adata.obsm`, `"cnv"` for `adata.obsm["X_cnv"]`, + as in the other functions; a full key of `adata.obsm` (e.g. `"X_cnv"`) still works with a :class:`FutureWarning` until 0.8.0. - .. versionadded:: 0.7.0 - `n_pairs` and `exclude_self_pairs`. + {added} Parameters ---------- @@ -288,36 +295,16 @@ def ithcna( Key under which the score will be stored in `adata.obs`. inplace If True, store the result in adata, otherwise return it. - n_jobs - Number of threads. Default: all cores (numba's setting, see :func:`numba.set_num_threads`). - n_pairs - Estimate the score of large groups from this many randomly drawn cell pairs instead of - all pairs, e.g. `1_000_000`. Pairs are drawn uniformly from the full correlation matrix, - so the estimate is consistent (it converges to the exact score as `n_pairs` grows); its - error shrinks like `1 / sqrt(n_pairs)` and the cost no longer grows with the group size. - Groups with fewer pairs are computed exactly. Default: all pairs (exact). - random_state - Seed for drawing the pairs with `n_pairs`. - exclude_self_pairs - Leave out the correlation of each cell with itself. These self-pairs are always 1; - including them (the default, as in earlier versions) raises the upper quartile for small - groups, while for large groups the difference is negligible. - - Returns - ------- - If `inplace=True`, returns `None` and sets the field below. - Otherwise returns a dictionary with one value per group. Groups with a - single cell have no cell pairs to correlate and get a score of `NaN`. - Sets the following fields: - - `adata.obs[key_added]` : float - Score of the group each cell belongs to, `NaN` for single-cell groups. + {params} """ key = f"X_{use_rep}" - if key not in adata.obsm and use_rep.startswith("X_") and use_rep in adata.obsm: + if key not in adata.obsm and use_rep in adata.obsm: + suffix = f"pass `use_rep={use_rep.removeprefix('X_')!r}`" + if not use_rep.startswith("X_"): + suffix = f"store the matrix as `adata.obsm[{key!r}]`" warnings.warn( f"`use_rep={use_rep!r}` is the full key; ithcna now takes the suffix, as the other functions, and " - f"the full key will stop working in infercnvpy 0.8.0. Pass `use_rep={use_rep.removeprefix('X_')!r}`.", + f"the full key will stop working in infercnvpy 0.8.0. To use `adata.obsm[{use_rep!r}]`, {suffix}.", FutureWarning, skip_file_prefixes=SKIP_OWN_FRAMES, ) diff --git a/src/infercnvpy/tl/_spatial.py b/src/infercnvpy/tl/_spatial.py index d000031..6d8ba29 100644 --- a/src/infercnvpy/tl/_spatial.py +++ b/src/infercnvpy/tl/_spatial.py @@ -1,11 +1,11 @@ -"""Spatial tools: smooth a CNV matrix over a spatial neighbour graph, and regularise labels with a Potts prior. +"""Spatial tools: smooth a CNV matrix over a spatial neighbor graph, and regularize labels with a Potts prior. Both read the graph from `adata.obsp`, e.g. `obsp["spatial_connectivities"]` from :func:`squidpy.gr.spatial_neighbors_knn`; infercnvpy itself does not depend on squidpy. `_hop_csr` / `_hop_dense` compute one smoothing step ``out = S @ X`` with ``S = diag(a) + diag((1 - a) / d) @ A`` (`A` the graph without its diagonal, `d` its row sums, `a` the self weight per cell) without forming `S`: the row -normalisation and the self weight are applied on the fly while each output row gathers its neighbours' rows, so a +normalization and the self weight are applied on the fly while each output row gathers its neighbors' rows, so a CSR input is read in place and never densified. `_icm` runs sequential iterated conditional modes sweeps. """ @@ -25,14 +25,14 @@ @njit def _row_weights(g_indptr, g_indices, g_data, i, alpha): - """Self weight `a` and neighbour scale `(1 - a) / d` of cell `i` (diagonal entries of the graph are ignored).""" + """Self weight `a` and neighbor scale `(1 - a) / d` of cell `i` (diagonal entries of the graph are ignored).""" d = 0.0 deg = 0 for e in range(g_indptr[i], g_indptr[i + 1]): if g_indices[e] != i: d += g_data[e] deg += 1 - if d <= 0: # no neighbours: keep the cell's own values + if d <= 0: # no neighbors: keep the cell's own values return 1.0, 0.0 a = alpha if alpha >= 0 else 1.0 / (deg + 1) return a, (1.0 - a) / d @@ -74,7 +74,7 @@ def _hop_dense(g_indptr, g_indices, g_data, alpha, X, out): # pragma: no cover @njit def _icm(indptr, indices, weights, unary, labels, beta, n_iter): # pragma: no cover - """Iterated conditional modes for a Potts model, in place on `labels` (-1 = unlabelled). + """Iterated conditional modes for a Potts model, in place on `labels` (-1 = unlabeled). Returns the sweeps run and whether the last one changed no label. @@ -95,7 +95,7 @@ def _icm(indptr, indices, weights, unary, labels, beta, n_iter): # pragma: no c for p in range(indptr[i], indptr[i + 1]): lab = labels[indices[p]] if lab >= 0 and indices[p] != i: - # penalty beta * w for every class but the neighbour's, up to a constant per cell + # penalty beta * w for every class but the neighbor's, up to a constant per cell e[lab] -= beta * weights[p] informative = True current = labels[i] @@ -117,7 +117,7 @@ def _graph(adata: AnnData, connectivity_key: str) -> sp.csr_matrix: if connectivity_key == "distances" or connectivity_key.endswith("_distances"): raise ValueError( f"`connectivity_key={connectivity_key!r}` looks like a distance matrix. The graph weights must be " - "similarities or binary connectivities, where a larger weight means a closer neighbour, e.g. " + "similarities or binary connectivities, where a larger weight means a closer neighbor, e.g. " f"{connectivity_key.removesuffix('distances') + 'connectivities'!r}." ) if connectivity_key not in adata.obsp: @@ -144,10 +144,10 @@ def smooth_cnv( inplace: bool = True, n_jobs: int | None = None, ) -> np.ndarray | None: - r"""Smooth a CNV matrix over a spatial neighbour graph. + r"""Smooth a CNV matrix over a spatial neighbor graph. Each cell's CNV profile is replaced by a weighted mean of its own profile and those of its - neighbours in the graph `adata.obsp[connectivity_key]`. One step is + neighbors in the graph `adata.obsp[connectivity_key]`. One step is .. math:: @@ -155,28 +155,28 @@ def smooth_cnv( with the graph weights :math:`w_{ij}` (the diagonal is ignored) and the self weight :math:`a_i` (`alpha`). With `n_hops > 1` the step is repeated on its own output, which - widens the neighbourhood to cells up to `n_hops` edges away. Cells without neighbours keep + widens the neighborhood to cells up to `n_hops` edges away. Cells without neighbors keep their values. On imaging-based spatial data (Xenium, CosMx, ...) a cell has a few hundred counts from a - gene panel, so single-cell CNV profiles are noisy. Neighbouring cells in a tumour usually + gene panel, so single-cell CNV profiles are noisy. Neighboring cells in a tumor usually share the same copy number state, and averaging over them recovers the signal, as kNN smoothing in InSituCNV :cite:`Jensen2025` and spatial pooling in SlideCNA :cite:`Zhang2025` do. Han et al. :cite:`Han2026` found that aggregating cells on a grid improved CNV recovery - on imaging platforms at native resolution. Smoothing also blurs the boundary between tumour + on imaging platforms at native resolution. Smoothing also blurs the boundary between tumor and normal tissue by about `n_hops` cells. By default, :func:`infercnvpy.tl.infercnv` sets the values in its noise band to 0, and this function averages those zeros. On simulated imaging data, smoothing the unfiltered matrix of `tl.infercnv(..., dynamic_threshold=None)` recovered the CNV profiles better and - separated tumour from normal cells better than smoothing the filtered one. + separated tumor from normal cells better than smoothing the filtered one. The graph usually comes from :func:`squidpy.gr.spatial_neighbors_knn`, e.g. `sq.gr.spatial_neighbors_knn(adata, n_neighs=10)`, or from :func:`squidpy.gr.spatial_neighbors_delaunay`. Any cells x cells graph with non-negative weights works, as long as the weights are similarities or binary connectivities (a larger - weight means a closer neighbour). Distances, e.g. `obsp["spatial_distances"]`, give the - farthest neighbours the largest weight and are rejected. + weight means a closer neighbor). Distances, e.g. `obsp["spatial_distances"]`, give the + farthest neighbors the largest weight and are rejected. Parameters ---------- @@ -192,8 +192,8 @@ def smooth_cnv( The result is stored in `adata.obsm["X_{key_added}"]` and `adata.uns[key_added]`. alpha Weight of the cell itself, between 0 and 1. Default: the cell counts as one more - neighbour, i.e. `1 / (n_neighbours + 1)`; with binary connectivities this is the plain - mean over the cell and its neighbours. + neighbor, i.e. `1 / (n_neighbors + 1)`; with binary connectivities this is the plain + mean over the cell and its neighbors. n_hops Number of smoothing steps. inplace @@ -244,7 +244,7 @@ def smooth_cnv( n_isolated = int(np.count_nonzero(np.diff(G.indptr) == 0)) log_result( "smooth_cnv", - "smoothed %d cells over %d hops" + (", %d cells without neighbours kept their values" if n_isolated else ""), + "smoothed %d cells over %d hops" + (", %d cells without neighbors kept their values" if n_isolated else ""), *([X.shape[0], n_hops] + ([n_isolated] if n_isolated else [])), time=start, ) @@ -269,9 +269,9 @@ def spatial_regions( key_added: str | None = None, inplace: bool = True, ) -> pd.Series | None: - r"""Make per-cell labels spatially coherent with a Potts prior on a neighbour graph. + r"""Make per-cell labels spatially coherent with a Potts prior on a neighbor graph. - Any categorical labels work, e.g. tumour/normal calls or CNV clusters from + Any categorical labels work, e.g. tumor/normal calls or CNV clusters from :func:`infercnvpy.tl.leiden`. The new labels are a local minimum of the energy of a hidden Markov random field with a Potts prior, @@ -281,18 +281,18 @@ def spatial_regions( where :math:`w_{ij}` are the graph weights, symmetrised and scaled so that the average cell's weights sum to 1. :math:`\beta` is the penalty for disagreeing with a whole - neighbourhood. The unary term :math:`U_i` is either + neighborhood. The unary term :math:`U_i` is either - a label-flip penalty: 0 for the original label of the cell and 1 for any other label (default). For binary connectivities, two classes and a cell with the average number of - neighbours, the cell changes its label when more than :math:`(1 + 1/\beta) / 2` of its - neighbours disagree with it (75% for `beta=2`, 62.5% for `beta=4`). + neighbors, the cell changes its label when more than :math:`(1 + 1/\beta) / 2` of its + neighbors disagree with it (75% for `beta=2`, 62.5% for `beta=4`). - the negative log of per-cell class probabilities (`probabilities`). - Cells with a missing label (`NaN`) have no unary term and take the label their neighbours - favour. The energy is minimised by iterated conditional modes (ICM, :cite:`Besag1986`): + Cells with a missing label (`NaN`) have no unary term and take the label their neighbors + favor. The energy is minimized by iterated conditional modes (ICM, :cite:`Besag1986`): cells are visited in a fixed order and each takes the label with the lowest energy given - its neighbours, until a sweep changes nothing or after `n_iter` sweeps. The result is + its neighbors, until a sweep changes nothing or after `n_iter` sweeps. The result is deterministic and stays close to the input labels; it is not a new clustering. STARCH :cite:`Elyanow2021` and CalicoST :cite:`Ma2024` use the same kind of prior on clone labels. @@ -308,7 +308,7 @@ def spatial_regions( probabilities Key in `adata.obsm` with per-cell class probabilities or non-negative scores, cells x classes, in the order of the categories of `adata.obs[key]` (or a - :class:`~pandas.DataFrame` with the categories as columns). Rows are normalised to sum + :class:`~pandas.DataFrame` with the categories as columns). Rows are normalized to sum to 1. Default: use the label-flip penalty. beta Strength of the spatial prior. `0` returns the input labels (or the most probable class). @@ -326,7 +326,7 @@ def spatial_regions( Sets the following fields: `adata.obs[key_added]` : :class:`pandas.Categorical` - Spatially regularised labels, with the categories of `adata.obs[key]`. + Spatially regularized labels, with the categories of `adata.obs[key]`. `adata.uns[key_added]["params"]` : :class:`dict` Parameters of this call and the number of sweeps run (`n_sweeps`). """ diff --git a/tests/conftest.py b/tests/conftest.py index ae650ef..6d09af9 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -61,6 +61,12 @@ def _oligodendroglioma() -> AnnData: return cnv.datasets.oligodendroglioma() +@pytest.fixture(scope="session") +def oligo(_oligodendroglioma) -> AnnData: + """The oligodendroglioma dataset, for tests that leave it as it is.""" + return _oligodendroglioma + + @pytest.fixture def adata_oligodendroglioma(_oligodendroglioma) -> AnnData: """A fresh copy of the oligodendroglioma dataset with raw counts in .X (csr).""" diff --git a/tests/test_datasets.py b/tests/test_datasets.py index c24de57..fd6ad2c 100644 --- a/tests/test_datasets.py +++ b/tests/test_datasets.py @@ -24,3 +24,25 @@ def test_maynard2020_3k(tmp_path, monkeypatch): adata = cnv.datasets.maynard2020_3k() assert adata.shape == tuple(_DATASETS["maynard2020_3k"].metadata["shape"]) assert (tmp_path / "maynard2020_3k.h5ad").is_file() + + +@pytest.mark.parametrize("cached", [True, False]) +def test_maynard2020_3k_reuses_the_file_of_infercnvpy_0_6(tmp_path, monkeypatch, cached): + """infercnvpy 0.6 downloaded it into scanpy's `datasetdir`: a file there is used instead of a new download.""" + import pooch + import scanpy as sc + + import infercnvpy.datasets + + monkeypatch.delenv(DATASET_ENV_VAR, raising=False) + if cached: + (tmp_path / "maynard2020_3k.h5ad").touch() + used = [] + monkeypatch.setattr(infercnvpy.datasets, "fetch", lambda entry, cache_dir, **kw: used.append(cache_dir)) + datasetdir = sc.settings.datasetdir + sc.settings.datasetdir = tmp_path + try: + cnv.datasets.maynard2020_3k() + finally: + sc.settings.datasetdir = datasetdir + assert used == [tmp_path if cached else pooch.os_cache("infercnvpy")] diff --git a/tests/test_infercnv_kernel.py b/tests/test_infercnv_kernel.py index a36b0d0..4eebe8b 100644 --- a/tests/test_infercnv_kernel.py +++ b/tests/test_infercnv_kernel.py @@ -1,5 +1,6 @@ """The numba kernels of `tl.infercnv` against a plain numpy implementation of the algorithm.""" +import math import threading import numba @@ -154,11 +155,6 @@ def test_median(n): assert _median(values, np.empty(n)) == np.median(values) -@pytest.fixture(scope="module") -def oligo(): - return cnv.datasets.oligodendroglioma() - - REF = {"reference_key": "cell_type", "reference_cat": ["Microglia/Macrophage", "Oligodendrocytes (non-malignant)"]} @@ -280,6 +276,27 @@ def test_csr_builder(preallocated): assert res.indptr.dtype == res.indices.dtype == np.int32 +@pytest.mark.parametrize("dtype", [np.float32, np.float64]) +@pytest.mark.parametrize("expm1", [False, True]) +def test_group_sums(dtype, expm1): + """Per chunk and group, the float64 sums of the rows (or of their `expm1`) in row order, CSR and dense alike.""" + from infercnvpy.tl._infercnv_kernels import _group_sums_csr, _group_sums_dense + + rng = np.random.default_rng(0) + D = (rng.random((300, 40)) * (rng.random((300, 40)) < 0.3) * 4).astype(dtype) + rows = np.sort(rng.choice(300, 200, replace=False)).astype(np.int64) + groups = rng.integers(0, 3, rows.size).astype(np.int64) + want = np.zeros((2, 3, 40)) + for i, (r, g) in enumerate(zip(rows, groups, strict=True)): + x = D[r].astype(np.float64) + want[i * 2 // rows.size, g] += [math.expm1(v) for v in x] if expm1 else x # libm, as numba + C = sp.csr_matrix(D) + for kernel, X in ((_group_sums_csr, (C.data, C.indices, C.indptr)), (_group_sums_dense, (D,))): + sums = np.zeros_like(want) + kernel(*X, rows, groups, expm1, sums) + np.testing.assert_array_equal(sums, want) + + @pytest.fixture(scope="module") def oligo_log(oligo): """The oligodendroglioma counts, log-normalised as infercnv expects (and as `min_mean` assumes).""" diff --git a/tests/test_io.py b/tests/test_io.py index 496ad4a..1732036 100644 --- a/tests/test_io.py +++ b/tests/test_io.py @@ -44,6 +44,24 @@ def test_get_genomic_position_from_gtf_without_adata(gtf, testdata): assert not gene_pos["gene_name"].duplicated().any() +def test_gtfparse_2_polars_frame(testdata, monkeypatch): + """gtfparse 2 returns a polars DataFrame: it is converted with its `to_pandas`.""" + import gtfparse + + expected = genomic_position_from_gtf(testdata / "chr21_gencode.gtf") + + class PolarsLike: + def __init__(self, df): + self.df = df + + def to_pandas(self): + return self.df + + read_gtf = gtfparse.read_gtf + monkeypatch.setattr(gtfparse, "read_gtf", lambda *args, **kwargs: PolarsLike(read_gtf(*args, **kwargs))) + pd.testing.assert_frame_equal(genomic_position_from_gtf(testdata / "chr21_gencode.gtf"), expected) + + def _genepos_message(caplog, func): [msg] = [r.getMessage() for r in caplog.records if r.getMessage().startswith(f"{func}:")] return msg diff --git a/tests/test_lazy.py b/tests/test_lazy.py index d8d2bb7..cc8ad4d 100644 --- a/tests/test_lazy.py +++ b/tests/test_lazy.py @@ -19,23 +19,18 @@ @pytest.fixture(scope="module") -def base(): - return cnv.datasets.oligodendroglioma() - - -@pytest.fixture(scope="module") -def expected(base): +def expected(oligo): """In-memory result, with gene values.""" - return cnv.tl.infercnv(base, **REF, calculate_gene_values=True, inplace=False) + return cnv.tl.infercnv(oligo, **REF, calculate_gene_values=True, inplace=False) -def _with_X(base, X): - return ad.AnnData(X=X, obs=base.obs.copy(), var=base.var.copy()) +def _with_X(oligo, X): + return ad.AnnData(X=X, obs=oligo.obs.copy(), var=oligo.var.copy()) -def _csr_blocks(base, block_rows, counter=None): +def _csr_blocks(oligo, block_rows, counter=None): """Dask array whose row blocks are CSR matrices loaded lazily, like a lazily loaded `spatialdata_io` table.""" - X = sp.csr_matrix(base.X) + X = sp.csr_matrix(oligo.X) n, g = X.shape def load(start, stop): @@ -68,27 +63,27 @@ def _check(adata, expected, gene_values=True, **kw): @pytest.mark.parametrize("block_rows", [40, 500]) -def test_dask_dense(base, expected, block_rows): - X = da.from_array(base.X.toarray(), chunks=(block_rows, 1000)) # also split along genes - _check(_with_X(base, X), expected) +def test_dask_dense(oligo, expected, block_rows): + X = da.from_array(oligo.X.toarray(), chunks=(block_rows, 1000)) # also split along genes + _check(_with_X(oligo, X), expected) @pytest.mark.parametrize("block_rows", [30, 184]) -def test_dask_sparse_blocks(base, expected, block_rows): - _check(_with_X(base, _csr_blocks(base, block_rows)), expected) +def test_dask_sparse_blocks(oligo, expected, block_rows): + _check(_with_X(oligo, _csr_blocks(oligo, block_rows)), expected) -def test_dask_layer(base, expected): - a = _with_X(base, sp.csr_matrix(base.X.shape, dtype=np.float32)) - a.layers["raw"] = _csr_blocks(base, 50) +def test_dask_layer(oligo, expected): + a = _with_X(oligo, sp.csr_matrix(oligo.X.shape, dtype=np.float32)) + a.layers["raw"] = _csr_blocks(oligo, 50) _check(a, expected, layer="raw") @pytest.mark.parametrize("fmt", ["csr", "dense"]) -def test_backed_h5ad(base, expected, tmp_path, fmt): - X = sp.csr_matrix(base.X) if fmt == "csr" else base.X.toarray() +def test_backed_h5ad(oligo, expected, tmp_path, fmt): + X = sp.csr_matrix(oligo.X) if fmt == "csr" else oligo.X.toarray() path = tmp_path / "a.h5ad" - _with_X(base, X).write_h5ad(path) + _with_X(oligo, X).write_h5ad(path) a = ad.read_h5ad(path, backed="r") assert a.isbacked _check(a, expected, chunksize=60) @@ -96,35 +91,57 @@ def test_backed_h5ad(base, expected, tmp_path, fmt): assert isinstance(a.obsm["X_cnv"], sp.csr_matrix) +def test_backed_csc_is_read_once(oligo, expected, tmp_path, monkeypatch): + """A row block of a backed CSC matrix would read all of it: the matrix is read into memory once instead.""" + import h5py + + path = tmp_path / "a.h5ad" + _with_X(oligo, sp.csc_matrix(oligo.X)).write_h5ad(path) + a = ad.read_h5ad(path, backed="r") + read = [] + getitem = h5py.Dataset.__getitem__ + + def counting(self, *args, **kwargs): + res = getitem(self, *args, **kwargs) + if self.name == "/X/data": + read.append(np.size(res)) + return res + + monkeypatch.setattr(h5py.Dataset, "__getitem__", counting) + with pytest.warns(UserWarning, match="CSC"): + _check(a, expected, chunksize=60) + assert sum(read) == a.X.group["data"].shape[0] + + @pytest.mark.parametrize("dynamic_threshold", [None, 1.5]) -def test_reads_each_block_once_per_pass(base, expected, dynamic_threshold): +def test_reads_each_block_once_per_pass(oligo, expected, dynamic_threshold): """Reference means and the noise band (blocks with reference cells), and the result: one read of a block each.""" loads = [] - a = _with_X(base, _csr_blocks(base, 46, counter=loads)) + a = _with_X(oligo, _csr_blocks(oligo, 46, counter=loads)) cnv.tl.infercnv(a, **REF, dynamic_threshold=dynamic_threshold) - n_blocks = -(-base.n_obs // 46) - is_ref = base.obs["cell_type"].isin(REF["reference_cat"]).to_numpy() - with_ref = sum(is_ref[s : s + 46].any() for s in range(0, base.n_obs, 46)) + n_blocks = -(-oligo.n_obs // 46) + is_ref = oligo.obs["cell_type"].isin(REF["reference_cat"]).to_numpy() + with_ref = sum(is_ref[s : s + 46].any() for s in range(0, oligo.n_obs, 46)) assert 0 < with_ref < n_blocks assert len(loads) == n_blocks + with_ref * (1 if dynamic_threshold is None else 2) assert max(stop - start for start, stop in loads) <= 46, "more than one block was loaded at once" -def test_streaming_reference_matches_in_memory(base): - lazy = _with_X(base, _csr_blocks(base, 25)) +def test_streaming_reference_matches_in_memory(oligo): + lazy = _with_X(oligo, _csr_blocks(oligo, 25)) for key, cats in ((None, None), (REF["reference_key"], REF["reference_cat"])): with pytest.warns(UserWarning) if key is None else nullcontext(): got = _get_reference(lazy, key, cats, None, None) - want = np.asarray(_get_reference(base, key, cats, None, None)) + want = np.asarray(_get_reference(oligo, key, cats, None, None)) np.testing.assert_allclose(got, want, rtol=1e-5, atol=1e-7) @pytest.mark.parametrize("scheduler", ["threads", "synchronous"]) -def test_lazy_output(base, expected, scheduler): +def test_lazy_output(oligo, expected, scheduler): """With `lazy=True`, the result stays a dask array of CSR blocks and equals the eager result exactly.""" - X = _csr_blocks(base, 50) - eager = cnv.tl.infercnv(_with_X(base, X), **REF, calculate_gene_values=True, inplace=False) - a = _with_X(base, X) + X = _csr_blocks(oligo, 50) + eager = cnv.tl.infercnv(_with_X(oligo, X), **REF, calculate_gene_values=True, inplace=False) + a = _with_X(oligo, X) with dask.config.set(scheduler=scheduler): cnv.tl.infercnv(a, **REF, lazy=True, calculate_gene_values=True) assert isinstance(a.obsm["X_cnv"], da.Array) @@ -132,20 +149,40 @@ def test_lazy_output(base, expected, scheduler): res, gv = dask.compute(a.obsm["X_cnv"], a.layers["gene_values_cnv"]) assert (sp.csr_matrix(res) != eager[1]).nnz == 0 np.testing.assert_array_equal(gv, eager[2]) - _check(_with_X(base, X), expected, lazy=True) + _check(_with_X(oligo, X), expected, lazy=True) + + +@pytest.mark.parametrize("scheduler", ["synchronous", "processes"]) +def test_kernels_in_dask_tasks_run_serially(oligo, expected, scheduler, monkeypatch): + """Inside a dask task, every kernel runs its serial twin, also on the main thread (synchronous, or processes).""" + from infercnvpy._numba import njit_parallel + from infercnvpy.tl import _infercnv_kernels + + parallel_calls = [] + def spy(name, kernel): + return lambda *args: (parallel_calls.append(name), kernel.serial(*args))[1] -def test_lazy_is_opt_in(base): + for name, kernel in vars(_infercnv_kernels).items(): + if isinstance(kernel, njit_parallel): + monkeypatch.setattr(kernel, "parallel", spy(name, kernel)) + a = _with_X(oligo, _csr_blocks(oligo, 50)) + with dask.config.set(scheduler=scheduler): + _check(a, expected, lazy=True) + assert parallel_calls == [] + + +def test_lazy_is_opt_in(oligo): """Without `lazy=True`, the result of a dask array is in memory, whatever the active dask scheduler.""" - a = _with_X(base, _csr_blocks(base, 100)) + a = _with_X(oligo, _csr_blocks(oligo, 100)) with dask.config.set(scheduler="threads"): cnv.tl.infercnv(a, **REF) assert isinstance(a.obsm["X_cnv"], sp.csr_matrix) -def test_tools_that_need_cnv_in_memory(base): +def test_tools_that_need_cnv_in_memory(oligo): """Tools without block-wise reading tell the user to load a lazy CNV matrix; cnv_score reads it block by block.""" - a = _with_X(base, _csr_blocks(base, 100)) + a = _with_X(oligo, _csr_blocks(oligo, 100)) cnv.tl.infercnv(a, **REF, lazy=True) eager = a.copy() eager.obsm["X_cnv"] = eager.obsm["X_cnv"].compute() @@ -164,30 +201,33 @@ def test_tools_that_need_cnv_in_memory(base): ) -def test_lazy_without_threshold(base, expected): - res = cnv.tl.infercnv(_with_X(base, _csr_blocks(base, 50)), **REF, dynamic_threshold=None, lazy=True, inplace=False) - want = cnv.tl.infercnv(base, **REF, dynamic_threshold=None, inplace=False)[1] +def test_lazy_without_threshold(oligo, expected): + res = cnv.tl.infercnv( + _with_X(oligo, _csr_blocks(oligo, 50)), **REF, dynamic_threshold=None, lazy=True, inplace=False + ) + want = cnv.tl.infercnv(oligo, **REF, dynamic_threshold=None, inplace=False)[1] np.testing.assert_allclose(res[1].compute().toarray(), want.toarray(), rtol=1e-5, atol=1e-5) -def test_block_size_does_not_change_the_result(base): +def test_block_size_does_not_change_the_result(oligo): results = [ - cnv.tl.infercnv(_with_X(base, _csr_blocks(base, rows)), **REF, inplace=False)[1].toarray() for rows in (10, 97) + cnv.tl.infercnv(_with_X(oligo, _csr_blocks(oligo, rows)), **REF, inplace=False)[1].toarray() + for rows in (10, 97) ] np.testing.assert_allclose(results[0], results[1], rtol=1e-5, atol=1e-6) -def test_lazily_normalised_table(base): +def test_lazily_normalised_table(oligo): """A dask array of CSR blocks, normalised lazily by scanpy, gives the eager result.""" - counts = sp.csr_matrix(base.X) + counts = sp.csr_matrix(oligo.X) def normalised(X): - a = _with_X(base, X) + a = _with_X(oligo, X) sc.pp.normalize_total(a, target_sum=1e4) sc.pp.log1p(a) return a - lazy_input = _with_X(base, counts) + lazy_input = _with_X(oligo, counts) lazy_input.X = da.from_array(counts, chunks=(40, -1), asarray=False) lazy = normalised(lazy_input.X) assert isinstance(lazy.X, da.Array) @@ -197,21 +237,21 @@ def normalised(X): @pytest.fixture(scope="module") -def zarr_store(base, tmp_path_factory): +def zarr_store(oligo, tmp_path_factory): path = tmp_path_factory.mktemp("lazy") / "oligo.zarr" - b = base.copy() + b = oligo.copy() b.X = sp.csr_matrix(b.X) b.write_zarr(path) return path -def test_read_elem_lazy(base, expected, zarr_store): +def test_read_elem_lazy(oligo, expected, zarr_store): """`obs`/`var` in memory and `X` lazy via `read_elem_lazy`, as in the scanpy dask tutorial.""" import zarr from anndata.experimental import read_elem_lazy f = zarr.open(zarr_store, mode="r") - X = read_elem_lazy(f["X"], chunks=(50, base.n_vars)) + X = read_elem_lazy(f["X"], chunks=(50, oligo.n_vars)) _check(ad.AnnData(X=X, obs=ad.io.read_elem(f["obs"]), var=ad.io.read_elem(f["var"])), expected) @@ -223,11 +263,11 @@ def test_read_lazy(expected, zarr_store): assert isinstance(a.obsm["X_cnv"], sp.csr_matrix) -def test_read_lazy_downstream(base, tmp_path): +def test_read_lazy_downstream(oligo, tmp_path): """Scores and heatmaps read the groups from the lazy `obs` of `anndata.experimental.read_lazy`.""" import matplotlib.pyplot as plt - b = base.copy() + b = oligo.copy() cnv.tl.infercnv(b, **REF) b.obs["cell_type"] = b.obs["cell_type"].astype("category") b.write_zarr(tmp_path / "cnv.zarr") @@ -249,16 +289,16 @@ def test_read_lazy_downstream(base, tmp_path): np.testing.assert_array_equal(list(lazy_ith.values()), list(ith.values())) -def test_unplaced_and_extra_features(base, expected): +def test_unplaced_and_extra_features(oligo, expected): """Features without a position (e.g. control probes of a spatial panel) are skipped in lazy input too.""" n_ctrl = 20 - ctrl = sp.csr_matrix(np.random.default_rng(0).poisson(0.5, (base.n_obs, n_ctrl)).astype(np.float32)) - X = sp.hstack([sp.csr_matrix(base.X), ctrl]).tocsr() + ctrl = sp.csr_matrix(np.random.default_rng(0).poisson(0.5, (oligo.n_obs, n_ctrl)).astype(np.float32)) + X = sp.hstack([sp.csr_matrix(oligo.X), ctrl]).tocsr() var = pd.concat( - [base.var[["chromosome", "start", "end"]], pd.DataFrame({"chromosome": [None] * n_ctrl, "start": np.nan})] + [oligo.var[["chromosome", "start", "end"]], pd.DataFrame({"chromosome": [None] * n_ctrl, "start": np.nan})] ) - var.index = list(base.var_names) + [f"NegControlProbe_{i}" for i in range(n_ctrl)] - a = ad.AnnData(X=da.from_array(X, chunks=(60, -1), asarray=False), obs=base.obs.copy(), var=var) + var.index = list(oligo.var_names) + [f"NegControlProbe_{i}" for i in range(n_ctrl)] + a = ad.AnnData(X=da.from_array(X, chunks=(60, -1), asarray=False), obs=oligo.obs.copy(), var=var) chr_pos, res, gv = cnv.tl.infercnv(a, **REF, calculate_gene_values=True, inplace=False) np.testing.assert_allclose(res.toarray(), expected[1].toarray(), rtol=1e-5, atol=1e-6) assert np.isnan(gv[:, -n_ctrl:]).all() diff --git a/tests/test_numba.py b/tests/test_numba.py new file mode 100644 index 0000000..6ab478d --- /dev/null +++ b/tests/test_numba.py @@ -0,0 +1,49 @@ +"""The parallel and serial twins of `njit_parallel`, and running without numba's JIT.""" + +import os +import subprocess +import sys +import textwrap + +import numpy as np + +from infercnvpy._numba import njit_parallel, prange + + +@njit_parallel +def _inverse(x, out): + for i in prange(x.size): + out[i] = 1.0 / x[i] + + +def test_twins_divide_by_zero_alike(): + """Both twins use numpy's error model: a division by zero gives inf, not a ZeroDivisionError in one of them.""" + x = np.array([0.0, 2.0]) + for twin in (_inverse.parallel, _inverse.serial): + out = np.empty(2) + twin(x, out) + np.testing.assert_array_equal(out, [np.inf, 0.5]) + + +def test_without_jit(): + """With `NUMBA_DISABLE_JIT=1` (e.g. to debug or measure coverage), the kernels run as Python.""" + script = """ + import numpy as np, pandas as pd, scipy.sparse as sp + from anndata import AnnData + import infercnvpy as cnv + + var = pd.DataFrame({"chromosome": ["chr1"] * 12, "start": np.arange(12), "end": np.arange(12) + 1}) + X = sp.random(6, 12, density=0.5, format="csr", random_state=0, dtype=np.float32) + adata = AnnData(X=X, var=var.set_index(var.index.astype(str))) + cnv.tl.infercnv(adata, window_size=4, step=2) + assert adata.obsm["X_cnv"].shape == (6, 5) + """ + env = {**os.environ, "NUMBA_DISABLE_JIT": "1"} + res = subprocess.run( + [sys.executable, "-W", "ignore", "-c", textwrap.dedent(script)], + env=env, + capture_output=True, + text=True, + check=False, + ) + assert res.returncode == 0, res.stderr diff --git a/tests/test_plotting.py b/tests/test_plotting.py index f5eddb3..80b5ad7 100644 --- a/tests/test_plotting.py +++ b/tests/test_plotting.py @@ -1,6 +1,5 @@ import warnings -import dask import dask.array as da import numpy as np import pandas as pd @@ -49,17 +48,7 @@ def _image(axes) -> np.ndarray: def _dask(X, block_rows=7): """A dask array with row blocks of `block_rows` (CSR blocks for sparse `X`).""" - if not sp.issparse(X): - return da.from_array(X, chunks=(block_rows, -1)) - X = sp.csr_matrix(X) - meta = sp.csr_matrix((0, X.shape[1]), dtype=X.dtype) - blocks = [ - da.from_delayed( - dask.delayed(X[s : s + block_rows]), shape=X[s : s + block_rows].shape, dtype=X.dtype, meta=meta - ) - for s in range(0, X.shape[0], block_rows) - ] - return da.concatenate(blocks, axis=0) + return da.from_array(sp.csr_matrix(X) if sp.issparse(X) else X, chunks=(block_rows, -1), asarray=False) # --- data --------------------------------------------------------------------------------------------------------- diff --git a/tests/test_scores.py b/tests/test_scores.py index 1df9f08..367f815 100644 --- a/tests/test_scores.py +++ b/tests/test_scores.py @@ -2,6 +2,7 @@ import dask.array as da import numpy as np +import pandas as pd import pytest import scipy.sparse as sp @@ -115,6 +116,24 @@ def test_ithcna_use_rep_full_key_is_deprecated(adata_ithgex): assert cnv.tl.ithcna(adata_ithgex, groupby="group", use_rep="X_cnv", inplace=False) == expected +def test_ithcna_use_rep_any_obsm_key_is_deprecated(adata_ithgex): + """As in 0.6, `use_rep` can still name an obsm key without the `X_` prefix.""" + expected = cnv.tl.ithcna(adata_ithgex, groupby="group", inplace=False) + adata_ithgex.obsm["custom"] = adata_ithgex.obsm["X_cnv"] + with pytest.warns(FutureWarning, match="use_rep='custom'"): + assert cnv.tl.ithcna(adata_ithgex, groupby="group", use_rep="custom", inplace=False) == expected + + +def test_dataframe_cnv_matrix(adata_ithgex): + """A CNV matrix stored as a DataFrame in obsm (as 0.6 accepted) works as its values.""" + funcs = (cnv.tl.ithcna, cnv.tl.cnv_score) + expected = [f(adata_ithgex, groupby="group", inplace=False) for f in funcs] + X = adata_ithgex.obsm["X_cnv"] + adata_ithgex.obsm["X_cnv"] = pd.DataFrame(X.toarray() if sp.issparse(X) else X, index=adata_ithgex.obs_names) + for f, want in zip(funcs, expected, strict=True): + assert f(adata_ithgex, groupby="group", inplace=False) == pytest.approx(want) + + @pytest.mark.parametrize("func", [cnv.tl.cnv_score, cnv.tl.ithgex, cnv.tl.ithcna]) def test_positional_groupby_is_deprecated(adata_ithgex, func): expected = func(adata_ithgex, groupby="group", inplace=False) diff --git a/tests/test_spatial.py b/tests/test_spatial.py index b4578c0..0a176a8 100644 --- a/tests/test_spatial.py +++ b/tests/test_spatial.py @@ -71,7 +71,7 @@ def test_serial_twin_matches(graph_and_cnv): for fn, x in ((_hop_csr, (Xc.indptr, Xc.indices, Xc.data)), (_hop_dense, (X,))): par, ser = np.empty_like(X), np.empty_like(X) fn(*args, *x, par) - fn(*args, *x, ser, serial=True) + fn.serial(*args, *x, ser) np.testing.assert_allclose(par, ser, rtol=1e-6)