From 6c253cb874d9109a697dfdef81eda902e0bd1dde Mon Sep 17 00:00:00 2001 From: rwbaber <97636743+rwbaber@users.noreply.github.com> Date: Sun, 14 Jun 2026 18:02:13 +0100 Subject: [PATCH 01/11] feat(preprocessing): add optimal shifted CLR (PFlog1pPF) normalization --- docs/api/preprocessing.md | 1 + docs/references.bib | 12 ++ docs/release-notes/+normalize-clr.feat.md | 1 + src/scanpy/preprocessing/__init__.py | 3 +- src/scanpy/preprocessing/_normalization.py | 203 +++++++++++++++++++++ tests/test_normalization.py | 164 ++++++++++++++++- 6 files changed, 382 insertions(+), 2 deletions(-) create mode 100644 docs/release-notes/+normalize-clr.feat.md 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..32f34308de 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -111,6 +111,18 @@ @article{Blondel2008 pages = {P10008}, } +@article{Booeshaghi2022, + author = {Booeshaghi, A. Sina and Hallgr{\'\i}msd{\'o}ttir, Ingileif B. and G{\'a}lvez-Merch{\'a}n, {\'A}ngel and Pachter, Lior}, + title = {Depth normalization for single-cell genomics count data}, + elocation-id = {2022.05.06.490859}, + year = {2026}, + doi = {10.1101/2022.05.06.490859}, + publisher = {Cold Spring Harbor Laboratory}, + url = {https://www.biorxiv.org/content/early/2026/06/09/2022.05.06.490859}, + eprint = {https://www.biorxiv.org/content/early/2026/06/09/2022.05.06.490859.full.pdf}, + 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/+normalize-clr.feat.md b/docs/release-notes/+normalize-clr.feat.md new file mode 100644 index 0000000000..02b93ce488 --- /dev/null +++ b/docs/release-notes/+normalize-clr.feat.md @@ -0,0 +1 @@ +Add {func}`scanpy.pp.normalize_clr` for shifted centered log-ratio (PFlog1pPF) normalization, a variance-stabilizing, depth-invariant and rank-preserving count transform {cite:p}`Booeshaghi2022` {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..fd21838a49 100644 --- a/src/scanpy/preprocessing/_normalization.py +++ b/src/scanpy/preprocessing/_normalization.py @@ -304,3 +304,206 @@ def normalize_total( # noqa: PLR0912 elif not inplace: return dat return None + + +# ----------------------------------------------------------------------------- +# Shifted Centered Log-Ratio (PFlogPF) Core Implementation +# +# The mathematical logic and sparse matrix optimization ("offset trick") +# implemented below are adapted from the Pachter Lab reference repository: +# https://github.com/pachterlab/bhgp_2022 +# +# Reference Publication: +# Booeshaghi et al. (2022, 2026) "Depth normalization for single-cell genomics +# count data" bioRxiv. doi: 10.1101/2022.05.06.490859 +# +# Copyright (c) 2022, Pachter Lab. All rights reserved. +# Licensed under the BSD 2-Clause License. +# ----------------------------------------------------------------------------- + + +def _normalize_clr_helper( + x: np.ndarray | CSBase, + *, + c: float, + target_sum: float | None, + alpha: float | None, + scale: float | None, +) -> tuple[np.ndarray, np.ndarray]: + """Compute the shifted CLR (PFlog1pPF) transform backend logic. + + This function acts as the internal processing engine for `normalize_clr`, + adapted from the reference implementations provided by the Pachter Lab + (Booeshaghi et al. 2022). It applies initial depth-scaling followed by + an optimized sparse log-shift operation before mapping coordinates onto + the zero-sum Aitchison hyperplane. + + Returns the dense, per-cell-centered matrix and the raw per-cell depths + (the latter is used by the caller to warn about empty cells). + + See `normalize_clr` for the meaning of the parameters. + """ + cell_depths = np.asarray(stats.sum(x, axis=1)).ravel() + + # Determine the target scale factor (sf, corresponding to K or s_f in the paper) + # for the initial proportional-fitting step. + if alpha is not None: + if scale is None: + scale = cell_depths.mean() + # Delta-method calibration: Sets the target scale factor (sf) to 4 * alpha * scale. + # This implicitly sets the count-scale pseudocount to the optimal + # variance-stabilizing value: y0 = 1 / (4 * alpha). + sf = 4.0 * alpha * scale + elif target_sum is not None: + sf = target_sum + else: + sf = cell_depths.mean() + + # Proportional fitting: rescale each cell so its counts sum to the target scale factor `sf`. + # Dividing by `cell_depths / sf` (with allow_divide_by_zero=False) leaves + # empty cells as all-zero rows instead of producing invalid infinities. + u = axis_mul_or_truediv( + x, cell_depths / sf, op=truediv, axis=0, allow_divide_by_zero=False + ) + + n_genes = x.shape[1] + log_c = np.log(c) + + if isinstance(u, CSBase): + # Sparse "offset trick": To avoid densifying the sparse matrix prematurely, + # we store log(u_nz + c) - log(c) at the nonzero entries only. + # Structural zeros keep their value of 0. The true value log(u + c) + # at any position can be recovered via: (stored value + log(c)). + # This defers matrix densification until the final centering step. + if c == 1.0: + # log(u + 1) - log(1) = log1p(u); log1p is more accurate near 0. + u.data = np.log1p(u.data) + else: + u.data = np.log(u.data + c) - log_c + # Per-cell mean over all genes, including the zeros: each of the + # (n_genes - nnz) zero positions contributes log(c), so the true row sum + # is (stored row sum) + n_genes * log(c). + # The true mean is therefore: row_sum / n_genes + log(c). + row_sum = np.asarray(u.sum(axis=1)).ravel() + per_cell_mean = row_sum / n_genes + log_c + # Materialize dense and add log(c) back so every position holds log(u + c). + dense_log = u.toarray() + log_c + else: + u = np.asarray(u) + dense_log = np.log1p(u) if c == 1.0 else np.log(u + c) + per_cell_mean = dense_log.mean(axis=1) + + # Additive centering (the CLR step): map each cell onto the zero-sum + # (Aitchison) hyperplane. This is what fills the zeros and densifies the matrix. + return dense_log - per_cell_mean[:, None], cell_depths + + +def normalize_clr( + adata: AnnData, + *, + c: float = 1.0, + target_sum: float | None = None, + alpha: float | None = None, + scale: float | None = None, + layer: str | None = None, + inplace: bool = True, + copy: bool = False, +) -> AnnData | dict[str, np.ndarray] | None: + r"""Normalize counts with the shifted centered log-ratio (PFlog1pPF) transform. + + Computes the shifted centered log-ratio (CLR) transform + + .. math:: + T(x)_i = \log(u_i + c) - \frac{1}{D} \sum_{j=1}^D \log(u_j + c), + + where :math:`u_i = s_f \, x_i / \sum_j x_j` are the depth-normalized counts + (proportional fitting to a target scale factor :math:`s_f`) and :math:`D` is the + number of genes. Equivalently this is proportional fitting, then + :func:`~numpy.log1p`, then per-cell mean-centering in log space (the + centered-log-ratio step). It is the unique count transform that is + simultaneously variance-stabilizing, depth-invariant, and rank-preserving + :cite:p:`Booeshaghi2022`. + + Because the per-cell centering fills the zero entries, the output is always + dense, and each cell sums to exactly zero. + + Parameters + ---------- + adata + The annotated data matrix of shape `n_obs` × `n_vars`. + Rows correspond to cells and columns to genes. + c + Pseudocount shift added inside the logarithm. The default ``1.0`` uses + :func:`~numpy.log1p` and keeps the sparse computation exact. + target_sum + Target scale factor :math:`s_f` for the first proportional-fitting step. If + `None` (and `alpha` is not given), the empirical mean cell depth is used. + alpha + Negative-binomial overdispersion of the dataset (``var = μ + α·μ²``). If + given, it overrides `target_sum` and sets :math:`s_f = 4 α \cdot` `scale` + by the delta method, calibrating the count-scale pseudocount to the + variance-stabilizing value ``y0 = 1/(4·α)``. + scale + Mean total counts per cell used in the delta-method rule. Only used when + `alpha` is given; defaults to the mean cell depth of the data. + 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 the normalized copy of `adata.X` (key `"X"`) or + updates `adata` with the normalized matrix, depending on `inplace`. + + Example + ------- + >>> import numpy as np + >>> from anndata import AnnData + >>> import scanpy as sc + >>> adata = AnnData(np.array([[1, 2, 3], [4, 5, 6]], dtype="float32")) + >>> sc.pp.normalize_clr(adata) + >>> np.allclose(adata.X.sum(axis=1), 0, atol=1e-6) + 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 issubclass(x.dtype.type, int | np.integer): + x = x.astype(np.float32) + + start = logg.info("normalizing counts per cell via shifted CLR") + + x, cell_depths = _normalize_clr_helper( + x, c=c, target_sum=target_sum, alpha=alpha, scale=scale + ) + + if not isinstance(cell_depths, DaskArray) and not np.all(cell_depths > 0): + warn("Some cells have zero counts", UserWarning) + + dat = dict(X=x) + if inplace: + _set_obs_rep(adata, dat["X"], layer=layer) + + logg.info( + " finished ({time_passed})", + time=start, + ) + + if copy: + return adata + elif not inplace: + return dat + return None diff --git a/tests/test_normalization.py b/tests/test_normalization.py index 50347c8ab9..815b975730 100644 --- a/tests/test_normalization.py +++ b/tests/test_normalization.py @@ -11,6 +11,7 @@ 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, @@ -19,7 +20,11 @@ ) # 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 +334,160 @@ 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 / PFlog1pPF) +# ------------------------------------------------------------------------------ + +# 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 _clr_reference(x, *, c=1.0, target_sum=None, alpha=None, scale=None) -> np.ndarray: + """Self-contained dense shifted-CLR, independent of the implementation. + + PF to a target depth, log(. + c), then subtract the per-cell mean. Empty + cells (zero depth) are left as all-zero rows, matching `normalize_clr`. + """ + x = np.asarray(to_ndarray(x), dtype=np.float64) + depths = x.sum(axis=1) + if alpha is not None: + if scale is None: + scale = depths.mean() + sf = 4.0 * alpha * scale + elif target_sum is not None: + sf = target_sum + else: + sf = depths.mean() + safe_depths = np.where(depths == 0, 1.0, depths) + u = x * (sf / safe_depths)[:, None] + log_u = np.log(u + c) + return log_u - log_u.mean(axis=1, keepdims=True) + + +@pytest.mark.parametrize("array_type", ARRAY_TYPES_MEM) +@pytest.mark.parametrize("dtype", ["float32", "int64"]) +def test_normalize_clr_values(array_type, dtype): + """Values match the reference and every cell sums to zero, for all layouts. + + Asserting that numpy / csr / csc / sparse-array inputs all equal the same + dense reference also proves the sparse "offset trick" matches the dense path. + """ + adata = AnnData(array_type(X_clr).astype(dtype)) + sc.pp.normalize_clr(adata) + result = to_ndarray(adata.X) + + 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) + + +@pytest.mark.parametrize("array_type", ARRAY_TYPES_MEM) +@pytest.mark.parametrize( + "kwargs", + [{}, {"target_sum": 1e4}, {"c": 2.0}, {"alpha": 0.5}], + ids=["default", "target_sum", "shift_c", "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( + to_ndarray(adata.X), _clr_reference(X_clr, **kwargs), rtol=1e-5, atol=1e-5 + ) + + +def test_normalize_clr_alpha_overrides_target_sum(): + """`alpha` sets sf = 4*alpha*scale and overrides any given `target_sum`.""" + alpha = 0.5 + scale = X_clr.sum(axis=1).mean() + + via_alpha = AnnData(sparse.csr_matrix(X_clr)) # noqa: TID251 + sc.pp.normalize_clr(via_alpha, alpha=alpha) + + via_target = AnnData(sparse.csr_matrix(X_clr)) # noqa: TID251 + sc.pp.normalize_clr(via_target, target_sum=4.0 * alpha * scale) + np.testing.assert_allclose( + to_ndarray(via_alpha.X), to_ndarray(via_target.X), rtol=1e-5, atol=1e-5 + ) + + # passing both -> alpha wins, target_sum ignored + both = AnnData(sparse.csr_matrix(X_clr)) # noqa: TID251 + sc.pp.normalize_clr(both, alpha=alpha, target_sum=999.0) + np.testing.assert_allclose( + to_ndarray(both.X), to_ndarray(via_alpha.X), rtol=1e-5, atol=1e-5 + ) + + +@pytest.mark.parametrize("array_type", ARRAY_TYPES_MEM) +def test_normalize_clr_zero_cell(array_type): + """An empty cell stays all-zero, stays finite, and triggers a warning.""" + 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 = to_ndarray(adata.X) + 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(out["X"], _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( + to_ndarray(returned.X), _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` targets that 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( + to_ndarray(adata.layers["counts"]), + _clr_reference(X_clr), + rtol=1e-5, + atol=1e-5, + ) + + +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 From eec2659114f4c495be364e9c50467ebdbe57c357 Mon Sep 17 00:00:00 2001 From: rwbaber <97636743+rwbaber@users.noreply.github.com> Date: Sun, 14 Jun 2026 18:31:41 +0100 Subject: [PATCH 02/11] docs: update release notes fragment with PR number --- docs/release-notes/{+normalize-clr.feat.md => 4160.feat.md} | 0 src/scanpy/preprocessing/_normalization.py | 5 ++--- 2 files changed, 2 insertions(+), 3 deletions(-) rename docs/release-notes/{+normalize-clr.feat.md => 4160.feat.md} (100%) diff --git a/docs/release-notes/+normalize-clr.feat.md b/docs/release-notes/4160.feat.md similarity index 100% rename from docs/release-notes/+normalize-clr.feat.md rename to docs/release-notes/4160.feat.md diff --git a/src/scanpy/preprocessing/_normalization.py b/src/scanpy/preprocessing/_normalization.py index fd21838a49..c683c49a11 100644 --- a/src/scanpy/preprocessing/_normalization.py +++ b/src/scanpy/preprocessing/_normalization.py @@ -420,9 +420,8 @@ def normalize_clr( (proportional fitting to a target scale factor :math:`s_f`) and :math:`D` is the number of genes. Equivalently this is proportional fitting, then :func:`~numpy.log1p`, then per-cell mean-centering in log space (the - centered-log-ratio step). It is the unique count transform that is - simultaneously variance-stabilizing, depth-invariant, and rank-preserving - :cite:p:`Booeshaghi2022`. + centered-log-ratio step). This count transform is simultaneously + variance-stabilizing, depth-invariant, and rank-preserving :cite:p:`Booeshaghi2022`. Because the per-cell centering fills the zero entries, the output is always dense, and each cell sums to exactly zero. From 32fb68be003d2ea32c3e8a145f2817f201494311 Mon Sep 17 00:00:00 2001 From: rwbaber <97636743+rwbaber@users.noreply.github.com> Date: Sun, 14 Jun 2026 19:21:59 +0100 Subject: [PATCH 03/11] test: register normalize_clr in copy_sigs signature conventions --- src/scanpy/preprocessing/_normalization.py | 4 ++-- tests/test_package_structure.py | 4 +++- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/scanpy/preprocessing/_normalization.py b/src/scanpy/preprocessing/_normalization.py index c683c49a11..30afb365e9 100644 --- a/src/scanpy/preprocessing/_normalization.py +++ b/src/scanpy/preprocessing/_normalization.py @@ -419,7 +419,7 @@ def normalize_clr( where :math:`u_i = s_f \, x_i / \sum_j x_j` are the depth-normalized counts (proportional fitting to a target scale factor :math:`s_f`) and :math:`D` is the number of genes. Equivalently this is proportional fitting, then - :func:`~numpy.log1p`, then per-cell mean-centering in log space (the + :obj:`~numpy.log1p`, then per-cell mean-centering in log space (the centered-log-ratio step). This count transform is simultaneously variance-stabilizing, depth-invariant, and rank-preserving :cite:p:`Booeshaghi2022`. @@ -433,7 +433,7 @@ def normalize_clr( Rows correspond to cells and columns to genes. c Pseudocount shift added inside the logarithm. The default ``1.0`` uses - :func:`~numpy.log1p` and keeps the sparse computation exact. + :obj:`~numpy.log1p` and keeps the sparse computation exact. target_sum Target scale factor :math:`s_f` for the first proportional-fitting step. If `None` (and `alpha` is not given), the empirical mean cell depth is used. diff --git a/tests/test_package_structure.py b/tests/test_package_structure.py index 9f08af6725..95789eca00 100644 --- a/tests/test_package_structure.py +++ b/tests/test_package_structure.py @@ -96,7 +96,9 @@ class ExpectedSig(TypedDict): # other partial exceptions copy_sigs["sc.pp.normalize_total"]["return_ann"] = copy_sigs[ "sc.experimental.pp.normalize_pearson_residuals" -]["return_ann"] = "AnnData | dict[str, np.ndarray] | None" +]["return_ann"] = copy_sigs["sc.pp.normalize_clr"]["return_ann"] = ( + "AnnData | dict[str, np.ndarray] | None" +) copy_sigs["sc.external.pp.magic"]["copy_default"] = None From 6bd44b0eec0bce76ef52e14ccea587df2548e808 Mon Sep 17 00:00:00 2001 From: rwbaber <97636743+rwbaber@users.noreply.github.com> Date: Sun, 14 Jun 2026 21:22:12 +0100 Subject: [PATCH 04/11] refactor(pp): unify normalize_clr math variables and clean up bibliography --- docs/references.bib | 16 +++---- src/scanpy/preprocessing/_normalization.py | 53 ++++++++++++---------- tests/test_normalization.py | 13 +++--- 3 files changed, 42 insertions(+), 40 deletions(-) diff --git a/docs/references.bib b/docs/references.bib index 32f34308de..906747e0c2 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -112,15 +112,13 @@ @article{Blondel2008 } @article{Booeshaghi2022, - author = {Booeshaghi, A. Sina and Hallgr{\'\i}msd{\'o}ttir, Ingileif B. and G{\'a}lvez-Merch{\'a}n, {\'A}ngel and Pachter, Lior}, - title = {Depth normalization for single-cell genomics count data}, - elocation-id = {2022.05.06.490859}, - year = {2026}, - doi = {10.1101/2022.05.06.490859}, - publisher = {Cold Spring Harbor Laboratory}, - url = {https://www.biorxiv.org/content/early/2026/06/09/2022.05.06.490859}, - eprint = {https://www.biorxiv.org/content/early/2026/06/09/2022.05.06.490859.full.pdf}, - journal = {bioRxiv}, + author = {Booeshaghi, A. Sina and Hallgrímsdóttir, Ingileif B. and Gálvez-Merchán, Ángel and Pachter, Lior}, + title = {Depth normalization for single-cell genomics 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, diff --git a/src/scanpy/preprocessing/_normalization.py b/src/scanpy/preprocessing/_normalization.py index 30afb365e9..c0ce810ce4 100644 --- a/src/scanpy/preprocessing/_normalization.py +++ b/src/scanpy/preprocessing/_normalization.py @@ -345,25 +345,26 @@ def _normalize_clr_helper( """ cell_depths = np.asarray(stats.sum(x, axis=1)).ravel() - # Determine the target scale factor (sf, corresponding to K or s_f in the paper) - # for the initial proportional-fitting step. + # Resolve `target_sum`: the depth that the first proportional-fitting step + # scales every cell to. This is the paper's PF target constant K (Booeshaghi + # et al. 2022); the three branches below are the three ways to choose it. if alpha is not None: if scale is None: scale = cell_depths.mean() - # Delta-method calibration: Sets the target scale factor (sf) to 4 * alpha * scale. - # This implicitly sets the count-scale pseudocount to the optimal - # variance-stabilizing value: y0 = 1 / (4 * alpha). - sf = 4.0 * alpha * scale - elif target_sum is not None: - sf = target_sum - else: - sf = cell_depths.mean() - - # Proportional fitting: rescale each cell so its counts sum to the target scale factor `sf`. - # Dividing by `cell_depths / sf` (with allow_divide_by_zero=False) leaves - # empty cells as all-zero rows instead of producing invalid infinities. + # Delta-method calibration: K = 4 * alpha * s, where s (`scale`) is the + # mean cell depth. This implicitly sets the count-scale pseudocount to the + # optimal variance-stabilizing value y0 = 1 / (4 * alpha), and overrides + # any value passed as `target_sum`. + target_sum = 4.0 * alpha * scale + elif target_sum is None: + target_sum = cell_depths.mean() + # else: use the `target_sum` passed by the caller. + + # Proportional fitting: rescale each cell so its counts sum to `target_sum` (K). + # Dividing by `cell_depths / target_sum` (with allow_divide_by_zero=False) + # leaves empty cells as all-zero rows instead of producing invalid infinities. u = axis_mul_or_truediv( - x, cell_depths / sf, op=truediv, axis=0, allow_divide_by_zero=False + x, cell_depths / target_sum, op=truediv, axis=0, allow_divide_by_zero=False ) n_genes = x.shape[1] @@ -416,8 +417,8 @@ def normalize_clr( .. math:: T(x)_i = \log(u_i + c) - \frac{1}{D} \sum_{j=1}^D \log(u_j + c), - where :math:`u_i = s_f \, x_i / \sum_j x_j` are the depth-normalized counts - (proportional fitting to a target scale factor :math:`s_f`) and :math:`D` is the + where :math:`u_i = K \, x_i / \sum_j x_j` are the depth-normalized counts + (proportional fitting to a target depth :math:`K`) and :math:`D` is the number of genes. Equivalently this is proportional fitting, then :obj:`~numpy.log1p`, then per-cell mean-centering in log space (the centered-log-ratio step). This count transform is simultaneously @@ -435,16 +436,20 @@ def normalize_clr( Pseudocount shift added inside the logarithm. The default ``1.0`` uses :obj:`~numpy.log1p` and keeps the sparse computation exact. target_sum - Target scale factor :math:`s_f` for the first proportional-fitting step. If - `None` (and `alpha` is not given), the empirical mean cell depth is used. + Target depth :math:`K` for the first proportional-fitting step. If `None` + (and `alpha` is not given), the empirical mean cell depth is used. Unlike + :func:`~scanpy.pp.normalize_total`, this is only an *intermediate* target: + cells sum to `target_sum` after proportional fitting, but the subsequent + log and centering steps make each cell sum to exactly zero in the output. alpha Negative-binomial overdispersion of the dataset (``var = μ + α·μ²``). If - given, it overrides `target_sum` and sets :math:`s_f = 4 α \cdot` `scale` - by the delta method, calibrating the count-scale pseudocount to the - variance-stabilizing value ``y0 = 1/(4·α)``. + given, it overrides `target_sum` and sets :math:`K = 4 \cdot α \cdot s` + by the delta method, where :math:`s` is `scale`. This calibrates the + count-scale pseudocount to the variance-stabilizing value ``y0 = 1/(4·α)``. scale - Mean total counts per cell used in the delta-method rule. Only used when - `alpha` is given; defaults to the mean cell depth of the data. + Mean total counts per cell (:math:`s` in :cite:p:`Booeshaghi2022`), used + in the delta-method rule :math:`K = 4 α s`. Only used when `alpha` is + given; defaults to the mean cell depth of the data. layer Layer to normalize instead of `X`. inplace diff --git a/tests/test_normalization.py b/tests/test_normalization.py index 815b975730..28884bb2c0 100644 --- a/tests/test_normalization.py +++ b/tests/test_normalization.py @@ -354,16 +354,15 @@ def _clr_reference(x, *, c=1.0, target_sum=None, alpha=None, scale=None) -> np.n """ x = np.asarray(to_ndarray(x), dtype=np.float64) depths = x.sum(axis=1) + # Resolve the PF target depth K (paper's notation; exposed as `target_sum`). if alpha is not None: if scale is None: scale = depths.mean() - sf = 4.0 * alpha * scale - elif target_sum is not None: - sf = target_sum - else: - sf = depths.mean() + target_sum = 4.0 * alpha * scale + elif target_sum is None: + target_sum = depths.mean() safe_depths = np.where(depths == 0, 1.0, depths) - u = x * (sf / safe_depths)[:, None] + u = x * (target_sum / safe_depths)[:, None] log_u = np.log(u + c) return log_u - log_u.mean(axis=1, keepdims=True) @@ -400,7 +399,7 @@ def test_normalize_clr_params(array_type, kwargs): def test_normalize_clr_alpha_overrides_target_sum(): - """`alpha` sets sf = 4*alpha*scale and overrides any given `target_sum`.""" + """`alpha` sets target_sum = 4*alpha*scale and overrides any given `target_sum`.""" alpha = 0.5 scale = X_clr.sum(axis=1).mean() From dd9aca0fd4b2589ebc4906ae20c42e90d20d256a Mon Sep 17 00:00:00 2001 From: rwbaber <97636743+rwbaber@users.noreply.github.com> Date: Sun, 14 Jun 2026 22:11:14 +0100 Subject: [PATCH 05/11] docs: polish normalize_clr docstring wording --- src/scanpy/preprocessing/_normalization.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/scanpy/preprocessing/_normalization.py b/src/scanpy/preprocessing/_normalization.py index c0ce810ce4..e5d7913467 100644 --- a/src/scanpy/preprocessing/_normalization.py +++ b/src/scanpy/preprocessing/_normalization.py @@ -447,9 +447,9 @@ def normalize_clr( by the delta method, where :math:`s` is `scale`. This calibrates the count-scale pseudocount to the variance-stabilizing value ``y0 = 1/(4·α)``. scale - Mean total counts per cell (:math:`s` in :cite:p:`Booeshaghi2022`), used - in the delta-method rule :math:`K = 4 α s`. Only used when `alpha` is - given; defaults to the mean cell depth of the data. + Mean total counts per cell, used as :math:`s` in the delta-method + rule :math:`K = 4 α s`. Only used when `alpha` is given; defaults + to the mean cell depth of the data. layer Layer to normalize instead of `X`. inplace From f652a7e0cb0341e9ba048eaa0dd295957bc94dcb Mon Sep 17 00:00:00 2001 From: rwbaber <97636743+rwbaber@users.noreply.github.com> Date: Mon, 15 Jun 2026 15:41:43 +0100 Subject: [PATCH 06/11] refactor(pp): finalize normalize_clr API per review Address maintainer review on #4160: - Drop c (redundant with target_sum) and scale (internal mean depth). - alpha now accepts 'auto' via closed-form OLS overdispersion estimation matching runorm. - Raise ValueError on non-positive alpha instead of clamping. - Cast integer input to float64 instead of float32. --- src/scanpy/preprocessing/_normalization.py | 127 ++++++++++++--------- tests/test_normalization.py | 51 +++++++-- 2 files changed, 117 insertions(+), 61 deletions(-) diff --git a/src/scanpy/preprocessing/_normalization.py b/src/scanpy/preprocessing/_normalization.py index e5d7913467..bda84c57fb 100644 --- a/src/scanpy/preprocessing/_normalization.py +++ b/src/scanpy/preprocessing/_normalization.py @@ -14,6 +14,8 @@ from ..get import _get_obs_rep, _set_obs_rep if TYPE_CHECKING: + from typing import Literal + from anndata import AnnData @@ -322,13 +324,47 @@ def normalize_total( # noqa: PLR0912 # ----------------------------------------------------------------------------- +def _estimate_overdispersion(x: np.ndarray | CSBase) -> 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 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. + """ + n_cells = x.shape[0] + col_sum = np.asarray(stats.sum(x, axis=0)).ravel() + if isinstance(x, CSBase): + x_sq = x.copy() + x_sq.data **= 2 + else: + x_sq = np.square(x) + col_sum_sq = np.asarray(stats.sum(x_sq, axis=0)).ravel() + + mu = col_sum / n_cells + var = col_sum_sq / n_cells - mu**2 + mu2 = mu**2 + denominator = float(np.sum(mu2 * mu2)) + if denominator == 0.0: + msg = ( + "Cannot estimate overdispersion: every gene has zero mean. " + "Pass `alpha` or `target_sum` explicitly." + ) + raise ValueError(msg) + return float(np.sum((var - mu) * mu2) / denominator) + + def _normalize_clr_helper( x: np.ndarray | CSBase, *, - c: float, target_sum: float | None, - alpha: float | None, - scale: float | None, + alpha: float | Literal["auto"] | None, ) -> tuple[np.ndarray, np.ndarray]: """Compute the shifted CLR (PFlog1pPF) transform backend logic. @@ -345,17 +381,23 @@ def _normalize_clr_helper( """ cell_depths = np.asarray(stats.sum(x, axis=1)).ravel() - # Resolve `target_sum`: the depth that the first proportional-fitting step - # scales every cell to. This is the paper's PF target constant K (Booeshaghi - # et al. 2022); the three branches below are the three ways to choose it. + # Resolve `target_sum` (the paper's PF target constant K, Booeshaghi et al. + # 2022): the depth the first proportional-fitting step scales every cell to. if alpha is not None: - if scale is None: - scale = cell_depths.mean() - # Delta-method calibration: K = 4 * alpha * s, where s (`scale`) is the - # mean cell depth. This implicitly sets the count-scale pseudocount to the - # optimal variance-stabilizing value y0 = 1 / (4 * alpha), and overrides - # any value passed as `target_sum`. - target_sum = 4.0 * alpha * scale + if alpha == "auto": + # Estimate the overdispersion from the raw counts (closed-form OLS, + # as in the reference implementation). + alpha = _estimate_overdispersion(x) + if not alpha > 0: + msg = ( + f"`alpha` must be positive to derive K = 4 * alpha * s, got {alpha}. " + "The data may be underdispersed; pass `target_sum` explicitly instead." + ) + raise ValueError(msg) + # Delta-method calibration: K = 4 * alpha * s, where s is the mean cell + # depth. This sets the count-scale pseudocount to the variance-stabilizing + # value y0 = 1 / (4 * alpha), and overrides any value passed as `target_sum`. + target_sum = 4.0 * alpha * cell_depths.mean() elif target_sum is None: target_sum = cell_depths.mean() # else: use the `target_sum` passed by the caller. @@ -368,30 +410,18 @@ def _normalize_clr_helper( ) n_genes = x.shape[1] - log_c = np.log(c) if isinstance(u, CSBase): - # Sparse "offset trick": To avoid densifying the sparse matrix prematurely, - # we store log(u_nz + c) - log(c) at the nonzero entries only. - # Structural zeros keep their value of 0. The true value log(u + c) - # at any position can be recovered via: (stored value + log(c)). - # This defers matrix densification until the final centering step. - if c == 1.0: - # log(u + 1) - log(1) = log1p(u); log1p is more accurate near 0. - u.data = np.log1p(u.data) - else: - u.data = np.log(u.data + c) - log_c - # Per-cell mean over all genes, including the zeros: each of the - # (n_genes - nnz) zero positions contributes log(c), so the true row sum - # is (stored row sum) + n_genes * log(c). - # The true mean is therefore: row_sum / n_genes + log(c). - row_sum = np.asarray(u.sum(axis=1)).ravel() - per_cell_mean = row_sum / n_genes + log_c - # Materialize dense and add log(c) back so every position holds log(u + c). - dense_log = u.toarray() + log_c + # Sparse "offset trick": store log1p only at the nonzero entries. The + # structural zeros already hold log1p(0) = 0, so the per-cell mean is just + # the stored row sum / n_genes, and no densification is needed until the + # final centering step. + u.data = np.log1p(u.data) + per_cell_mean = np.asarray(u.sum(axis=1)).ravel() / n_genes + dense_log = u.toarray() else: u = np.asarray(u) - dense_log = np.log1p(u) if c == 1.0 else np.log(u + c) + dense_log = np.log1p(u) per_cell_mean = dense_log.mean(axis=1) # Additive centering (the CLR step): map each cell onto the zero-sum @@ -402,10 +432,8 @@ def _normalize_clr_helper( def normalize_clr( adata: AnnData, *, - c: float = 1.0, target_sum: float | None = None, - alpha: float | None = None, - scale: float | None = None, + alpha: float | Literal["auto"] | None = None, layer: str | None = None, inplace: bool = True, copy: bool = False, @@ -415,7 +443,7 @@ def normalize_clr( Computes the shifted centered log-ratio (CLR) transform .. math:: - T(x)_i = \log(u_i + c) - \frac{1}{D} \sum_{j=1}^D \log(u_j + c), + T(x)_i = \log(u_i + 1) - \frac{1}{D} \sum_{j=1}^D \log(u_j + 1), where :math:`u_i = K \, x_i / \sum_j x_j` are the depth-normalized counts (proportional fitting to a target depth :math:`K`) and :math:`D` is the @@ -432,9 +460,6 @@ def normalize_clr( adata The annotated data matrix of shape `n_obs` × `n_vars`. Rows correspond to cells and columns to genes. - c - Pseudocount shift added inside the logarithm. The default ``1.0`` uses - :obj:`~numpy.log1p` and keeps the sparse computation exact. target_sum Target depth :math:`K` for the first proportional-fitting step. If `None` (and `alpha` is not given), the empirical mean cell depth is used. Unlike @@ -442,14 +467,14 @@ def normalize_clr( cells sum to `target_sum` after proportional fitting, but the subsequent log and centering steps make each cell sum to exactly zero in the output. alpha - Negative-binomial overdispersion of the dataset (``var = μ + α·μ²``). If - given, it overrides `target_sum` and sets :math:`K = 4 \cdot α \cdot s` - by the delta method, where :math:`s` is `scale`. This calibrates the - count-scale pseudocount to the variance-stabilizing value ``y0 = 1/(4·α)``. - scale - Mean total counts per cell, used as :math:`s` in the delta-method - rule :math:`K = 4 α s`. Only used when `alpha` is given; defaults - to the mean cell depth of the data. + Negative-binomial overdispersion of the dataset (``var = μ + α·μ²``). When + given, it overrides `target_sum` and sets :math:`K = 4 \cdot α \cdot s` by + the delta method, where :math:`s` is the mean cell depth, calibrating the + count-scale pseudocount to the variance-stabilizing value ``y0 = 1/(4·α)`` + :cite:p:`Booeshaghi2022`. Pass ``"auto"`` to estimate :math:`α` from the + data (closed-form least squares of ``var = μ + α·μ²`` across genes). A + :class:`ValueError` is raised if the estimated or supplied :math:`α` is + not positive (e.g. underdispersed data); pass `target_sum` instead. layer Layer to normalize instead of `X`. inplace @@ -486,13 +511,11 @@ def normalize_clr( if isinstance(x, CSCBase): x = x.tocsr() if issubclass(x.dtype.type, int | np.integer): - x = x.astype(np.float32) + x = x.astype(np.float64) start = logg.info("normalizing counts per cell via shifted CLR") - x, cell_depths = _normalize_clr_helper( - x, c=c, target_sum=target_sum, alpha=alpha, scale=scale - ) + x, cell_depths = _normalize_clr_helper(x, target_sum=target_sum, alpha=alpha) if not isinstance(cell_depths, DaskArray) and not np.all(cell_depths > 0): warn("Some cells have zero counts", UserWarning) diff --git a/tests/test_normalization.py b/tests/test_normalization.py index 28884bb2c0..752a54f80e 100644 --- a/tests/test_normalization.py +++ b/tests/test_normalization.py @@ -346,24 +346,33 @@ def test_compute_nnz_median(array_type, dtype): ) -def _clr_reference(x, *, c=1.0, target_sum=None, alpha=None, scale=None) -> np.ndarray: +def _estimate_alpha_reference(x) -> float: + """Closed-form OLS overdispersion, independent of the implementation.""" + 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_sum=None, alpha=None) -> np.ndarray: """Self-contained dense shifted-CLR, independent of the implementation. - PF to a target depth, log(. + c), then subtract the per-cell mean. Empty - cells (zero depth) are left as all-zero rows, matching `normalize_clr`. + PF to a target depth, log1p, then subtract the per-cell mean. Empty cells + (zero depth) are left as all-zero rows, matching `normalize_clr`. """ x = np.asarray(to_ndarray(x), dtype=np.float64) depths = x.sum(axis=1) # Resolve the PF target depth K (paper's notation; exposed as `target_sum`). if alpha is not None: - if scale is None: - scale = depths.mean() - target_sum = 4.0 * alpha * scale + if alpha == "auto": + alpha = _estimate_alpha_reference(x) + target_sum = 4.0 * alpha * depths.mean() elif target_sum is None: target_sum = depths.mean() safe_depths = np.where(depths == 0, 1.0, depths) u = x * (target_sum / safe_depths)[:, None] - log_u = np.log(u + c) + log_u = np.log1p(u) return log_u - log_u.mean(axis=1, keepdims=True) @@ -387,8 +396,8 @@ def test_normalize_clr_values(array_type, dtype): @pytest.mark.parametrize("array_type", ARRAY_TYPES_MEM) @pytest.mark.parametrize( "kwargs", - [{}, {"target_sum": 1e4}, {"c": 2.0}, {"alpha": 0.5}], - ids=["default", "target_sum", "shift_c", "alpha"], + [{}, {"target_sum": 1e4}, {"alpha": 0.5}, {"alpha": "auto"}], + ids=["default", "target_sum", "alpha", "alpha_auto"], ) def test_normalize_clr_params(array_type, kwargs): adata = AnnData(array_type(X_clr).astype("float32")) @@ -420,6 +429,30 @@ def test_normalize_clr_alpha_overrides_target_sum(): ) +@pytest.mark.parametrize("array_type", ARRAY_TYPES_MEM) +def test_normalize_clr_alpha_auto(array_type): + """`alpha="auto"` estimates the overdispersion and matches an explicit alpha.""" + estimated = _estimate_alpha_reference(X_clr) + assert estimated > 0 + + auto = AnnData(array_type(X_clr).astype("float32")) + sc.pp.normalize_clr(auto, alpha="auto") + + explicit = AnnData(array_type(X_clr).astype("float32")) + sc.pp.normalize_clr(explicit, alpha=estimated) + np.testing.assert_allclose( + to_ndarray(auto.X), to_ndarray(explicit.X), 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): + """A non-positive `alpha` cannot derive K = 4*alpha*s and raises.""" + adata = AnnData(sparse.csr_matrix(X_clr)) # noqa: TID251 + with pytest.raises(ValueError, match=r"alpha.*positive"): + sc.pp.normalize_clr(adata, alpha=alpha) + + @pytest.mark.parametrize("array_type", ARRAY_TYPES_MEM) def test_normalize_clr_zero_cell(array_type): """An empty cell stays all-zero, stays finite, and triggers a warning.""" From 6eee9543727782f6dd7a2d204509bd6fd1347077 Mon Sep 17 00:00:00 2001 From: rwbaber <97636743+rwbaber@users.noreply.github.com> Date: Mon, 15 Jun 2026 21:37:46 +0100 Subject: [PATCH 07/11] feat(pp): support dask arrays in normalize_clr - Switch overdispersion estimation to dask-aware mean_var - Use map_blocks with gene rechunking for cell centering - Expand test suites to cover dense and sparse Dask ARRAY_TYPES --- src/scanpy/preprocessing/_normalization.py | 110 +++++++++++++-------- tests/test_normalization.py | 6 +- 2 files changed, 74 insertions(+), 42 deletions(-) diff --git a/src/scanpy/preprocessing/_normalization.py b/src/scanpy/preprocessing/_normalization.py index bda84c57fb..a16fe0ae1e 100644 --- a/src/scanpy/preprocessing/_normalization.py +++ b/src/scanpy/preprocessing/_normalization.py @@ -7,6 +7,7 @@ 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 .. import logging as logg from .._compat import CSBase, CSCBase, CSRBase, DaskArray, warn @@ -324,48 +325,72 @@ def normalize_total( # noqa: PLR0912 # ----------------------------------------------------------------------------- -def _estimate_overdispersion(x: np.ndarray | CSBase) -> float: +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 variance - over cells. The model is linear in :math:`α`, so the ordinary least-squares - solution is closed form + :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. + 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. """ - n_cells = x.shape[0] - col_sum = np.asarray(stats.sum(x, axis=0)).ravel() - if isinstance(x, CSBase): - x_sq = x.copy() - x_sq.data **= 2 - else: - x_sq = np.square(x) - col_sum_sq = np.asarray(stats.sum(x_sq, axis=0)).ravel() - - mu = col_sum / n_cells - var = col_sum_sq / n_cells - mu**2 + mu, var = mean_var(x, axis=0, correction=0) mu2 = mu**2 - denominator = float(np.sum(mu2 * mu2)) + 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 `alpha` or `target_sum` explicitly." ) raise ValueError(msg) - return float(np.sum((var - mu) * mu2) / denominator) + return float(numerator / denominator) + + +def _clr_log_center(u: np.ndarray | CSBase) -> np.ndarray: + """Apply log1p then per-cell mean-centering to PF-scaled counts; return dense. + + Operates on a full-width block of cells: every gene must be present so the + per-cell mean spans the whole row. Used directly for in-memory input, and + per row-chunk under :meth:`~dask.array.Array.map_blocks` (after the genes + have been rechunked into a single chunk). + """ + n_genes = u.shape[1] + if isinstance(u, CSBase): + # Sparse "offset trick": store log1p only at the nonzero entries. The + # structural zeros already hold log1p(0) = 0, so the per-cell mean is just + # the stored row sum / n_genes, and no densification is needed until the + # final centering step. + u.data = np.log1p(u.data) + per_cell_mean = np.asarray(u.sum(axis=1)).ravel() / n_genes + dense_log = u.toarray() + else: + u = np.asarray(u) + dense_log = np.log1p(u) + per_cell_mean = dense_log.mean(axis=1) + + # Additive centering (the CLR step): map each cell onto the zero-sum + # (Aitchison) hyperplane. This is what fills the zeros and densifies the matrix. + return dense_log - per_cell_mean[:, None] def _normalize_clr_helper( - x: np.ndarray | CSBase, + x: np.ndarray | CSBase | DaskArray, *, target_sum: float | None, alpha: float | Literal["auto"] | None, -) -> tuple[np.ndarray, np.ndarray]: +) -> tuple[np.ndarray | DaskArray, np.ndarray | DaskArray]: """Compute the shifted CLR (PFlog1pPF) transform backend logic. This function acts as the internal processing engine for `normalize_clr`, @@ -379,7 +404,10 @@ def _normalize_clr_helper( See `normalize_clr` for the meaning of the parameters. """ - cell_depths = np.asarray(stats.sum(x, axis=1)).ravel() + # 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() # Resolve `target_sum` (the paper's PF target constant K, Booeshaghi et al. # 2022): the depth the first proportional-fitting step scales every cell to. @@ -397,9 +425,10 @@ def _normalize_clr_helper( # Delta-method calibration: K = 4 * alpha * s, where s is the mean cell # depth. This sets the count-scale pseudocount to the variance-stabilizing # value y0 = 1 / (4 * alpha), and overrides any value passed as `target_sum`. - target_sum = 4.0 * alpha * cell_depths.mean() + target_sum = 4.0 * alpha * float(cell_depths.mean()) elif target_sum is None: - target_sum = cell_depths.mean() + # `float(...)` triggers a single blocking `.compute()` for dask input. + target_sum = float(cell_depths.mean()) # else: use the `target_sum` passed by the caller. # Proportional fitting: rescale each cell so its counts sum to `target_sum` (K). @@ -409,24 +438,18 @@ def _normalize_clr_helper( x, cell_depths / target_sum, op=truediv, axis=0, allow_divide_by_zero=False ) - n_genes = x.shape[1] - - if isinstance(u, CSBase): - # Sparse "offset trick": store log1p only at the nonzero entries. The - # structural zeros already hold log1p(0) = 0, so the per-cell mean is just - # the stored row sum / n_genes, and no densification is needed until the - # final centering step. - u.data = np.log1p(u.data) - per_cell_mean = np.asarray(u.sum(axis=1)).ravel() / n_genes - dense_log = u.toarray() + if isinstance(x, DaskArray): + # Per-cell centering needs the whole row, so collapse the genes into a + # single chunk; the transform is then embarrassingly parallel over the + # row-chunks. The output is always dense, hence the dense numpy `meta`. + u = u.rechunk({1: -1}) + dense_log = u.map_blocks( + _clr_log_center, dtype=np.float64, meta=np.array([], dtype=np.float64) + ) else: - u = np.asarray(u) - dense_log = np.log1p(u) - per_cell_mean = dense_log.mean(axis=1) + dense_log = _clr_log_center(u) - # Additive centering (the CLR step): map each cell onto the zero-sum - # (Aitchison) hyperplane. This is what fills the zeros and densifies the matrix. - return dense_log - per_cell_mean[:, None], cell_depths + return dense_log, cell_depths def normalize_clr( @@ -455,6 +478,15 @@ def normalize_clr( Because the per-cell centering fills the zero entries, the output is always dense, and each cell sums to exactly zero. + .. note:: + When used with a :class:`~dask.array.Array` in `adata.X`, deriving the + proportional-fitting target :math:`K` from the data requires a global + reduction and therefore triggers a blocking `.compute()`: once for the + default mean-depth target, and once more for `alpha` (including + ``alpha="auto"``). Only the scalar reduction is materialized, not the + matrix. Passing `target_sum` explicitly avoids this and keeps the whole + operation lazy. + Parameters ---------- adata diff --git a/tests/test_normalization.py b/tests/test_normalization.py index 752a54f80e..01a4856703 100644 --- a/tests/test_normalization.py +++ b/tests/test_normalization.py @@ -376,7 +376,7 @@ def _clr_reference(x, *, target_sum=None, alpha=None) -> np.ndarray: return log_u - log_u.mean(axis=1, keepdims=True) -@pytest.mark.parametrize("array_type", ARRAY_TYPES_MEM) +@pytest.mark.parametrize("array_type", ARRAY_TYPES) @pytest.mark.parametrize("dtype", ["float32", "int64"]) def test_normalize_clr_values(array_type, dtype): """Values match the reference and every cell sums to zero, for all layouts. @@ -393,7 +393,7 @@ def test_normalize_clr_values(array_type, dtype): np.testing.assert_allclose(result.sum(axis=1), 0.0, atol=1e-5) -@pytest.mark.parametrize("array_type", ARRAY_TYPES_MEM) +@pytest.mark.parametrize("array_type", ARRAY_TYPES) @pytest.mark.parametrize( "kwargs", [{}, {"target_sum": 1e4}, {"alpha": 0.5}, {"alpha": "auto"}], @@ -429,7 +429,7 @@ def test_normalize_clr_alpha_overrides_target_sum(): ) -@pytest.mark.parametrize("array_type", ARRAY_TYPES_MEM) +@pytest.mark.parametrize("array_type", ARRAY_TYPES) def test_normalize_clr_alpha_auto(array_type): """`alpha="auto"` estimates the overdispersion and matches an explicit alpha.""" estimated = _estimate_alpha_reference(X_clr) From b284927f63ae20b001b8d121611ae4f0dcc6f7d2 Mon Sep 17 00:00:00 2001 From: rwbaber <97636743+rwbaber@users.noreply.github.com> Date: Mon, 15 Jun 2026 23:03:12 +0100 Subject: [PATCH 08/11] test(pp): cover alpha="auto" zero-mean overdispersion error --- tests/test_normalization.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/tests/test_normalization.py b/tests/test_normalization.py index 01a4856703..0fc4caca5d 100644 --- a/tests/test_normalization.py +++ b/tests/test_normalization.py @@ -453,6 +453,13 @@ def test_normalize_clr_nonpositive_alpha_raises(alpha): sc.pp.normalize_clr(adata, alpha=alpha) +def test_normalize_clr_alpha_auto_zero_mean_raises(): + """`alpha="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, alpha="auto") + + @pytest.mark.parametrize("array_type", ARRAY_TYPES_MEM) def test_normalize_clr_zero_cell(array_type): """An empty cell stays all-zero, stays finite, and triggers a warning.""" From 3549f1921bfa5b1a7c51ac00487c3cce16ca6cc9 Mon Sep 17 00:00:00 2001 From: rwbaber <97636743+rwbaber@users.noreply.github.com> Date: Wed, 8 Jul 2026 07:48:41 +0100 Subject: [PATCH 09/11] Add sparse PFlog normalization and PCA support --- docs/references.bib | 2 +- docs/release-notes/4160.feat.md | 2 +- src/scanpy/preprocessing/_normalization.py | 292 +++++++++++---------- src/scanpy/preprocessing/_pca/__init__.py | 32 ++- src/scanpy/preprocessing/_pca/_compat.py | 81 ++++++ tests/test_normalization.py | 152 +++++++---- tests/test_package_structure.py | 7 +- tests/test_pca.py | 81 ++++++ 8 files changed, 458 insertions(+), 191 deletions(-) diff --git a/docs/references.bib b/docs/references.bib index 906747e0c2..0b1453f989 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -113,7 +113,7 @@ @article{Blondel2008 @article{Booeshaghi2022, author = {Booeshaghi, A. Sina and Hallgrímsdóttir, Ingileif B. and Gálvez-Merchán, Ángel and Pachter, Lior}, - title = {Depth normalization for single-cell genomics count data}, + 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}, diff --git a/docs/release-notes/4160.feat.md b/docs/release-notes/4160.feat.md index 02b93ce488..2c12d631db 100644 --- a/docs/release-notes/4160.feat.md +++ b/docs/release-notes/4160.feat.md @@ -1 +1 @@ -Add {func}`scanpy.pp.normalize_clr` for shifted centered log-ratio (PFlog1pPF) normalization, a variance-stabilizing, depth-invariant and rank-preserving count transform {cite:p}`Booeshaghi2022` {smaller}`R Baber` +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/_normalization.py b/src/scanpy/preprocessing/_normalization.py index a16fe0ae1e..51da1f9354 100644 --- a/src/scanpy/preprocessing/_normalization.py +++ b/src/scanpy/preprocessing/_normalization.py @@ -8,6 +8,7 @@ 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 @@ -309,22 +310,6 @@ def normalize_total( # noqa: PLR0912 return None -# ----------------------------------------------------------------------------- -# Shifted Centered Log-Ratio (PFlogPF) Core Implementation -# -# The mathematical logic and sparse matrix optimization ("offset trick") -# implemented below are adapted from the Pachter Lab reference repository: -# https://github.com/pachterlab/bhgp_2022 -# -# Reference Publication: -# Booeshaghi et al. (2022, 2026) "Depth normalization for single-cell genomics -# count data" bioRxiv. doi: 10.1101/2022.05.06.490859 -# -# Copyright (c) 2022, Pachter Lab. All rights reserved. -# Licensed under the BSD 2-Clause License. -# ----------------------------------------------------------------------------- - - def _estimate_overdispersion(x: np.ndarray | CSBase | DaskArray) -> float: r"""Estimate the negative-binomial overdispersion :math:`α` from raw counts. @@ -352,161 +337,168 @@ def _estimate_overdispersion(x: np.ndarray | CSBase | DaskArray) -> float: if denominator == 0.0: msg = ( "Cannot estimate overdispersion: every gene has zero mean. " - "Pass `alpha` or `target_sum` explicitly." + "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 float(numerator / denominator) + return alpha -def _clr_log_center(u: np.ndarray | CSBase) -> np.ndarray: - """Apply log1p then per-cell mean-centering to PF-scaled counts; return dense. +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) - Operates on a full-width block of cells: every gene must be present so the - per-cell mean spans the whole row. Used directly for in-memory input, and - per row-chunk under :meth:`~dask.array.Array.map_blocks` (after the genes - have been rechunked into a single chunk). - """ - n_genes = u.shape[1] - if isinstance(u, CSBase): - # Sparse "offset trick": store log1p only at the nonzero entries. The - # structural zeros already hold log1p(0) = 0, so the per-cell mean is just - # the stored row sum / n_genes, and no densification is needed until the - # final centering step. - u.data = np.log1p(u.data) - per_cell_mean = np.asarray(u.sum(axis=1)).ravel() / n_genes - dense_log = u.toarray() - else: - u = np.asarray(u) - dense_log = np.log1p(u) - per_cell_mean = dense_log.mean(axis=1) - # Additive centering (the CLR step): map each cell onto the zero-sum - # (Aitchison) hyperplane. This is what fills the zeros and densifies the matrix. - return dense_log - per_cell_mean[:, None] +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( +def _normalize_clr_helper( # noqa: PLR0912 x: np.ndarray | CSBase | DaskArray, *, - target_sum: float | None, - alpha: float | Literal["auto"] | None, -) -> tuple[np.ndarray | DaskArray, np.ndarray | DaskArray]: - """Compute the shifted CLR (PFlog1pPF) transform backend logic. - - This function acts as the internal processing engine for `normalize_clr`, - adapted from the reference implementations provided by the Pachter Lab - (Booeshaghi et al. 2022). It applies initial depth-scaling followed by - an optimized sparse log-shift operation before mapping coordinates onto - the zero-sum Aitchison hyperplane. - - Returns the dense, per-cell-centered matrix and the raw per-cell depths - (the latter is used by the caller to warn about empty cells). - - See `normalize_clr` for the meaning of the parameters. - """ + 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() - # Resolve `target_sum` (the paper's PF target constant K, Booeshaghi et al. - # 2022): the depth the first proportional-fitting step scales every cell to. + mean_depth = ( + float(cell_depths.mean().compute()) + if isinstance(cell_depths, DaskArray) + else float(cell_depths.mean()) + ) if alpha is not None: - if alpha == "auto": - # Estimate the overdispersion from the raw counts (closed-form OLS, - # as in the reference implementation). - alpha = _estimate_overdispersion(x) if not alpha > 0: msg = ( - f"`alpha` must be positive to derive K = 4 * alpha * s, got {alpha}. " - "The data may be underdispersed; pass `target_sum` explicitly instead." + f"`alpha` must be positive to compute PFlog, got {alpha}. " + "The data may be underdispersed." ) raise ValueError(msg) - # Delta-method calibration: K = 4 * alpha * s, where s is the mean cell - # depth. This sets the count-scale pseudocount to the variance-stabilizing - # value y0 = 1 / (4 * alpha), and overrides any value passed as `target_sum`. - target_sum = 4.0 * alpha * float(cell_depths.mean()) - elif target_sum is None: - # `float(...)` triggers a single blocking `.compute()` for dask input. - target_sum = float(cell_depths.mean()) - # else: use the `target_sum` passed by the caller. - - # Proportional fitting: rescale each cell so its counts sum to `target_sum` (K). - # Dividing by `cell_depths / target_sum` (with allow_divide_by_zero=False) - # leaves empty cells as all-zero rows instead of producing invalid infinities. - u = axis_mul_or_truediv( - x, cell_depths / target_sum, op=truediv, axis=0, allow_divide_by_zero=False - ) + 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": + depths = ( + cell_depths.compute() + if isinstance(cell_depths, DaskArray) + else cell_depths + ) + target_sum = float(np.median(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 isinstance(x, DaskArray): - # Per-cell centering needs the whole row, so collapse the genes into a - # single chunk; the transform is then embarrassingly parallel over the - # row-chunks. The output is always dense, hence the dense numpy `meta`. - u = u.rechunk({1: -1}) - dense_log = u.map_blocks( - _clr_log_center, dtype=np.float64, meta=np.array([], dtype=np.float64) + if alpha is not None: + x = x * scale + log_values = x.map_blocks( + _log1p_sparse_block, dtype=np.float64, meta=x._meta.astype(np.float64) ) else: - dense_log = _clr_log_center(u) - - return dense_log, cell_depths + if alpha is not None: + x = x * scale + 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_sum: float | None = None, - alpha: float | Literal["auto"] | None = None, + 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] | None: - r"""Normalize counts with the shifted centered log-ratio (PFlog1pPF) transform. +) -> AnnData | dict[str, np.ndarray | CSBase | DaskArray | dict[str, object]] | None: + r"""Normalize counts with the shifted centered log-ratio (PFlog) transform. - Computes the shifted centered log-ratio (CLR) transform + With `target="auto"` (or explicit `alpha`), computes PFlog as .. math:: - T(x)_i = \log(u_i + 1) - \frac{1}{D} \sum_{j=1}^D \log(u_j + 1), + T(x)_i = \log(1 + 4 α x_i) + - \frac{1}{D} \sum_{j=1}^D \log(1 + 4 α x_j), - where :math:`u_i = K \, x_i / \sum_j x_j` are the depth-normalized counts - (proportional fitting to a target depth :math:`K`) and :math:`D` is the - number of genes. Equivalently this is proportional fitting, then - :obj:`~numpy.log1p`, then per-cell mean-centering in log space (the - centered-log-ratio step). This count transform is simultaneously - variance-stabilizing, depth-invariant, and rank-preserving :cite:p:`Booeshaghi2022`. - - Because the per-cell centering fills the zero entries, the output is always - dense, and each cell sums to exactly zero. + 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:: - When used with a :class:`~dask.array.Array` in `adata.X`, deriving the - proportional-fitting target :math:`K` from the data requires a global - reduction and therefore triggers a blocking `.compute()`: once for the - default mean-depth target, and once more for `alpha` (including - ``alpha="auto"``). Only the scalar reduction is materialized, not the - matrix. Passing `target_sum` explicitly avoids this and keeps the whole - operation lazy. + `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_sum - Target depth :math:`K` for the first proportional-fitting step. If `None` - (and `alpha` is not given), the empirical mean cell depth is used. Unlike - :func:`~scanpy.pp.normalize_total`, this is only an *intermediate* target: - cells sum to `target_sum` after proportional fitting, but the subsequent - log and centering steps make each cell sum to exactly zero in the output. + 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_sum` and sets :math:`K = 4 \cdot α \cdot s` by - the delta method, where :math:`s` is the mean cell depth, calibrating the - count-scale pseudocount to the variance-stabilizing value ``y0 = 1/(4·α)`` - :cite:p:`Booeshaghi2022`. Pass ``"auto"`` to estimate :math:`α` from the - data (closed-form least squares of ``var = μ + α·μ²`` across genes). A - :class:`ValueError` is raised if the estimated or supplied :math:`α` is - not positive (e.g. underdispersed data); pass `target_sum` instead. + 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 @@ -518,17 +510,17 @@ def normalize_clr( Returns ------- - Returns a dictionary with the normalized copy of `adata.X` (key `"X"`) or - updates `adata` with the normalized matrix, depending on `inplace`. + 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, 3], [4, 5, 6]], dtype="float32")) - >>> sc.pp.normalize_clr(adata) - >>> np.allclose(adata.X.sum(axis=1), 0, atol=1e-6) + >>> 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: @@ -542,19 +534,49 @@ def normalize_clr( 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 shifted CLR") + start = logg.info("normalizing counts per cell via PFlog") - x, cell_depths = _normalize_clr_helper(x, target_sum=target_sum, alpha=alpha) + 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) - dat = dict(X=x) + 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( + report, + method="PFlog", + encoding_type="shifted_clr", + row_center_key=row_center_key, + layer=layer, + densify=densify, + ) + dat = dict( + X=dense_x if densify else log_values, + row_center=row_center_obs, + metadata=metadata, + ) if inplace: - _set_obs_rep(adata, dat["X"], layer=layer) + 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})", diff --git a/src/scanpy/preprocessing/_pca/__init__.py b/src/scanpy/preprocessing/_pca/__init__.py index 50fb6f4c1a..c3cdca883b 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,17 @@ 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): + row_center_key = adata.uns[layer].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 +282,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 +329,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 0fc4caca5d..fc4eb94a4a 100644 --- a/tests/test_normalization.py +++ b/tests/test_normalization.py @@ -18,6 +18,7 @@ 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 ( @@ -337,7 +338,7 @@ def test_compute_nnz_median(array_type, dtype): # ------------------------------------------------------------------------------ -# normalize_clr (shifted CLR / PFlog1pPF) +# normalize_clr (shifted CLR / PFlog) # ------------------------------------------------------------------------------ # A small count matrix with no empty cells, used for the value/equivalence tests. @@ -347,7 +348,7 @@ def test_compute_nnz_median(array_type, dtype): def _estimate_alpha_reference(x) -> float: - """Closed-form OLS overdispersion, independent of the implementation.""" + """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 @@ -355,120 +356,132 @@ def _estimate_alpha_reference(x) -> float: return float(np.sum((var - mu) * mu2) / np.sum(mu2 * mu2)) -def _clr_reference(x, *, target_sum=None, alpha=None) -> np.ndarray: - """Self-contained dense shifted-CLR, independent of the implementation. +def _clr_reference(x, *, target="auto", alpha=None) -> np.ndarray: + """Calculate shifted CLR densely for reference. - PF to a target depth, log1p, then subtract the per-cell mean. Empty cells - (zero depth) are left as all-zero rows, matching `normalize_clr`. + 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) - # Resolve the PF target depth K (paper's notation; exposed as `target_sum`). if alpha is not None: - if alpha == "auto": - alpha = _estimate_alpha_reference(x) - target_sum = 4.0 * alpha * depths.mean() - elif target_sum is None: - target_sum = depths.mean() - safe_depths = np.where(depths == 0, 1.0, depths) - u = x * (target_sum / safe_depths)[:, None] - log_u = np.log1p(u) + 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) -@pytest.mark.parametrize("array_type", ARRAY_TYPES) +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): - """Values match the reference and every cell sums to zero, for all layouts. - - Asserting that numpy / csr / csc / sparse-array inputs all equal the same - dense reference also proves the sparse "offset trick" matches the dense path. - """ + """Check values against the reference and zero-sum cells.""" adata = AnnData(array_type(X_clr).astype(dtype)) sc.pp.normalize_clr(adata) - result = to_ndarray(adata.X) + 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" -@pytest.mark.parametrize("array_type", ARRAY_TYPES) +@pytest.mark.parametrize("array_type", ARRAY_TYPES_MEM) @pytest.mark.parametrize( "kwargs", - [{}, {"target_sum": 1e4}, {"alpha": 0.5}, {"alpha": "auto"}], - ids=["default", "target_sum", "alpha", "alpha_auto"], + [{}, {"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( - to_ndarray(adata.X), _clr_reference(X_clr, **kwargs), rtol=1e-5, atol=1e-5 + _reconstruct_clr(adata), _clr_reference(X_clr, **kwargs), rtol=1e-5, atol=1e-5 ) -def test_normalize_clr_alpha_overrides_target_sum(): - """`alpha` sets target_sum = 4*alpha*scale and overrides any given `target_sum`.""" +def test_normalize_clr_alpha_is_constant_scale_and_overrides_target(): + """`alpha` stores uncentered log1p(4 * alpha * x), independent of depth.""" alpha = 0.5 - scale = X_clr.sum(axis=1).mean() via_alpha = AnnData(sparse.csr_matrix(X_clr)) # noqa: TID251 sc.pp.normalize_clr(via_alpha, alpha=alpha) - - via_target = AnnData(sparse.csr_matrix(X_clr)) # noqa: TID251 - sc.pp.normalize_clr(via_target, target_sum=4.0 * alpha * scale) np.testing.assert_allclose( - to_ndarray(via_alpha.X), to_ndarray(via_target.X), rtol=1e-5, atol=1e-5 + via_alpha.layers["pflog"].toarray(), + np.log1p(4.0 * alpha * X_clr), + rtol=1e-5, + atol=1e-5, ) - # passing both -> alpha wins, target_sum ignored both = AnnData(sparse.csr_matrix(X_clr)) # noqa: TID251 - sc.pp.normalize_clr(both, alpha=alpha, target_sum=999.0) + sc.pp.normalize_clr(both, alpha=alpha, target=999.0) np.testing.assert_allclose( - to_ndarray(both.X), to_ndarray(via_alpha.X), rtol=1e-5, atol=1e-5 + _reconstruct_clr(both), _reconstruct_clr(via_alpha), rtol=1e-5, atol=1e-5 ) -@pytest.mark.parametrize("array_type", ARRAY_TYPES) +@pytest.mark.parametrize("array_type", ARRAY_TYPES_MEM) def test_normalize_clr_alpha_auto(array_type): - """`alpha="auto"` estimates the overdispersion and matches an explicit alpha.""" + """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, alpha="auto") + sc.pp.normalize_clr(auto) explicit = AnnData(array_type(X_clr).astype("float32")) sc.pp.normalize_clr(explicit, alpha=estimated) np.testing.assert_allclose( - to_ndarray(auto.X), to_ndarray(explicit.X), rtol=1e-5, atol=1e-5 + _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): - """A non-positive `alpha` cannot derive K = 4*alpha*s and raises.""" + """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(): - """`alpha="auto"` cannot estimate overdispersion when every gene mean is zero.""" + """`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, alpha="auto") + sc.pp.normalize_clr(adata) @pytest.mark.parametrize("array_type", ARRAY_TYPES_MEM) def test_normalize_clr_zero_cell(array_type): - """An empty cell stays all-zero, stays finite, and triggers a warning.""" + """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 = to_ndarray(adata.X) + result = _reconstruct_clr(adata) assert np.isfinite(result).all() np.testing.assert_allclose(result[1], 0.0, atol=1e-6) @@ -479,7 +492,12 @@ def test_normalize_clr_inplace_false(): out = sc.pp.normalize_clr(adata, inplace=False) assert isinstance(out, dict) - np.testing.assert_allclose(out["X"], _clr_reference(X_clr), rtol=1e-5, atol=1e-5) + 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) @@ -492,7 +510,7 @@ def test_normalize_clr_copy(): assert isinstance(returned, AnnData) assert returned is not adata np.testing.assert_allclose( - to_ndarray(returned.X), _clr_reference(X_clr), rtol=1e-5, atol=1e-5 + _reconstruct_clr(returned), _clr_reference(X_clr), rtol=1e-5, atol=1e-5 ) # original is left untouched assert isinstance(adata.X, CSBase) @@ -507,7 +525,7 @@ def test_normalize_clr_copy_inplace_error(): def test_normalize_clr_layer(): - """`layer` targets that layer and leaves `X` untouched.""" + """`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 @@ -517,13 +535,49 @@ def test_normalize_clr_layer(): np.testing.assert_array_equal(to_ndarray(adata.X), x_before) np.testing.assert_allclose( - to_ndarray(adata.layers["counts"]), + _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 + ) + + def test_normalize_clr_view(): adata = AnnData(X_clr.copy()) v = adata[:, :] diff --git a/tests/test_package_structure.py b/tests/test_package_structure.py index 95789eca00..c077dd32a0 100644 --- a/tests/test_package_structure.py +++ b/tests/test_package_structure.py @@ -94,11 +94,12 @@ 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"] = copy_sigs["sc.pp.normalize_clr"]["return_ann"] = ( - "AnnData | dict[str, np.ndarray] | None" -) +]["return_ann"] = "AnnData | dict[str, np.ndarray] | None" copy_sigs["sc.external.pp.magic"]["copy_default"] = 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"] From 30a7f4f0172d94773a9d381bd8ff99cc42723761 Mon Sep 17 00:00:00 2001 From: rwbaber <97636743+rwbaber@users.noreply.github.com> Date: Wed, 8 Jul 2026 19:53:58 +0100 Subject: [PATCH 10/11] Refine PFlog metadata and Dask target handling --- src/scanpy/preprocessing/_normalization.py | 26 ++++++++++------------ src/scanpy/preprocessing/_pca/__init__.py | 4 +++- tests/test_normalization.py | 11 +++++++++ 3 files changed, 26 insertions(+), 15 deletions(-) diff --git a/src/scanpy/preprocessing/_normalization.py b/src/scanpy/preprocessing/_normalization.py index 51da1f9354..a330203a67 100644 --- a/src/scanpy/preprocessing/_normalization.py +++ b/src/scanpy/preprocessing/_normalization.py @@ -408,12 +408,14 @@ def _normalize_clr_helper( # noqa: PLR0912 if target == "mean": target_sum = mean_depth elif target == "median": - depths = ( - cell_depths.compute() - if isinstance(cell_depths, DaskArray) - else cell_depths - ) - target_sum = float(np.median(depths)) + 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) @@ -427,15 +429,14 @@ def _normalize_clr_helper( # noqa: PLR0912 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): - if alpha is not None: - x = x * scale log_values = x.map_blocks( _log1p_sparse_block, dtype=np.float64, meta=x._meta.astype(np.float64) ) else: - if alpha is not None: - x = x * scale log_values = _log1p_sparse_block(x) row_center = stats.sum(log_values, axis=1) / x.shape[1] @@ -559,12 +560,9 @@ def normalize_clr( dense_x = _densify_shifted_clr(log_values, row_center) if densify else None row_center_key = f"{key_added}_center" metadata = dict( - report, - method="PFlog", encoding_type="shifted_clr", row_center_key=row_center_key, - layer=layer, - densify=densify, + params=report | {"layer": layer, "densify": densify}, ) dat = dict( X=dense_x if densify else log_values, diff --git a/src/scanpy/preprocessing/_pca/__init__.py b/src/scanpy/preprocessing/_pca/__init__.py index c3cdca883b..8f5a4b86e1 100644 --- a/src/scanpy/preprocessing/_pca/__init__.py +++ b/src/scanpy/preprocessing/_pca/__init__.py @@ -245,7 +245,9 @@ def pca( # noqa: PLR0912, PLR0913, PLR0915 row_center_key = None if return_anndata and layer is not None and isinstance(adata.uns.get(layer), dict): - row_center_key = adata.uns[layer].get("row_center_key") + 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`." diff --git a/tests/test_normalization.py b/tests/test_normalization.py index fc4eb94a4a..a032f60b6a 100644 --- a/tests/test_normalization.py +++ b/tests/test_normalization.py @@ -406,6 +406,7 @@ def test_normalize_clr_values(array_type, dtype): 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) @@ -578,6 +579,16 @@ def test_normalize_clr_dask(sparse_blocks, densify): ) +@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[:, :] From 3fd7086f66476e2a14e85ce37bc503929b390c05 Mon Sep 17 00:00:00 2001 From: rwbaber <97636743+rwbaber@users.noreply.github.com> Date: Wed, 8 Jul 2026 20:54:23 +0100 Subject: [PATCH 11/11] Fix stale AnnData concatenation tutorial link --- docs/release-notes/1.6.0.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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`