diff --git a/docs/api/preprocessing.md b/docs/api/preprocessing.md index d820337987..a6acbfdcc7 100644 --- a/docs/api/preprocessing.md +++ b/docs/api/preprocessing.md @@ -29,6 +29,7 @@ For visual quality control, see {func}`~scanpy.pl.highest_expr_genes` and pp.log1p pp.pca pp.normalize_total + pp.normalize_clr pp.regress_out pp.scale pp.sample diff --git a/docs/references.bib b/docs/references.bib index aa2261fe7b..0b1453f989 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -111,6 +111,16 @@ @article{Blondel2008 pages = {P10008}, } +@article{Booeshaghi2022, + author = {Booeshaghi, A. Sina and Hallgrímsdóttir, Ingileif B. and Gálvez-Merchán, Ángel and Pachter, Lior}, + title = {Normalization for sampled count data}, + year = {2026}, + url = {https://doi.org/10.1101/2022.05.06.490859}, + doi = {10.1101/2022.05.06.490859}, + publisher = {Cold Spring Harbor Laboratory}, + journal = {bioRxiv}, +} + @article{Burczynski2006, author = {Burczynski, Michael E. and Peterson, Ron L. and Twine, Natalie C. and Zuberek, Krystyna A. and Brodeur, Brendan J. and Casciotti, Lori and Maganti, Vasu and Reddy, Padma S. and Strahs, Andrew and Immermann, Fred and Spinelli, Walter and Schwertschlag, Ulrich and Slager, Anna M. and Cotreau, Monette M. and Dorner, Andrew J.}, title = {Molecular Classification of Crohn’s Disease and Ulcerative Colitis Patients Using Transcriptional Profiles in Peripheral Blood Mononuclear Cells}, diff --git a/docs/release-notes/1.6.0.md b/docs/release-notes/1.6.0.md index 19b227fc05..6118ec3ee1 100644 --- a/docs/release-notes/1.6.0.md +++ b/docs/release-notes/1.6.0.md @@ -40,7 +40,7 @@ This release includes an overhaul of {func}`~scanpy.pl.dotplot`, {func}`~scanpy. #### Additions -- {func}`~anndata.concat` is now exported from scanpy, see {doc}`anndata:concatenation` for more info. {pr}`1338` {smaller}`I Virshup` +- {func}`~anndata.concat` is now exported from scanpy, see {doc}`anndata:tutorials/concatenation` for more info. {pr}`1338` {smaller}`I Virshup` - Added highly variable gene selection strategy from Seurat v3 {pr}`1204` {smaller}`A Gayoso` - Added [CellRank](https://github.com/theislab/cellrank/) to scanpy ecosystem {pr}`1304` {smaller}`giovp` - Added `backup_url` param to {func}`~scanpy.read_10x_h5` {pr}`1296` {smaller}`A Gayoso` diff --git a/docs/release-notes/4160.feat.md b/docs/release-notes/4160.feat.md new file mode 100644 index 0000000000..2c12d631db --- /dev/null +++ b/docs/release-notes/4160.feat.md @@ -0,0 +1 @@ +Add {func}`scanpy.pp.normalize_clr` for PFlog shifted centered log-ratio normalization, a variance-stabilizing, depth-invariant and rank-preserving count transform {cite:p}`Booeshaghi2022`. The default representation stores sparse shifted-log values together with per-cell centering offsets, and {func}`scanpy.pp.pca` can use this representation directly for memory-efficient PCA without materializing the dense centered matrix. {smaller}`R Baber` diff --git a/src/scanpy/preprocessing/__init__.py b/src/scanpy/preprocessing/__init__.py index cb1aedb90c..974e133124 100644 --- a/src/scanpy/preprocessing/__init__.py +++ b/src/scanpy/preprocessing/__init__.py @@ -7,7 +7,7 @@ from ._deprecated.sampling import subsample from ._harmony import harmony_integrate from ._highly_variable_genes import highly_variable_genes -from ._normalization import normalize_total +from ._normalization import normalize_clr, normalize_total from ._pca import pca from ._qc import calculate_qc_metrics from ._recipes import recipe_seurat, recipe_weinreb17, recipe_zheng17 @@ -33,6 +33,7 @@ "highly_variable_genes", "log1p", "neighbors", + "normalize_clr", "normalize_total", "pca", "recipe_seurat", diff --git a/src/scanpy/preprocessing/_normalization.py b/src/scanpy/preprocessing/_normalization.py index 79eb5cf0d1..a330203a67 100644 --- a/src/scanpy/preprocessing/_normalization.py +++ b/src/scanpy/preprocessing/_normalization.py @@ -7,6 +7,8 @@ import numpy as np from fast_array_utils import stats from fast_array_utils.numba import njit +from fast_array_utils.stats import mean_var +from scipy import sparse from .. import logging as logg from .._compat import CSBase, CSCBase, CSRBase, DaskArray, warn @@ -14,6 +16,8 @@ from ..get import _get_obs_rep, _set_obs_rep if TYPE_CHECKING: + from typing import Literal + from anndata import AnnData @@ -304,3 +308,281 @@ def normalize_total( # noqa: PLR0912 elif not inplace: return dat return None + + +def _estimate_overdispersion(x: np.ndarray | CSBase | DaskArray) -> float: + r"""Estimate the negative-binomial overdispersion :math:`α` from raw counts. + + Fits :math:`\mathrm{Var}_g = μ_g + α \cdot μ_g^2` across genes, where + :math:`μ_g` and :math:`\mathrm{Var}_g` are the per-gene mean and (population) + variance over cells. The model is linear in :math:`α`, so the ordinary + least-squares solution is closed form + + .. math:: + α = \frac{\sum_g (\mathrm{Var}_g - μ_g) \, μ_g^2}{\sum_g μ_g^4}, + + which is exactly the minimizer a non-linear `curve_fit` would converge to, + but without the dependency. :func:`~fast_array_utils.stats.mean_var` is + dispatched for dense, sparse and dask input alike, so the only dask-specific + step is computing the two final scalar sums. + """ + mu, var = mean_var(x, axis=0, correction=0) + mu2 = mu**2 + numerator = np.sum((var - mu) * mu2) + denominator = np.sum(mu2 * mu2) + if isinstance(x, DaskArray): + import dask + + numerator, denominator = dask.compute(numerator, denominator) + if denominator == 0.0: + msg = ( + "Cannot estimate overdispersion: every gene has zero mean. " + "Pass a positive `alpha` explicitly." + ) + raise ValueError(msg) + alpha = float(numerator / denominator) + if not alpha > 0: + msg = ( + f"Estimated overdispersion is non-positive (alpha = {alpha}); " + "pass a positive `alpha` explicitly." + ) + raise ValueError(msg) + return alpha + + +def _log1p_sparse_block(x: np.ndarray | CSBase) -> np.ndarray | CSBase: + """Apply log1p to a dense or sparse block while preserving sparse zeros.""" + if isinstance(x, CSBase): + x = x.copy() + x.data = np.log1p(x.data) + return x + return np.log1p(x) + + +def _densify_shifted_clr( + log_values: np.ndarray | CSBase | DaskArray, row_center: np.ndarray | DaskArray +) -> np.ndarray | DaskArray: + dense_log = log_values.toarray() if isinstance(log_values, CSBase) else log_values + return dense_log - row_center[:, None] + + +def _normalize_clr_helper( # noqa: PLR0912 + x: np.ndarray | CSBase | DaskArray, + *, + target: Literal["auto", "mean", "median"] | float, + alpha: float | None, +) -> tuple[ + np.ndarray | CSBase | DaskArray, + np.ndarray | DaskArray, + np.ndarray | DaskArray, + dict[str, object], +]: + """Compute sparse PFlog / shifted-CLR values and row centers.""" + # Keep the depths lazy for dask; `.ravel()` would otherwise materialize them. + cell_depths = stats.sum(x, axis=1) + if not isinstance(x, DaskArray): + cell_depths = np.asarray(cell_depths).ravel() + + mean_depth = ( + float(cell_depths.mean().compute()) + if isinstance(cell_depths, DaskArray) + else float(cell_depths.mean()) + ) + if alpha is not None: + if not alpha > 0: + msg = ( + f"`alpha` must be positive to compute PFlog, got {alpha}. " + "The data may be underdispersed." + ) + raise ValueError(msg) + scale: float | np.ndarray | DaskArray = 4.0 * float(alpha) + target_sum = 4.0 * float(alpha) * mean_depth + resolved_target = "alpha" + elif target == "auto": + alpha = _estimate_overdispersion(x) + scale = 4.0 * alpha + target_sum = 4.0 * alpha * mean_depth + resolved_target = "auto" + else: + alpha = None + if target == "mean": + target_sum = mean_depth + elif target == "median": + if isinstance(cell_depths, DaskArray): + msg = ( + "`target='median'` is not supported for dask input; pass " + "`target='auto'`, `target='mean'`, a positive numeric target, " + "or explicit `alpha`." + ) + raise NotImplementedError(msg) + target_sum = float(np.median(cell_depths)) + elif isinstance(target, bool) or not isinstance(target, int | float): + msg = "`target` must be 'auto', 'mean', 'median', or a positive number." + raise ValueError(msg) + else: + target_sum = float(target) + if not target_sum > 0: + msg = f"`target` must resolve to a positive depth, got {target_sum}." + raise ValueError(msg) + resolved_target = str(target) + x = axis_mul_or_truediv( + x, cell_depths / target_sum, op=truediv, axis=0, allow_divide_by_zero=False + ) + + if alpha is not None: + x = x * scale + + if isinstance(x, DaskArray): + log_values = x.map_blocks( + _log1p_sparse_block, dtype=np.float64, meta=x._meta.astype(np.float64) + ) + else: + log_values = _log1p_sparse_block(x) + + row_center = stats.sum(log_values, axis=1) / x.shape[1] + if not isinstance(row_center, DaskArray): + row_center = np.asarray(row_center).ravel() + report = dict( + target=resolved_target, + alpha=None if alpha is None else float(alpha), + k=float(target_sum), + mean_depth=mean_depth, + ) + return log_values, row_center, cell_depths, report + + +def normalize_clr( + adata: AnnData, + *, + target: Literal["auto", "mean", "median"] | float = "auto", + alpha: float | None = None, + key_added: str = "pflog", + densify: bool = False, + layer: str | None = None, + inplace: bool = True, + copy: bool = False, +) -> AnnData | dict[str, np.ndarray | CSBase | DaskArray | dict[str, object]] | None: + r"""Normalize counts with the shifted centered log-ratio (PFlog) transform. + + With `target="auto"` (or explicit `alpha`), computes PFlog as + + .. math:: + T(x)_i = \log(1 + 4 α x_i) + - \frac{1}{D} \sum_{j=1}^D \log(1 + 4 α x_j), + + which is equivalent to centering + :math:`\log(x_i + 1 / (4 α))` because the constant :math:`\log(4 α)` + cancels during CLR centering. For memory efficiency, the sparse log values + are stored separately from the per-cell centering vector. + + .. note:: + `target="auto"` is the recommended PFlog default. Fixed targets + (``"mean"``, ``"median"``, or a positive number) apply shifted CLR + after rescaling each cell to the same total count. + + Parameters + ---------- + adata + The annotated data matrix of shape `n_obs` × `n_vars`. + Rows correspond to cells and columns to genes. + target + ``"auto"`` estimates the negative-binomial overdispersion and applies + PFlog. ``"mean"``, ``"median"``, or a positive numeric value set a + fixed total count used before applying the shifted CLR transform. + alpha + Negative-binomial overdispersion of the dataset (``var = μ + α·μ²``). + When given, it overrides `target` and applies PFlog with sparse values + ``log1p(4 * alpha * x)``. + key_added + Layer key used to store the sparse log values. The row-centering vector + is stored in ``adata.obs[f"{key_added}_center"]``. + densify + If `True`, also materialize the dense centered matrix into `X` or the + selected `layer`. This is intended for small data or generic downstream + tools that cannot consume the sparse-plus-row-center representation. + layer + Layer to normalize instead of `X`. + inplace + Whether to update `adata` or return a dictionary with the normalized + matrix. + copy + Whether to modify a copied input object. Not compatible with + `inplace=False`. + + Returns + ------- + Returns a dictionary with normalized values and metadata or updates `adata`, + depending on `inplace`. + + Example + ------- + >>> import numpy as np + >>> from anndata import AnnData + >>> import scanpy as sc + >>> adata = AnnData(np.array([[1, 2, 30], [4, 50, 6]], dtype="float32")) + >>> sc.pp.normalize_clr(adata, alpha=0.5) + >>> "pflog" in adata.layers and "pflog_center" in adata.obs + True + """ + if copy: + if not inplace: + msg = "`copy=True` cannot be used with `inplace=False`." + raise ValueError(msg) + adata = adata.copy() + + view_to_actual(adata) + + x = _get_obs_rep(adata, layer=layer) + if isinstance(x, CSCBase): + x = x.tocsr() + if not inplace: + x = x.copy() + if issubclass(x.dtype.type, int | np.integer): + x = x.astype(np.float64) + + start = logg.info("normalizing counts per cell via PFlog") + + log_values, row_center, cell_depths, report = _normalize_clr_helper( + x, target=target, alpha=alpha + ) + + if not isinstance(cell_depths, DaskArray) and not np.all(cell_depths > 0): + warn("Some cells have zero counts", UserWarning) + + row_center_obs = ( + np.asarray(row_center.compute()).ravel() + if isinstance(row_center, DaskArray) + else np.asarray(row_center).ravel() + ) + if not isinstance(log_values, DaskArray | CSBase): + log_values = sparse.csr_matrix(log_values) # noqa: TID251 + + dense_x = _densify_shifted_clr(log_values, row_center) if densify else None + row_center_key = f"{key_added}_center" + metadata = dict( + encoding_type="shifted_clr", + row_center_key=row_center_key, + params=report | {"layer": layer, "densify": densify}, + ) + dat = dict( + X=dense_x if densify else log_values, + row_center=row_center_obs, + metadata=metadata, + ) + if inplace: + adata.layers[key_added] = log_values + adata.obs[row_center_key] = row_center_obs + adata.uns[key_added] = metadata + if densify: + _set_obs_rep(adata, dense_x, layer=layer) + + logg.info( + " finished ({time_passed})", + time=start, + ) + + if copy: + return adata + elif not inplace: + return dat + return None diff --git a/src/scanpy/preprocessing/_pca/__init__.py b/src/scanpy/preprocessing/_pca/__init__.py index 50fb6f4c1a..8f5a4b86e1 100644 --- a/src/scanpy/preprocessing/_pca/__init__.py +++ b/src/scanpy/preprocessing/_pca/__init__.py @@ -14,7 +14,7 @@ from ..._utils.random import _accepts_legacy_random_state, _legacy_random_state from ...get import _check_mask, _get_obs_rep from .._docs import doc_mask_var -from ._compat import _pca_compat_sparse +from ._compat import _pca_compat_sparse, _pca_compat_sparse_row_center if TYPE_CHECKING: from collections.abc import Container @@ -108,6 +108,10 @@ def pca( # noqa: PLR0912, PLR0913, PLR0915 or 1 - minimum dimension size of selected representation. layer If provided, which element of :attr:`~anndata.AnnData.layers` to use for PCA instead of `X`. + PFlog layers written by :func:`~scanpy.pp.normalize_clr` are represented + as sparse values plus a row-centering vector; PCA applies this centering + implicitly for in-memory sparse layers. Dask-backed PFlog layers are not + supported. obsm If provided, which element of :attr:`~anndata.AnnData.obsm` to use for PCA instead of `X`. zero_center @@ -239,6 +243,19 @@ def pca( # noqa: PLR0912, PLR0913, PLR0915 msg = f"PCA is not implemented for matrices of type {type(x)} from layers/obsm" raise NotImplementedError(msg) + row_center_key = None + if return_anndata and layer is not None and isinstance(adata.uns.get(layer), dict): + layer_metadata = adata.uns[layer] + if layer_metadata.get("encoding_type") == "shifted_clr": + row_center_key = layer_metadata.get("row_center_key") + if row_center_key is not None: + if row_center_key not in adata_comp.obs: + msg = f"{row_center_key!r} not found in `adata.obs`." + raise KeyError(msg) + row_center = np.asarray(adata_comp.obs[row_center_key], dtype=np.float64) + else: + row_center = None + if chunked: if not zero_center or svd_solver not in {None, "arpack"}: logg.debug("Ignoring zero_center, rng, svd_solver") @@ -267,7 +284,17 @@ def pca( # noqa: PLR0912, PLR0913, PLR0915 chunk_dense = chunk.toarray() if isinstance(chunk, CSBase) else chunk x_pca[start:end] = pca_.transform(chunk_dense) elif zero_center: - if isinstance(x, CSBase) and svd_solver == "lobpcg": + if row_center is not None: + if not isinstance(x, CSBase): + msg = "Sparse PFlog PCA is only supported for in-memory sparse input." + raise NotImplementedError(msg) + if svd_solver not in {None, "arpack"}: + msg = "Ignoring svd_solver for sparse PFlog PCA; using 'arpack'." + warn(msg, UserWarning) + x_pca, pca_ = _pca_compat_sparse_row_center( + x, row_center, n_comps, solver="arpack", rng=rng + ) + elif isinstance(x, CSBase) and svd_solver == "lobpcg": msg = ( f"{svd_solver=} for sparse relies on legacy code and will not be supported in the future. " "Also the lobpcg solver has been observed to be inaccurate. Please use 'arpack' instead." @@ -304,6 +331,9 @@ def pca( # noqa: PLR0912, PLR0913, PLR0915 ) x_pca = pca_.fit_transform(x) else: + if row_center is not None: + msg = "Sparse PFlog PCA requires `zero_center=True`." + raise NotImplementedError(msg) if isinstance(x, DaskArray): if isinstance(x._meta, CSBase): msg = ( diff --git a/src/scanpy/preprocessing/_pca/_compat.py b/src/scanpy/preprocessing/_pca/_compat.py index 064e791565..be3c98ef05 100644 --- a/src/scanpy/preprocessing/_pca/_compat.py +++ b/src/scanpy/preprocessing/_pca/_compat.py @@ -83,3 +83,84 @@ def rmat_op(v: NDArray[np.floating]): pca.explained_variance_ratio_ = ev_ratio pca.components_ = v return x_pca, pca + + +@_accepts_legacy_random_state(None) +def _pca_compat_sparse_row_center( + x: CSBase, + row_center: NDArray[np.floating], + n_pcs: int, + *, + solver: Literal["arpack", "lobpcg"], + rng: SeedLike | RNGLike | None = None, +) -> tuple[NDArray[np.floating], PCA]: + """Sparse PCA for matrices represented as ``x - row_center[:, None]``.""" + rng = np.random.default_rng(rng) + random_state_meta = _legacy_random_state(rng) + + x = check_array(x, accept_sparse=["csr", "csc"]) + row_center = np.asarray(row_center, dtype=np.float64).ravel() + if row_center.shape[0] != x.shape[0]: + msg = ( + f"`row_center` length must match the number of observations, " + f"got {row_center.shape[0]} for {x.shape[0]} observations." + ) + raise ValueError(msg) + + n_obs, _ = x.shape + row_center_sum = row_center.sum() + mu = (np.asarray(x.sum(axis=0)).ravel() - row_center_sum) / n_obs + mu = mu[None, :] + ones = np.ones(n_obs)[None, :].dot + + def mat_op(v: NDArray[np.floating]): + row_scale = np.sum(v, axis=0) + row_offset = ( + row_center * row_scale + if np.ndim(row_scale) == 0 + else row_center[:, None] * row_scale + ) + return (x @ v) - row_offset - (mu @ v) + + def rmat_op(v: NDArray[np.floating]): + row_offset = row_center @ v + if np.ndim(row_offset) > 0: + row_offset = row_offset[None, :] + return (x.T.conj() @ v) - row_offset - (mu.T @ ones(v)) + + linop = LinearOperator( + dtype=x.dtype, + shape=x.shape, + matvec=mat_op, + matmat=mat_op, + rmatvec=rmat_op, + rmatmat=rmat_op, + ) + + random_init = rng.uniform(size=np.min(x.shape)) + kw: Any = ( + dict(rng=rng) if SCIPY_1_15 else dict(random_state=_legacy_random_state(rng)) + ) + u, s, v = svds(linop, solver=solver, k=n_pcs, v0=random_init, **kw) + u, v = svd_flip(u, v, u_based_decision=not SCIPY_1_15) + idx = np.argsort(-s) + v = v[idx, :] + + x_pca = (u * s)[:, idx] + ev = s[idx] ** 2 / (n_obs - 1) + + x_sum = np.asarray(x.sum(axis=0)).ravel() + x_sq_sum = np.asarray(x.multiply(x).sum(axis=0)).ravel() + x_row_center_sum = np.asarray(x.T @ row_center).ravel() + y_sum = x_sum - row_center_sum + y_sq_sum = x_sq_sum - 2.0 * x_row_center_sum + np.sum(row_center**2) + total_var = np.sum((y_sq_sum - y_sum**2 / n_obs) / (n_obs - 1)) + ev_ratio = ev / total_var + + from sklearn.decomposition import PCA + + pca = PCA(n_components=n_pcs, svd_solver=solver, random_state=random_state_meta) + pca.explained_variance_ = ev + pca.explained_variance_ratio_ = ev_ratio + pca.components_ = v + return x_pca, pca diff --git a/tests/test_normalization.py b/tests/test_normalization.py index 50347c8ab9..a032f60b6a 100644 --- a/tests/test_normalization.py +++ b/tests/test_normalization.py @@ -11,15 +11,21 @@ from scipy import sparse import scanpy as sc +from scanpy._compat import CSBase from scanpy.preprocessing._normalization import _compute_nnz_median from testing.scanpy._helpers import ( _check_check_values_warnings, check_rep_mutation, check_rep_results, ) +from testing.scanpy._pytest.marks import needs # TODO: Add support for sparse-in-dask -from testing.scanpy._pytest.params import ARRAY_TYPES, ARRAY_TYPES_DENSE +from testing.scanpy._pytest.params import ( + ARRAY_TYPES, + ARRAY_TYPES_DENSE, + ARRAY_TYPES_MEM, +) if TYPE_CHECKING: from collections.abc import Callable @@ -329,3 +335,263 @@ def test_compute_nnz_median(array_type, dtype): data = np.array([0, 0, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9], dtype=dtype) data = array_type(data) np.testing.assert_allclose(_compute_nnz_median(data), 5) + + +# ------------------------------------------------------------------------------ +# normalize_clr (shifted CLR / PFlog) +# ------------------------------------------------------------------------------ + +# A small count matrix with no empty cells, used for the value/equivalence tests. +X_clr = np.array( + [[5, 0, 3, 2], [1, 1, 0, 4], [0, 7, 2, 1], [3, 3, 3, 3]], dtype="float32" +) + + +def _estimate_alpha_reference(x) -> float: + """Calculate OLS overdispersion for reference.""" + x = np.asarray(to_ndarray(x), dtype=np.float64) + mu = x.mean(axis=0) + var = (x**2).mean(axis=0) - mu**2 + mu2 = mu**2 + return float(np.sum((var - mu) * mu2) / np.sum(mu2 * mu2)) + + +def _clr_reference(x, *, target="auto", alpha=None) -> np.ndarray: + """Calculate shifted CLR densely for reference. + + PFlog uses a constant log1p(4 * alpha * x) scale. Depth targets use the + fixed-target scale log1p(K * x / depth). Empty cells are left as all-zero rows. + """ + x = np.asarray(to_ndarray(x), dtype=np.float64) + depths = x.sum(axis=1) + if alpha is not None: + log_u = np.log1p(4.0 * alpha * x) + elif target == "auto": + log_u = np.log1p(4.0 * _estimate_alpha_reference(x) * x) + else: + if target == "mean": + target_sum = depths.mean() + elif target == "median": + target_sum = np.median(depths) + else: + target_sum = float(target) + safe_depths = np.where(depths == 0, 1.0, depths) + log_u = np.log1p(x * (target_sum / safe_depths)[:, None]) + return log_u - log_u.mean(axis=1, keepdims=True) + + +def _reconstruct_clr(adata, key="pflog") -> np.ndarray: + return ( + to_ndarray(adata.layers[key]) - adata.obs[f"{key}_center"].to_numpy()[:, None] + ) + + +def _materialize(x): + if hasattr(x, "compute"): + x = x.compute() + return to_ndarray(x) + + +@pytest.mark.parametrize("array_type", ARRAY_TYPES_MEM) +@pytest.mark.parametrize("dtype", ["float32", "int64"]) +def test_normalize_clr_values(array_type, dtype): + """Check values against the reference and zero-sum cells.""" + adata = AnnData(array_type(X_clr).astype(dtype)) + sc.pp.normalize_clr(adata) + result = _reconstruct_clr(adata) + + np.testing.assert_allclose(result, _clr_reference(X_clr), rtol=1e-5, atol=1e-5) + # zero-sum (Aitchison) hyperplane + np.testing.assert_allclose(result.sum(axis=1), 0.0, atol=1e-5) + assert adata.layers["pflog"].nnz == sparse.csr_matrix(X_clr).nnz # noqa: TID251 + assert adata.uns["pflog"]["encoding_type"] == "shifted_clr" + assert adata.uns["pflog"]["row_center_key"] == "pflog_center" + assert adata.uns["pflog"]["params"]["target"] == "auto" + + +@pytest.mark.parametrize("array_type", ARRAY_TYPES_MEM) +@pytest.mark.parametrize( + "kwargs", + [{}, {"target": 1e4}, {"target": "mean"}, {"target": "median"}, {"alpha": 0.5}], + ids=["default", "fixed_target", "mean_target", "median_target", "alpha"], +) +def test_normalize_clr_params(array_type, kwargs): + adata = AnnData(array_type(X_clr).astype("float32")) + sc.pp.normalize_clr(adata, **kwargs) + np.testing.assert_allclose( + _reconstruct_clr(adata), _clr_reference(X_clr, **kwargs), rtol=1e-5, atol=1e-5 + ) + + +def test_normalize_clr_alpha_is_constant_scale_and_overrides_target(): + """`alpha` stores uncentered log1p(4 * alpha * x), independent of depth.""" + alpha = 0.5 + + via_alpha = AnnData(sparse.csr_matrix(X_clr)) # noqa: TID251 + sc.pp.normalize_clr(via_alpha, alpha=alpha) + np.testing.assert_allclose( + via_alpha.layers["pflog"].toarray(), + np.log1p(4.0 * alpha * X_clr), + rtol=1e-5, + atol=1e-5, + ) + + both = AnnData(sparse.csr_matrix(X_clr)) # noqa: TID251 + sc.pp.normalize_clr(both, alpha=alpha, target=999.0) + np.testing.assert_allclose( + _reconstruct_clr(both), _reconstruct_clr(via_alpha), rtol=1e-5, atol=1e-5 + ) + + +@pytest.mark.parametrize("array_type", ARRAY_TYPES_MEM) +def test_normalize_clr_alpha_auto(array_type): + """Check that `target="auto"` matches explicit alpha.""" + estimated = _estimate_alpha_reference(X_clr) + assert estimated > 0 + + auto = AnnData(array_type(X_clr).astype("float32")) + sc.pp.normalize_clr(auto) + + explicit = AnnData(array_type(X_clr).astype("float32")) + sc.pp.normalize_clr(explicit, alpha=estimated) + np.testing.assert_allclose( + _reconstruct_clr(auto), _reconstruct_clr(explicit), rtol=1e-5, atol=1e-5 + ) + + +@pytest.mark.parametrize("alpha", [0.0, -0.5], ids=["zero", "negative"]) +def test_normalize_clr_nonpositive_alpha_raises(alpha): + """Raise for non-positive `alpha`.""" + adata = AnnData(sparse.csr_matrix(X_clr)) # noqa: TID251 + with pytest.raises(ValueError, match=r"alpha.*positive"): + sc.pp.normalize_clr(adata, alpha=alpha) + + +def test_normalize_clr_alpha_auto_zero_mean_raises(): + """`target="auto"` cannot estimate overdispersion when every gene mean is zero.""" + adata = AnnData(np.zeros((3, 4), dtype="float32")) + with pytest.raises(ValueError, match="Cannot estimate overdispersion"): + sc.pp.normalize_clr(adata) + + +@pytest.mark.parametrize("array_type", ARRAY_TYPES_MEM) +def test_normalize_clr_zero_cell(array_type): + """Keep an empty cell finite and all-zero.""" + x = X_clr.copy() + x[1] = 0 # make the second cell empty + adata = AnnData(array_type(x)) + with pytest.warns(UserWarning, match="Some cells have zero counts"): + sc.pp.normalize_clr(adata) + result = _reconstruct_clr(adata) + assert np.isfinite(result).all() + np.testing.assert_allclose(result[1], 0.0, atol=1e-6) + + +def test_normalize_clr_inplace_false(): + adata = AnnData(sparse.csr_matrix(X_clr)) # noqa: TID251 + x_before = to_ndarray(adata.X).copy() + out = sc.pp.normalize_clr(adata, inplace=False) + + assert isinstance(out, dict) + np.testing.assert_allclose( + to_ndarray(out["X"]) - out["row_center"][:, None], + _clr_reference(X_clr), + rtol=1e-5, + atol=1e-5, + ) + # input is left untouched + assert isinstance(adata.X, CSBase) + np.testing.assert_array_equal(to_ndarray(adata.X), x_before) + + +def test_normalize_clr_copy(): + adata = AnnData(sparse.csr_matrix(X_clr)) # noqa: TID251 + returned = sc.pp.normalize_clr(adata, copy=True) + + assert isinstance(returned, AnnData) + assert returned is not adata + np.testing.assert_allclose( + _reconstruct_clr(returned), _clr_reference(X_clr), rtol=1e-5, atol=1e-5 + ) + # original is left untouched + assert isinstance(adata.X, CSBase) + + +def test_normalize_clr_copy_inplace_error(): + adata = AnnData(sparse.csr_matrix(X_clr)) # noqa: TID251 + with pytest.raises( + ValueError, match="`copy=True` cannot be used with `inplace=False`" + ): + sc.pp.normalize_clr(adata, copy=True, inplace=False) + + +def test_normalize_clr_layer(): + """`layer` selects the input layer and leaves `X` untouched.""" + adata = AnnData( + sparse.csr_matrix(X_clr), # noqa: TID251 + layers={"counts": sparse.csr_matrix(X_clr)}, # noqa: TID251 + ) + x_before = to_ndarray(adata.X).copy() + sc.pp.normalize_clr(adata, layer="counts") + + np.testing.assert_array_equal(to_ndarray(adata.X), x_before) + np.testing.assert_allclose( + _reconstruct_clr(adata), + _clr_reference(X_clr), + rtol=1e-5, + atol=1e-5, + ) + + +def test_normalize_clr_densify(): + adata = AnnData(sparse.csr_matrix(X_clr)) # noqa: TID251 + sc.pp.normalize_clr(adata, densify=True) + np.testing.assert_allclose(to_ndarray(adata.X), _clr_reference(X_clr), rtol=1e-5) + assert "pflog" in adata.layers + + +@needs.dask +@pytest.mark.parametrize("densify", [False, True], ids=["sparse_encoded", "densified"]) +@pytest.mark.parametrize( + "sparse_blocks", [False, True], ids=["dense_dask", "sparse_dask"] +) +def test_normalize_clr_dask(sparse_blocks, densify): + import dask.array as da + + chunks = (2, X_clr.shape[1]) + x = ( + da.from_array(sparse.csr_matrix(X_clr), chunks=chunks, asarray=False) # noqa: TID251 + if sparse_blocks + else da.from_array(X_clr, chunks=chunks) + ) + adata = AnnData(x.astype("float32")) + + sc.pp.normalize_clr(adata, densify=densify) + + result = ( + _materialize(adata.layers["pflog"]) + - adata.obs["pflog_center"].to_numpy()[:, None] + ) + np.testing.assert_allclose(result, _clr_reference(X_clr), rtol=1e-5, atol=1e-5) + if densify: + np.testing.assert_allclose( + _materialize(adata.X), _clr_reference(X_clr), rtol=1e-5, atol=1e-5 + ) + + +@needs.dask +def test_normalize_clr_dask_median_raises(): + import dask.array as da + + adata = AnnData(da.from_array(X_clr, chunks=(2, X_clr.shape[1]))) + + with pytest.raises(NotImplementedError, match=r"target='median'.*dask"): + sc.pp.normalize_clr(adata, target="median") + + +def test_normalize_clr_view(): + adata = AnnData(X_clr.copy()) + v = adata[:, :] + with pytest.warns(UserWarning, match=r"Received a view"): + sc.pp.normalize_clr(v) + assert not v.is_view diff --git a/tests/test_package_structure.py b/tests/test_package_structure.py index 9f08af6725..c077dd32a0 100644 --- a/tests/test_package_structure.py +++ b/tests/test_package_structure.py @@ -94,6 +94,9 @@ class ExpectedSig(TypedDict): copy_sigs["sc.pp.scale"]["first_name"] = "data" copy_sigs["sc.pp.sqrt"]["first_name"] = "data" # other partial exceptions +copy_sigs["sc.pp.normalize_clr"]["return_ann"] = ( + "AnnData | dict[str, np.ndarray | CSBase | DaskArray | dict[str, object]] | None" +) copy_sigs["sc.pp.normalize_total"]["return_ann"] = copy_sigs[ "sc.experimental.pp.normalize_pearson_residuals" ]["return_ann"] = "AnnData | dict[str, np.ndarray] | None" diff --git a/tests/test_pca.py b/tests/test_pca.py index 3611dabbf6..7c009ad207 100644 --- a/tests/test_pca.py +++ b/tests/test_pca.py @@ -11,6 +11,7 @@ from anndata.tests.helpers import assert_equal from packaging.version import Version from scipy import sparse +from sklearn.decomposition import PCA import scanpy as sc from scanpy._compat import CSBase, DaskArray, pkg_version @@ -331,6 +332,86 @@ def test_pca_sparse(key_added: str | None, keys: _PcaKeys): np.testing.assert_allclose(implicit.varm["PCs"], explicit.varm[keys.varm]) +def test_pca_sparse_shifted_clr_layer(): + counts = np.array( + [ + [20, 1, 0, 3, 0, 0], + [22, 2, 0, 4, 0, 0], + [1, 18, 0, 0, 3, 0], + [0, 20, 1, 0, 4, 0], + [0, 0, 15, 1, 0, 5], + [0, 0, 18, 2, 0, 4], + [12, 0, 0, 7, 1, 0], + [0, 12, 0, 1, 7, 0], + [0, 0, 12, 0, 1, 7], + ], + dtype=np.float64, + ) + adata = AnnData(sparse.csr_matrix(counts)) # noqa: TID251 + sc.pp.normalize_clr(adata, alpha=0.5) + sc.pp.pca(adata, n_comps=3, layer="pflog", dtype=np.float64) + + dense = ( + adata.layers["pflog"].toarray() - adata.obs["pflog_center"].to_numpy()[:, None] + ) + ref = PCA(n_components=3, svd_solver="full").fit(dense) + np.testing.assert_allclose(adata.uns["pca"]["variance"], ref.explained_variance_) + np.testing.assert_allclose( + adata.uns["pca"]["variance_ratio"], ref.explained_variance_ratio_ + ) + ref_scores = ref.transform(dense) + for i in range(3): + corr = abs(np.corrcoef(adata.obsm["X_pca"][:, i], ref_scores[:, i])[0, 1]) + assert corr > 0.999 + + +def test_pca_sparse_shifted_clr_layer_matches_dense_wide(): + rng = np.random.default_rng(0) + counts = rng.negative_binomial(4, 0.35, size=(15, 80)).astype(np.float64) + counts[rng.random(counts.shape) < 0.65] = 0 + counts[:, 0] += 1 # avoid all-zero rows + + adata = AnnData(sparse.csr_matrix(counts)) # noqa: TID251 + sc.pp.normalize_clr(adata, alpha=0.5) + sc.pp.pca(adata, n_comps=5, layer="pflog", dtype=np.float64, random_state=0) + + dense = ( + adata.layers["pflog"].toarray() - adata.obs["pflog_center"].to_numpy()[:, None] + ) + ref = PCA(n_components=5, svd_solver="full").fit(dense) + ref_scores = ref.transform(dense) + + scores = adata.obsm["X_pca"].copy() + loadings = adata.varm["PCs"].copy() + for i in range(5): + if np.corrcoef(scores[:, i], ref_scores[:, i])[0, 1] < 0: + scores[:, i] *= -1 + loadings[:, i] *= -1 + + np.testing.assert_allclose( + adata.uns["pca"]["variance"], ref.explained_variance_, rtol=1e-10 + ) + np.testing.assert_allclose( + adata.uns["pca"]["variance_ratio"], + ref.explained_variance_ratio_, + rtol=1e-10, + ) + np.testing.assert_allclose(scores, ref_scores, rtol=1e-7, atol=1e-7) + np.testing.assert_allclose(loadings, ref.components_.T, rtol=1e-10, atol=1e-10) + + +@needs.dask +def test_pca_dask_shifted_clr_layer_raises(): + import dask.array as da + + counts = da.from_array(A_list.astype(np.float64), chunks=(2, A_list.shape[1])) + adata = AnnData(counts) + sc.pp.normalize_clr(adata, alpha=0.5) + + with pytest.raises(NotImplementedError, match="in-memory sparse input"): + sc.pp.pca(adata, n_comps=2, layer="pflog") + + @pytest.mark.parametrize("rng_arg", ["rng", "random_state"]) def test_pca_reproducible( subtests: pytest.Subtests, array_type, rng_arg: Literal["rng", "random_state"]