From d8c20d7481bb13dd39dff0e7f0e055f54e3cc9a5 Mon Sep 17 00:00:00 2001 From: giovp Date: Fri, 9 Apr 2021 22:46:50 +0200 Subject: [PATCH] add moran from scanpy (#317) * add moran from scanpy * add abs to pvalue * fix pvalues * add standard approx from perm pvalue (pysal dep return) * change req * update tests and req * add gearys' and unify function * add analytic pvalues * fix tests and docstring * explicilty add leiden * explicilty add leiden * Fix tests * resolve comments * more fixes * Add CI job * Pass mode instead of function * Force restore tox cache, fix not passing env var * Pass backend as env * Fix doc linting condition * remove pysal from setup * remaining comments, modify return excluding nan * fix tests * Fix docs, mypy, and some pet peeves * Fix wrong requirement * [ci skip] Remove doc linting Co-authored-by: Michal Klein --- .github/workflows/ci.yml | 4 +- .mypy.ini | 2 +- .pre-commit-config.yaml | 12 +- docs/source/api.rst | 2 +- docs/source/conf.py | 2 +- docs/source/utils.py | 4 + requirements.txt | 3 +- setup.py | 2 +- squidpy/_constants/_constants.py | 6 + squidpy/_utils.py | 2 +- squidpy/gr/__init__.py | 2 +- squidpy/gr/_ppatterns.py | 299 +++++++++++++----- squidpy/im/_container.py | 4 +- squidpy/im/_feature_mixin.py | 6 +- squidpy/im/_utils.py | 4 +- tests/_images/Graph_interaction.png | Bin 1064 -> 1338 bytes tests/_images/Graph_interaction_dendro.png | Bin 1282 -> 1252 bytes .../_images/Graph_nhood_enrichment_dendro.png | Bin 1279 -> 1252 bytes tests/conftest.py | 2 - tests/graph/test_ppatterns.py | 60 ++-- tests/plotting/test_interactive.py | 2 +- tox.ini | 8 +- 22 files changed, 299 insertions(+), 127 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 20f12ce0..1ec78dcb 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -112,14 +112,12 @@ jobs: with: path: .tox key: tox-${{ runner.os }}-${{ env.pythonLocation }}-${{ hashFiles('**/tox.ini') }} - restore-keys: | - tox-${{ runner.os }}-${{ env.pythonLocation }}- - tox-${{ runner.os }}- - name: Testing run: | tox -vv env: + MPLBACKEND: agg PLATFORM: ${{ matrix.os }} DISPLAY: :42 diff --git a/.mypy.ini b/.mypy.ini index c2c41bb4..2ce1a9a9 100644 --- a/.mypy.ini +++ b/.mypy.ini @@ -14,7 +14,7 @@ disallow_any_generics = True strict_optional = True strict_equality = True warn_return_any = True -warn_unreachable = True +warn_unreachable = False check_untyped_defs = True ; because of docrep allow_untyped_decorators = True diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index afc9851b..33e0a26b 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -7,7 +7,7 @@ default_stages: minimum_pre_commit_version: 2.9.3 repos: - repo: https://github.com/pre-commit/mirrors-mypy - rev: v0.800 + rev: v0.812 hooks: - id: mypy additional_dependencies: [numpy, scipy, pandas] @@ -17,7 +17,7 @@ repos: - id: black additional_dependencies: [toml] - repo: https://github.com/timothycrosley/isort - rev: 5.7.0 + rev: 5.8.0 hooks: - id: isort additional_dependencies: [toml] @@ -61,12 +61,12 @@ repos: - id: autoflake args: [--in-place, --remove-all-unused-imports, --remove-unused-variable, --ignore-init-module-imports] - repo: https://gitlab.com/pycqa/flake8 - rev: 3.8.4 + rev: 3.9.0 hooks: - id: flake8 additional_dependencies: [flake8-tidy-imports, flake8-docstrings, flake8-rst-docstrings, flake8-comprehensions, flake8-bugbear, flake8-logging-format, flake8-blind-except, flake8-builtins, flake8-pytest-style, flake8-mock, flake8-string-format] - repo: https://github.com/jumanjihouse/pre-commit-hooks - rev: 2.1.4 + rev: 2.1.5 hooks: - id: script-must-have-extension name: Check executable files use .sh extension @@ -77,7 +77,7 @@ repos: hooks: - id: rstcheck - repo: https://github.com/asottile/blacken-docs - rev: v1.9.2 + rev: v1.10.0 hooks: - id: blacken-docs additional_dependencies: [black==20.8b1] @@ -87,7 +87,7 @@ repos: - id: pyupgrade args: [--py3-plus, --py37-plus] - repo: https://github.com/pre-commit/pygrep-hooks - rev: v1.7.0 + rev: v1.8.0 hooks: - id: python-no-eval - id: python-use-type-annotations diff --git a/docs/source/api.rst b/docs/source/api.rst index 0bda4c6c..d3cffa3c 100644 --- a/docs/source/api.rst +++ b/docs/source/api.rst @@ -19,7 +19,7 @@ Graph gr.centrality_scores gr.interaction_matrix gr.ligrec - gr.moran + gr.spatial_autocorr gr.ripley_k gr.co_occurrence diff --git a/docs/source/conf.py b/docs/source/conf.py index 98d6ef18..15305e9b 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -67,7 +67,7 @@ pandas=("https://pandas.pydata.org/pandas-docs/stable/", None), anndata=("https://anndata.readthedocs.io/en/stable/", None), scanpy=("https://scanpy.readthedocs.io/en/stable/", None), - matplotlib=("https://matplotlib.org/", None), + matplotlib=("https://matplotlib.org/stable/", None), seaborn=("https://seaborn.pydata.org/", None), joblib=("https://joblib.readthedocs.io/en/latest/", None), networkx=("https://networkx.org/documentation/stable/", None), diff --git a/docs/source/utils.py b/docs/source/utils.py index 4e26ba7b..c71ec2da 100644 --- a/docs/source/utils.py +++ b/docs/source/utils.py @@ -93,6 +93,10 @@ def _download_dir(url: str, *, path: Union[str, Path], depth: int) -> Tuple[bool def _download_notebooks(org: str, repo: str, raise_exc: bool = False) -> None: + if not int(os.environ.get("SQUIDPY_DOWNLOAD_NOTEBOOKS", 1)): + info("Not downloading notebooks because a flag is set") + return + ep = ENDPOINT_FMT.format(org=org, repo=repo) for path in [FIXED_TUTORIALS_DIR, TUTORIALS_DIR, EXAMPLES_DIR, GENMOD_DIR]: ok, reason = _download_dir(urljoin(ep, path), path=path, depth=DEPTH) diff --git a/requirements.txt b/requirements.txt index 76b06e86..b77c97df 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,10 +2,11 @@ anndata>=0.7.4 dask>=2.30 docrep>=0.3.1 ipywidgets>=7.5.1 # progress bar +leidenalg>=0.8.2 omnipath>=1.0.4 pandas>=1.2.0 rasterio>=1.1.8 -scanpy[leiden]>=1.7.0 +scanpy @ git+https://github.com/theislab/scanpy@master # scanpy[leiden]>=1.7.0 scikit-image>=0.17.1 statsmodels>=0.12.0 tifffile>=2020.10.1 diff --git a/setup.py b/setup.py index 8ddd853f..668630ac 100644 --- a/setup.py +++ b/setup.py @@ -45,7 +45,7 @@ if not l.startswith("-r") ], interactive=["PyQt5>=5.15.0", "napari>=0.4.2"], - all=["PyQt5>=5.15.0", "napari>=0.4.2", "esda>=2.3.1", "libpysal>=4.3.0", "astropy>=4.1"], + all=["PyQt5>=5.15.0", "napari>=0.4.2", "astropy>=4.1"], ), classifiers=[ "Development Status :: 5 - Production/Stable", diff --git a/squidpy/_constants/_constants.py b/squidpy/_constants/_constants.py index b968f3cd..5a2956cf 100644 --- a/squidpy/_constants/_constants.py +++ b/squidpy/_constants/_constants.py @@ -88,3 +88,9 @@ class DendrogramAxis(ModeEnum): # noqa: D101 class Symbol(ModeEnum): # noqa: D101 DISC = "disc" SQUARE = "square" + + +@unique +class SpatialAutocorr(ModeEnum): # noqa: D101 + MORAN = "moran" + GEARY = "geary" diff --git a/squidpy/_utils.py b/squidpy/_utils.py index f856d31c..1a9bf4cd 100644 --- a/squidpy/_utils.py +++ b/squidpy/_utils.py @@ -57,7 +57,7 @@ def _unique_order_preserving(iterable: Iterable[Hashable]) -> Tuple[List[Hashabl class Signal(Enum): - """Signalling values when informing parallelizer.""" + """Signaling values when informing parallelizer.""" NONE = 0 UPDATE = 1 diff --git a/squidpy/gr/__init__.py b/squidpy/gr/__init__.py index f9bea64b..18efcd01 100644 --- a/squidpy/gr/__init__.py +++ b/squidpy/gr/__init__.py @@ -2,4 +2,4 @@ from squidpy.gr._build import spatial_neighbors from squidpy.gr._nhood import nhood_enrichment, centrality_scores, interaction_matrix from squidpy.gr._ligrec import ligrec -from squidpy.gr._ppatterns import moran, ripley_k, co_occurrence +from squidpy.gr._ppatterns import ripley_k, co_occurrence, spatial_autocorr diff --git a/squidpy/gr/_ppatterns.py b/squidpy/gr/_ppatterns.py index 33a820c0..9ff27fad 100644 --- a/squidpy/gr/_ppatterns.py +++ b/squidpy/gr/_ppatterns.py @@ -1,15 +1,20 @@ """Functions for point patterns spatial statistics.""" -from typing import Tuple, Union, Iterable, Optional, Sequence +from typing import Any, Dict, Tuple, Union, Iterable, Optional, Sequence, TYPE_CHECKING from itertools import chain from typing_extensions import Literal # < 3.8 -import warnings from scanpy import logging as logg from anndata import AnnData +from scanpy.get import _get_obs_rep +from scanpy.metrics._gearys_c import _gearys_c +from scanpy.metrics._morans_i import _morans_i from numba import njit -from scipy.sparse import isspmatrix_lil +from scipy import stats +from numpy.random import default_rng +from scipy.sparse import spmatrix from sklearn.metrics import pairwise_distances +from sklearn.preprocessing import normalize from statsmodels.stats.multitest import multipletests import numpy as np import pandas as pd @@ -25,21 +30,10 @@ _assert_connectivity_key, _assert_non_empty_sequence, ) +from squidpy._constants._constants import SpatialAutocorr from squidpy._constants._pkg_constants import Key -try: - with warnings.catch_warnings(): - warnings.simplefilter("ignore", UserWarning) - from libpysal.weights import W - import esda - import libpysal -except ImportError: - esda = None - libpysal = None - W = None - - -__all__ = ["ripley_k", "moran", "co_occurrence"] +__all__ = ["ripley_k", "spatial_autocorr", "co_occurrence"] it = nt.int32 @@ -132,37 +126,52 @@ def ripley_k( @d.dedent -@inject_docs(key=Key.obsp.spatial_conn()) -def moran( +@inject_docs(key=Key.obsp.spatial_conn(), sp=SpatialAutocorr) +def spatial_autocorr( adata: AnnData, connectivity_key: str = Key.obsp.spatial_conn(), genes: Optional[Union[str, Sequence[str]]] = None, - transformation: Literal["r", "B", "D", "U", "V"] = "r", - n_perms: int = 1000, + mode: Literal["moran", "geary"] = SpatialAutocorr.MORAN.s, # type: ignore[assignment] + transformation: bool = True, + n_perms: Optional[int] = None, + two_tailed: bool = False, corr_method: Optional[str] = "fdr_bh", layer: Optional[str] = None, seed: Optional[int] = None, + use_raw: bool = False, copy: bool = False, n_jobs: Optional[int] = None, backend: str = "loky", show_progress_bar: bool = True, ) -> Optional[pd.DataFrame]: """ - Calculate Moran’s I Global Autocorrelation Statistic. + Calculate Global Autocorrelation Statistic (Moran’s I or Geary's C). + + See :cite:`pysal` for reference. Parameters ---------- %(adata)s %(conn_key)s genes - List of gene names, as stored in :attr:`anndata.AnnData.var_names`, used to compute Moran's I statistics - :cite:`pysal`. + List of gene names, as stored in :attr:`anndata.AnnData.var_names`, used to compute global + spatial autocorrelation statistic. If `None`, it's computed :attr:`anndata.AnnData.var` ``['highly_variable']``, if present. Otherwise, it's computed for all genes. + mode + Mode of score calculation: + + - `{sp.MORAN.s!r}` - `Moran's I autocorrelation `_. + - `{sp.GEARY.s!r}` - `Geary's C autocorrelation `_. + transformation - Transformation to be used, as reported in :class:`esda.Moran`. Default is `"r"`, row-standardized. + If `True`, weights in :attr:`anndata.AnnData.obsp` ``['{key}']`` are row-normalized, + advised for analytic p-value calculation. %(n_perms)s + If `None`, only p-values under normality assumption are computed. + two_tailed + If `True`, p-values are two-tailed, otherwise they are one-tailed. %(corr_method)s layer Layer in :attr:`anndata.AnnData.layers` to use. If `None`, use :attr:`anndata.AnnData.X`. @@ -174,19 +183,22 @@ def moran( ------- If ``copy = True``, returns a :class:`pandas.DataFrame` with the following keys: - - `'I'` - Moran's I statistic. + - `'I' or 'C'` - Moran's I or Geary's C statistic. + - `'pval_norm'` - p-value under normality assumption. + - `'var_norm'` - variance of `'score'` under normality assumption. + - `'{{p_val}}_{{corr_method}}'` - the corrected p-values if ``corr_method != None`` . + + If ``n_perms != None`` is not None, additionally returns the following columns: + + - `'pval_z_sim'` - p-value based on standard normal approximation from permutations. - `'pval_sim'` - p-value based on permutations. - - `'VI_sim'` - variance of `'I'` from permutations. - - `'pval_sim_{{corr_method}}'` - the corrected p-values if ``corr_method != None`` . + - `'var_sim'` - variance of `'score'` from permutations. Otherwise, modifies the ``adata`` with the following key: - - :attr:`anndata.AnnData.uns` ``['moranI']`` - the above mentioned dataframe. + - :attr:`anndata.AnnData.uns` ``['moranI']`` - the above mentioned dataframe, if ``mode = {sp.MORAN.s!r}``. + - :attr:`anndata.AnnData.uns` ``['gearyC']`` - the above mentioned dataframe, if ``mode = {sp.GEARY.s!r}``. """ - if esda is None or libpysal is None: - raise ImportError("Please install `esda` and `libpysal` as `pip install esda libpysal`.") - - _assert_positive(n_perms, name="n_perms") _assert_connectivity_key(adata, connectivity_key) if genes is None: @@ -196,51 +208,89 @@ def moran( genes = adata.var_names.values genes = _assert_non_empty_sequence(genes, name="genes") + mode = SpatialAutocorr(mode) # type: ignore[assignment] + if TYPE_CHECKING: + assert isinstance(mode, SpatialAutocorr) + params = {"mode": mode.s, "transformation": transformation, "two_tailed": two_tailed} + + if mode == SpatialAutocorr.MORAN: + params["func"] = _morans_i + params["stat"] = "I" + params["expected"] = -1.0 / (adata.shape[0] - 1) # expected score + params["ascending"] = False + elif mode == SpatialAutocorr.GEARY: + params["func"] = _gearys_c + params["stat"] = "C" + params["expected"] = 1.0 + params["ascending"] = True + else: + raise NotImplementedError(f"Mode `{mode}` is not yet implemented.") + n_jobs = _get_n_cores(n_jobs) - start = logg.info(f"Calculating for `{len(genes)}` genes using `{n_jobs}` core(s)") - - w = _set_weight_class(adata, key=connectivity_key) # init weights - df = parallelize( - _moran_helper, - collection=genes, - extractor=pd.concat, - use_ixs=True, - n_jobs=n_jobs, - backend=backend, - show_progress_bar=show_progress_bar, - )(adata=adata, weights=w, transformation=transformation, permutations=n_perms, layer=layer, seed=seed) + + vals = _get_obs_rep(adata[:, genes], use_raw=use_raw, layer=layer).T + g = adata.obsp[connectivity_key].copy() + # row-normalize + if transformation: + normalize(g, norm="l1", axis=1, copy=False) + + score = params["func"](g, vals) + + start = logg.info(f"Calculating {mode}'s statistic for `{n_perms}` permutations using `{n_jobs}` core(s)") + if n_perms is not None: + _assert_positive(n_perms, name="n_perms") + perms = np.arange(n_perms) + + score_perms = parallelize( + _score_helper, + collection=perms, + extractor=np.concatenate, + use_ixs=True, + n_jobs=n_jobs, + backend=backend, + show_progress_bar=show_progress_bar, + )(mode=mode, g=g, vals=vals, seed=seed) + else: + score_perms = None + + with np.errstate(divide="ignore"): + pval_results = _p_value_calc(score, score_perms, g, params) + + results = {params["stat"]: score} + results.update(pval_results) + + df = pd.DataFrame(results, index=genes) if corr_method is not None: - _, pvals_adj, _, _ = multipletests(df["pval_sim"].values, alpha=0.05, method=corr_method) - df[f"pval_sim_{corr_method}"] = pvals_adj + for pv in filter(lambda x: "pval" in x, df.columns): + _, pvals_adj, _, _ = multipletests(df[pv].values, alpha=0.05, method=corr_method) + df[f"{pv}_{corr_method}"] = pvals_adj - df.sort_values(by="I", ascending=False, inplace=True) + df.sort_values(by=params["stat"], ascending=params["ascending"], inplace=True) if copy: logg.info("Finish", time=start) return df - _save_data(adata, attr="uns", key="moranI", data=df, time=start) + _save_data(adata, attr="uns", key=params["mode"] + params["stat"], data=df, time=start) -def _moran_helper( +def _score_helper( ix: int, - gen: Iterable[str], - adata: AnnData, - weights: W, - transformation: Literal["r", "B", "D", "U", "V"] = "r", - permutations: int = 1000, - layer: Optional[str] = None, + perms: Sequence[int], + mode: SpatialAutocorr, + g: spmatrix, + vals: np.ndarray, seed: Optional[int] = None, queue: Optional[SigQueue] = None, ) -> pd.DataFrame: - if seed is not None: - np.random.seed(seed + ix) + score_perms = np.empty((len(perms), vals.shape[0])) + rng = default_rng(None if seed is None else ix + seed) + func = _morans_i if mode == SpatialAutocorr.MORAN else _gearys_c - moran_list = [] - for g in gen: - mi = _compute_moran(adata.obs_vector(g, layer=layer), weights, transformation, permutations) - moran_list.append(mi) + for i in range(len(perms)): + idx_shuffle = rng.permutation(g.shape[0]) + score_perms[i, :] = func(g[idx_shuffle, :], vals) if queue is not None: queue.put(Signal.UPDATE) @@ -248,23 +298,7 @@ def _moran_helper( if queue is not None: queue.put(Signal.FINISH) - return pd.DataFrame(moran_list, columns=["I", "pval_sim", "VI_sim"], index=gen) - - -def _compute_moran(y: np.ndarray, w: W, transformation: str, permutations: int) -> Tuple[float, float, float]: - mi = esda.moran.Moran(y, w, transformation=transformation, permutations=permutations) - return mi.I, mi.p_z_sim, mi.VI_sim - - -def _set_weight_class(adata: AnnData, key: str) -> W: - X = adata.obsp[key] - if not isspmatrix_lil(X): - X = X.tolil() - - neighbors = dict(enumerate(X.rows)) - weights = dict(enumerate(X.data)) - - return libpysal.weights.W(neighbors, weights, ids=adata.obs.index.values) + return score_perms @njit( @@ -471,3 +505,112 @@ def _find_min_max(spatial: np.ndarray) -> Tuple[float, float]: ].astype(fp) return thres_min, thres_max + + +def _p_value_calc( + score: np.ndarray, + sims: Optional[np.ndarray], + weights: Union[spmatrix, np.ndarray], + params: Dict[str, Any], +) -> Dict[str, Any]: + """ + Handle p-value calculation for spatial autocorrelation function. + + Parameters + ---------- + score + (n_features,). + sims + (n_simulations, n_features). + params + Object to store relevant function parameters. + + Returns + ------- + pval_norm + p-value under normality assumption + pval_sim + p-values based on permutations + pval_z_sim + p-values based on standard normal approximation from permutations + """ + p_norm, var_norm = _analytic_pval(score, weights, params) + results = {"pval_norm": p_norm, "var_norm": var_norm} + + if sims is None: + return results + + n_perms = sims.shape[0] + large_perm = (sims >= score).sum(axis=0) + # subtract total perm for negative values + large_perm[(n_perms - large_perm) < large_perm] = n_perms - large_perm[(n_perms - large_perm) < large_perm] + # get p-value based on permutation + p_sim: np.ndarray = (large_perm + 1) / (n_perms + 1) + + # get p-value based on standard normal approximation from permutations + e_score_sim = sims.sum(axis=0) / n_perms + se_score_sim = sims.std(axis=0) + z_sim = (score - e_score_sim) / se_score_sim + p_z_sim = np.empty(z_sim.shape) + + p_z_sim[z_sim > 0] = 1 - stats.norm.cdf(z_sim[z_sim > 0]) + p_z_sim[z_sim <= 0] = stats.norm.cdf(z_sim[z_sim <= 0]) + + var_sim = np.var(sims, axis=0) + + results["pval_z_sim"] = p_z_sim + results["pval_sim"] = p_sim + results["var_sim"] = var_sim + + return results + + +def _analytic_pval( + score: np.ndarray, g: Union[spmatrix, np.ndarray], params: Dict[str, Any] +) -> Tuple[np.ndarray, float]: + """ + Analytic p-value computation. + + See `Moran's I `_ and + `Geary's C `_ implementation. + """ + s0, s1, s2 = _g_moments(g) + n = g.shape[0] + s02 = s0 * s0 + n2 = n * n + v_num = n2 * s1 - n * s2 + 3 * s02 + v_den = (n - 1) * (n + 1) * s02 + + Vscore_norm = v_num / v_den - (1.0 / (n - 1)) ** 2 + seScore_norm = Vscore_norm ** (1 / 2.0) + + z_norm = (score - params["expected"]) / seScore_norm + p_norm = np.empty(score.shape) + p_norm[z_norm > 0] = 1 - stats.norm.cdf(z_norm[z_norm > 0]) + p_norm[z_norm <= 0] = stats.norm.cdf(z_norm[z_norm <= 0]) + + if params["two_tailed"]: + p_norm *= 2.0 + + return p_norm, Vscore_norm + + +def _g_moments(w: Union[spmatrix, np.ndarray]) -> Tuple[np.float_, np.float_, np.float_]: + """ + Compute moments of adjacency matrix for analytic p-value calculation. + + See `pysal `_ implementation. + """ + # s0 + s0 = w.sum() + + # s1 + t = w.transpose() + w + t2 = t.multiply(t) + s1 = t2.sum() / 2.0 + + # s2 + s2array = np.array(w.sum(1) + w.sum(0).transpose()) ** 2 + s2 = s2array.sum() + + return s0, s1, s2 diff --git a/squidpy/im/_container.py b/squidpy/im/_container.py index 771cf908..b2a59788 100644 --- a/squidpy/im/_container.py +++ b/squidpy/im/_container.py @@ -930,10 +930,10 @@ def _get_size(self, size: Optional[Union[FoI_t, Tuple[Optional[FoI_t], Optional[ res = list(size) if size[0] is None: - res[0] = self.shape[0] # type: ignore[unreachable] + res[0] = self.shape[0] if size[1] is None: - res[1] = self.shape[1] # type: ignore[unreachable] + res[1] = self.shape[1] return tuple(res) # type: ignore[return-value] diff --git a/squidpy/im/_feature_mixin.py b/squidpy/im/_feature_mixin.py index 2afcf0a9..d6855cdb 100644 --- a/squidpy/im/_feature_mixin.py +++ b/squidpy/im/_feature_mixin.py @@ -271,7 +271,7 @@ def features_segmentation( Only relevant for features that use the ``intensity_layer``. props Segmentation features that are calculated. See `properties` in :func:`skimage.measure.regionprops_table`. - Each feature is calculated for each segment (e.g., nucleous) and mean and std values are returned, except + Each feature is calculated for each segment (e.g., nucleus) and mean and std values are returned, except for `'centroid'` and `'label'`. Valid options are: - `'area'` - number of pixels of segment. @@ -280,7 +280,7 @@ def features_segmentation( - `'convex_area'` - number of pixels in convex hull of segment. - `'eccentricity'` - eccentricity of ellipse with same second moments as segment. - `'equivalent_diameter'` - diameter of circles with same area as segment. - - `'euler_number'` - euler characteristic of segment. + - `'euler_number'` - Euler characteristic of segment. - `'extent'` - ratio of pixels in segment to its bounding box. - `'feret_diameter_max'` - longest distance between points around convex hull of segment. - `'filled_area'` - number of pixels of segment with all holes filled in. @@ -292,7 +292,7 @@ def features_segmentation( - `'minor_axis_length'` - length of minor axis of ellipse with same second moments as segment. - `'orientation'` - angle of major axis of ellipse with same second moments as segment. - `'perimeter'` - perimeter of segment using 4-connectivity. - - `'perimeter_crofton'` - perimeter of segmeent approximated by the Crofton formula. + - `'perimeter_crofton'` - perimeter of segment approximated by the Crofton formula. - `'solidity'` - ratio of pixels in the segment to the convex hull of the segment. Returns diff --git a/squidpy/im/_utils.py b/squidpy/im/_utils.py index 2d1ef89e..8463bc31 100644 --- a/squidpy/im/_utils.py +++ b/squidpy/im/_utils.py @@ -104,7 +104,7 @@ def to_tuple(self) -> Tuple[float, float, float, float]: def __add__(self, other: "CropPadding") -> "CropCoords": if not isinstance(other, CropPadding): - return NotImplemented # type: ignore[unreachable] + return NotImplemented return CropCoords( x0=self.x0 - other.x_pre, y0=self.y0 - other.y_pre, x1=self.x1 + other.x_post, y1=self.y1 + other.y_post @@ -112,7 +112,7 @@ def __add__(self, other: "CropPadding") -> "CropCoords": def __sub__(self, other: "CropCoords") -> "CropPadding": if not isinstance(other, CropCoords): - return NotImplemented # type: ignore[unreachable] + return NotImplemented return CropPadding( x_pre=abs(self.x0 - other.x0), diff --git a/tests/_images/Graph_interaction.png b/tests/_images/Graph_interaction.png index 26e9b7984ce8eed7648d1b66839df5169b466d50..c317a1118a5452423b3a511723ff66579f1efa60 100644 GIT binary patch delta 1284 zcmV+f1^fD_2)YW8Mt?LeF)ScxbaZfYIxjD6VRUe8Z***FVlHoTXD`{*51;@5010qN zS#tmY1}6Xj1}6bcRM^J=00gH=L_t(|obBDgOB`1m2k_tMvLc1jCA1XV(7Uh{JaoHA z7K1}~upp(ig}EY*sqUOlB!mUii7bAK>+)Kr_(hLTcE*mdW9 zdk8x_8;2BNIpr4q)*#!xPoyFM3&A*QCLFflQKg@uJ^%%!CzjE;_Ca&i&?u)e;I z>FH_g?d@THet#aPr>8hNIzp*bLZMJVxm@m=%k#(2&v|uq6_rY*>vJnBD<~8SSYBR^ z&L#a^0NC8z#Magpc6WE9F=uCIad>#>?RWft0st_f)z|9@NV+odFB*-8)oQg~%|%|) z$jFEd3=AZGuPeZ{TFq|W-m&bvA9Mq9v!!1*=#l~o6TAP0M%;MZr1@5WW@gBQJm(E@Pu;23eycLT@?V=eu zbJ4K^Oy~j%FtyM+wJQafkc%BFz=SSh#b*VWP-sDPtN;^oP=Kjj#KPKFtyb;e;A#~M z;lSQYVl2$ z0!*l-4CbZ`CKswy$t~0p%`IN4jGFG+eapA*ex6E$+~}=zdW|34@!Q09)d1n8b6TRw zwPKa5%nLpzS5$~VwP*PSY4HTT3VvH#WT3m zRe$+Mk7gs@=^9MVxrC=l!o#T+%}q-*Ij6p8ZpvVCP=L8PrvOunzaLrwCKO+yqx<0Y zTU&EKJxt~GfNy&DV_T5%TY3mw|D|)#5=~C8QR1duG`pm0un>W%7R@a7vRMU~QJ`8h zIqjm!DZu2Mx+=LTgSiEXU6p_C^=LZcoqsZzobH3U1h+UO(;+;np>c2FgciN z(bTlb;-+0RyQE!o`~TLj*$1zEK5r7VT7Jh%cvq!1S?sh#bJI1LodV2F*I;siYSH9W zi{_?H7CWa}G_|e_R#zni&(&&jA258V#$}YF4k;EtAPuE|=@YT;(O4pP$?C@URUH4Vjr) uJC@C6)3Vv@vjF#aj*=o>gUKntMt?IdGb|uzbaZfYIxjD6VRUe8Z***FVlHoTXD`T?x10a~010qN zS#tmY1}6Xj1}6bcRM^J=00WgtL_t(|obBDQYg17e1@QCIkTfW2EyY^JDoz%P(#Z^6 z3yOb02Uq_C_aX>(?XDmaI=Z>K6q@Z-Xjd~-(3Uu)#575hdw)NN7)eMM>AeXa&TmoQ zaNmWGbKkwMU!RGY8H845(i4!lvV;FpYf;ACEbT=Zw=!FbGVWpH0C*1+i^ZYq*}>HA?CeCjxr|4(TK)U? zeqEE@-CgUawtszG!P6@ru==){8asz4SGXS7m zE@S)Mlim`Yn_EQw)&l^*>cawFJiiw7<^XW~GwiZIsPp!2w6SCHXSJZ?)m(Mt+MU~3 z?V_1^comqLGgEHnnR)jGkB7qF^O9jTB51xMF%;jXD~NsIwnzbHZ|NA z4fMJnV>HNpYs~~7v#S_*HBoTN3(N;%UiW5-+okU>#)r9blCA+ZseL0 z%+Bi>Oic-9r(HBTC77HN%+0(KOf6A@$*I8PynnhXxhcWi5(p0&~-=N^)AF*{Q(Xlwfuq>Z)X>1hX@BRZ>%e+36h2En}tvQ&WPu z>3(CVY-P5D?A0PFCuP|8m+)t`2~X|cg#NH+U+C7> zvC|UGP0wI*IwrBxyJ&7&qS>jya)&mS}R$sle2fU~Z;O7CG;GO?QWl18MaCJheJ%t{VPc=X^^sl=J@g>kwkjIkaCm zK>?SBt~Fg9NvqYeG)(~j)oK;lY!)*!Gn1I3TtK7Iz|ztZ78Vu&0PS`gX__K2Gn??- eBqBY7$^8a>En|*mz$|?L0000 diff --git a/tests/_images/Graph_interaction_dendro.png b/tests/_images/Graph_interaction_dendro.png index 063ad6388b475d6b66a570287661e1269441e1e2..633f8c4495987c1bd8738b6bbb6b0ab9c3ad8d18 100644 GIT binary patch delta 1198 zcmV;f1X25f3giipMt?LeF)ScxbaZfYIxjD6VRUe8Z***FVlHoTXD`{*51;@5010qN zS#tmY1}6Xj1}6bcRM^J=00dDwK+v6NrL0!Cdb-roLCzZ>eq&}s!Um`o;7KBP&K1Q|!?U~+R6 zJe$oh8jV(V@PFxaiZo4EcDxv3{S8+ccW-YGo12?Ue>QgxPSX@$-TMQ#ufH26%7)#) z?*<Jo_p5m;i7mdl-DoYb`|%zyhRU(iBb3 z^na{MYEBhQZrQF%oc!2QG!mRDn48lonw%<_+^cT*%ldQbMRU_-5<7!>(acoA>|V9V zGVh~(ckViQ``56~lq@~^ZV+_5^QTW1-fLb5R|zh^xBer%6Epy6fu?A3N-#IID*s26 zU_z>3a!$Qya;9E1HFdJsdG(^H%~h~=ReuuF6wPg}7hUbDJZvWz2~UeGa<1|3g-S3Z zb+WiQuL|a-b1*mcqRBb+qRACH2UBy6ebFm>R%PgKzK$*Q?#I&O`(FedFa2%$@2k)2 zHyhW%RVMFb4*>vud>b?XSv~FOAg4|iJ5?|@XLW{^U`DE7aykcd)6E8Ue#v#P5PvbN zd(ql-N62dx&COfgi>@}^9Xx6`C<#ZaXmXmO$>}nQoLZISRKetQ4ko7rb2BBFTqD6% zFFG5a007SIC0_XF!H?~#3cDv~J`$2<{c+SqwkR%Cqb|-cK0A((-EJLr?yG(G)F-k=L0HAEP1{%)VrvLx| M07*qoM6N<$f(6_}X8-^I delta 1228 zcmV;-1T*{O34#icMt?IdGb|uzbaZfYIxjD6VRUe8Z***FVlHoTXD`T?x10a~010qN zS#tmY1}6Xj1}6bcRM^J=00eJIL_t(|obBDuOB`ny2k`H#qef+|OA(4m1C~Fn==6^>=M_6239J+pRaDY+Qe${Ffv$L}!eed-2G}hPGU#pL!Xti3{*x1Ou z-PP4qG@H$lzISP936qnPL)X2FtLSF4i63tNit95s^6Mpl?{0mYYkQgd4|bfi3jlX- z-f`0I!-;#jw;8q!8ecrd%F4>n*ZhlFyScduV6<$IH-AgJ-9G<3i|Jr6fVZ`^WwWgP zy|aTy``=&=^SO7x^T%^<698;<=A>QV?dfwD$C$%BcI^`Y$w<8M{Vg;;`#PC}ug~1T zNAD!*J%@Q*^FGM60{~DNzmj{K0AS*tlXii(r_V3#i=k`pUmHWC(a<)DRb&kS?SpB> zP=Kk$NPjj*3NWcC-TW%Rq~i3W0!%7~2S5R)71<=T3NWej3=YGv0XF%)2GMj6cQ3}XNkU|O-&Hi=p`L=PQ&_xy1u zXOT?L-PtKuDz!N$?V3NIIcfKz$R*uc!_FZ{U~4bj}JGMJl9w@E{%Ds6P;^zCs{40AES_aE%w)2S)!cDwiz z?|+usG0BHjp@UrVqB{kcl`@!}R?*~Ci{@rE2a~g^MN{)M2a{8P$r+uhS0C4x_9ZX*R1b_WMenaiYpK`sqlte=`H&2BQa?0TU2L+gv zhG=rCMYGc!%uTgucG(=PQBr2|UM<)5o)4m-YuPH=2Y`nY_poxsM}Ko7 z!sl@FOG&X*s${217B^)uJ7q98J#1j7eK0u%n45-ZcE%>PXq~Dgr3@x#Y;y2jrz%OA zq|%`bW~CvToHCf40!&UrG&yB3Ih888S!FP}G{8Opge?OgLfA6s;utUOi!!_R#<=Vi zfqA&>klcT9xc@h1uU22V)-}fn!$lBYuh+9d5C8z$?KXVhN3E9csXuY+h0;HMowYFqVA8b!0000wK+v6NrL0!Cdb-roLCzZ>eq&}s!Um`o;7KBP&K1Q|!?U~+R6 zJe$oh8jV(V@PFxaiZo4EcDxv3{S8+ccW-YGo12?Ue>QgxPSX@$-TMQ#ufH26%7)#) z?*<Jo_p5m;i7mdl-DoYb`|%zyhRU(iBb3 z^na{MYEBhQZrQF%oc!2QG!mRDn48lonw%<_+^cT*%ldQbMRU_-5<7!>(acoA>|V9V zGVh~(ckViQ``56~lq@~^ZV+_5^QTW1-fLb5R|zh^xBer%6Epy6fu?A3N-#IID*s26 zU_z>3a!$Qya;9E1HFdJsdG(^H%~h~=ReuuF6wPg}7hUbDJZvWz2~UeGa<1|3g-S3Z zb+WiQuL|a-b1*mcqRBb+qRACH2UBy6ebFm>R%PgKzK$*Q?#I&O`(FedFa2%$@2k)2 zHyhW%RVMFb4*>vud>b?XSv~FOAg4|iJ5?|@XLW{^U`DE7aykcd)6E8Ue#v#P5PvbN zd(ql-N62dx&COfgi>@}^9Xx6`C<#ZaXmXmO$>}nQoLZISRKetQ4ko7rb2BBFTqD6% zFFG5a007SIC0_XF!H?~#3cDv~J`$2<{c+SqwkR%Cqb|-cK0A((-EJLr?yG(G)F-k=L0HAEP1{%)VrvLx| M07*qoM6N<$f=WX~;{X5v delta 1225 zcmV;)1UCER3I7R@Mt?IdGb|uzbaZfYIxjD6VRUe8Z***FVlHoTXD`T?x10a~010qN zS#tmY1}6Xj1}6bcRM^J=00eAFL_t(|obBDuYaC@52k`H^8Io?o#v~Guv=$;Fx@4i` zUc-XG7rnK&~u6l7&(1+7-AsM{B!QmJ5TYinfpb43QY zTCHMdXQ!y!=jLTTs0RlJn46m$>U;C}Fhm#nX#oJVTCE@-)M~YQ;WkP+^g}K&1(@8q zA-dD)pwVay{eQ!2Yil?;Ir+~wKRi6d;^N}K^}W4447)U_R4SOBo*wFZCnqPdzP|oa zeH=!+-A27$FSOm|-hGg$GASVlwU6ae0}T7T-)orcYixby8v+K z=ItQuJ{-TBdzGgX1_ZxGNGuq73Ht%*GJMF-o;s zt?0f;0A7nMZV6P2rlwUiIR%(p0tJ|wR|a!?hA{vNFs;P2O`_Hxq6ZEgxl$>}StRR) zJ3AGME0=?`Ya{rDyPl;(bOfa5mdkIyf%fiSn16a4{N6GG8luT1D5bUNqB|mCsL(;q zTS>KOW_lh>P5~yTRWv!(qRAR z3V$8sviF@Tz@(JH?6is|r&=^O(;Q6BREwtOGzXJYfXR8Cs^q2&=H{}mI#tQatCGb{ zLo_#Slh~<}#myO1i>9U_nwvAI7EMhV%uTCkcG)U=WG^EP^r^Et)uP$;1FTb(oU~0Mry-i0%ig1;0JG8%&CPjbFgMM?+_cCdr-u#f zGzXK*9yW~ZWu&9{v*vlr`}AYq2On!v=QR2a{8PxoL=I=Pj*8>r^EvWiUB!>A`oMswCypN{2F- zm4;|?%3yK|FgXp;)<3+Cs%)^Tg z$z8<3?mw8mS~+*E?T--0G2H3tsX{dx4FEv1*+djYn4O)?wF45?@$oUL)hZSi761U9 nP6v%f1Kn;H&MaF`3#ESn=0LornZZFB00000NkvXXu0mjf3s+PN diff --git a/tests/conftest.py b/tests/conftest.py index f0d7ac0d..7697a27f 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -15,14 +15,12 @@ import pandas as pd from matplotlib.testing.compare import compare_images -import matplotlib as mpl import matplotlib.pyplot as plt from squidpy.im._container import ImageContainer from squidpy._constants._pkg_constants import Key import squidpy as sp -mpl.use("agg") HERE: Path = Path(__file__).parent EXPECTED = HERE / "_images" diff --git a/tests/graph/test_ppatterns.py b/tests/graph/test_ppatterns.py index c8cbf97a..2edb9113 100644 --- a/tests/graph/test_ppatterns.py +++ b/tests/graph/test_ppatterns.py @@ -5,9 +5,10 @@ from pandas.testing import assert_frame_equal import numpy as np -from squidpy.gr import moran, ripley_k, co_occurrence +from squidpy.gr import ripley_k, co_occurrence, spatial_autocorr MORAN_K = "moranI" +GEARY_C = "gearyC" def test_ripley_k(adata: AnnData): @@ -22,48 +23,67 @@ def test_ripley_k(adata: AnnData): assert cat_ripley.isdisjoint(cat_adata) is False -def test_moran_seq_par(dummy_adata: AnnData): - """Check whether moran results are the same for seq. and parallel computation.""" - moran(dummy_adata) +@pytest.mark.parametrize("mode", ["moran", "geary"]) +def test_spatial_autocorr_seq_par(dummy_adata: AnnData, mode: str): + """Check whether spatial autocorr results are the same for seq. and parallel computation.""" + spatial_autocorr(dummy_adata, mode=mode) dummy_adata.var["highly_variable"] = np.random.choice([True, False], size=dummy_adata.var_names.shape) - df = moran(dummy_adata, copy=True, n_jobs=1, seed=42, n_perms=50) - df_parallel = moran(dummy_adata, copy=True, n_jobs=2, seed=42, n_perms=50) + df = spatial_autocorr(dummy_adata, mode=mode, copy=True, n_jobs=1, seed=42, n_perms=50) + df_parallel = spatial_autocorr(dummy_adata, mode=mode, copy=True, n_jobs=2, seed=42, n_perms=50) idx_df = df.index.values idx_adata = dummy_adata[:, dummy_adata.var.highly_variable.values].var_names.values - assert MORAN_K in dummy_adata.uns.keys() - assert "pval_sim_fdr_bh" in dummy_adata.uns[MORAN_K] - assert dummy_adata.uns[MORAN_K].columns.shape == (4,) + if mode == "moran": + UNS_KEY = MORAN_K + elif mode == "geary": + UNS_KEY = GEARY_C + assert UNS_KEY in dummy_adata.uns.keys() + assert "pval_sim_fdr_bh" in df + assert "pval_norm_fdr_bh" in dummy_adata.uns[UNS_KEY] + assert dummy_adata.uns[UNS_KEY].columns.shape == (4,) + assert df.columns.shape == (9,) + # test pval_norm same + np.testing.assert_array_equal(df["pval_norm"].values, df["pval_norm"].values) # test highly variable - assert dummy_adata.uns[MORAN_K].shape != df.shape + assert dummy_adata.uns[UNS_KEY].shape != df.shape # assert idx are sorted and contain same elements assert not np.array_equal(idx_df, idx_adata) np.testing.assert_array_equal(sorted(idx_df), sorted(idx_adata)) # check parallel gives same results - with pytest.raises(AssertionError, match=r'.*\(column name="pval_sim"\) are different.*'): + with pytest.raises(AssertionError, match=r'.*\(column name="pval_z_sim"\) are different.*'): # because the seeds will be different, we don't expect the pval_sim values to be the same assert_frame_equal(df, df_parallel) +@pytest.mark.parametrize("mode", ["moran", "geary"]) @pytest.mark.parametrize("n_jobs", [1, 2]) -def test_moran_reproducibility(dummy_adata: AnnData, n_jobs: int): - """Check moran reproducibility results.""" - moran(dummy_adata) +def test_spatial_autocorr_reproducibility(dummy_adata: AnnData, n_jobs: int, mode: str): + """Check spatial autocorr reproducibility results.""" + spatial_autocorr(dummy_adata, mode=mode) dummy_adata.var["highly_variable"] = np.random.choice([True, False], size=dummy_adata.var_names.shape) # seed will work only when multiprocessing/loky - df_1 = moran(dummy_adata, copy=True, n_jobs=n_jobs, seed=42, n_perms=50) - df_2 = moran(dummy_adata, copy=True, n_jobs=n_jobs, seed=42, n_perms=50) + df_1 = spatial_autocorr(dummy_adata, mode=mode, copy=True, n_jobs=n_jobs, seed=42, n_perms=50) + df_2 = spatial_autocorr(dummy_adata, mode=mode, copy=True, n_jobs=n_jobs, seed=42, n_perms=50) idx_df = df_1.index.values idx_adata = dummy_adata[:, dummy_adata.var.highly_variable.values].var_names.values - assert MORAN_K in dummy_adata.uns.keys() + if mode == "moran": + UNS_KEY = MORAN_K + elif mode == "geary": + UNS_KEY = GEARY_C + assert UNS_KEY in dummy_adata.uns.keys() # assert fdr correction in adata.uns - assert "pval_sim_fdr_bh" in dummy_adata.uns[MORAN_K] - assert dummy_adata.uns[MORAN_K].columns.shape == (4,) + assert "pval_sim_fdr_bh" in df_1 + assert "pval_norm_fdr_bh" in dummy_adata.uns[UNS_KEY] + # test pval_norm same + np.testing.assert_array_equal(df_1["pval_norm"].values, df_2["pval_norm"].values) + np.testing.assert_array_equal(df_1["var_norm"].values, df_2["var_norm"].values) + assert dummy_adata.uns[UNS_KEY].columns.shape == (4,) + assert df_2.columns.shape == (9,) # test highly variable - assert dummy_adata.uns[MORAN_K].shape != df_1.shape + assert dummy_adata.uns[UNS_KEY].shape != df_1.shape # assert idx are sorted and contain same elements assert not np.array_equal(idx_df, idx_adata) np.testing.assert_array_equal(sorted(idx_df), sorted(idx_adata)) diff --git a/tests/plotting/test_interactive.py b/tests/plotting/test_interactive.py index f03dc26c..f8b56d68 100644 --- a/tests/plotting/test_interactive.py +++ b/tests/plotting/test_interactive.py @@ -11,7 +11,7 @@ from tests.conftest import DPI, PlotTester, PlotTesterMeta -@pytest.mark.qt +@pytest.mark.qt() class TestNapari(PlotTester, metaclass=PlotTesterMeta): def test_add_same_layer(self, qtbot, adata: AnnData, napari_cont: ImageContainer, capsys): from napari.layers import Points diff --git a/tox.ini b/tox.ini index fbeb9171..f29e095c 100644 --- a/tox.ini +++ b/tox.ini @@ -129,7 +129,7 @@ deps = # see: https://github.com/numba/llvmlite/issues/669 extras = all setenv = linux: PYTEST_FLAGS=--test-napari -passenv = TOXENV CI CODECOV_* GITHUB_ACTIONS PYTEST_FLAGS DISPLAY XAUTHORITY +passenv = TOXENV CI CODECOV_* GITHUB_ACTIONS PYTEST_FLAGS DISPLAY XAUTHORITY MPLBACKEND usedevelop = true commands = pytest --cov --cov-append --cov-report=term-missing --cov-config={toxinidir}/tox.ini --ignore docs/ {posargs:-vv} {env:PYTEST_FLAGS:} @@ -175,9 +175,11 @@ basepython = python3.8 deps = -r{toxinidir}/docs/requirements.txt skip_install = true whitelist_externals = sphinx-build +passenv = SQUIDPY_DOWNLOAD_NOTEBOOKS commands = - sphinx-build -W -b spelling {toxinidir}/docs/source {toxinidir}/docs/build/spellcheck - # sphinx-build -q -W --keep-going -b linkcheck {toxinidir}/docs/source {toxinidir}/docs/build/linkcheck + # do not exit on warning, in the CI, we disable the notebook download + sphinx-build -b spelling {toxinidir}/docs/source {toxinidir}/docs/build/spellcheck + sphinx-build -q --keep-going -b linkcheck {toxinidir}/docs/source {toxinidir}/docs/build/linkcheck [testenv:docs] description = Build the documentation.