diff --git a/docs/api/preprocessing.md b/docs/api/preprocessing.md index 1ac8181ba2..ebd56654ca 100644 --- a/docs/api/preprocessing.md +++ b/docs/api/preprocessing.md @@ -29,6 +29,7 @@ For visual quality control, see {func}`~scanpy.pl.highest_expr_genes` and pp.log1p pp.pca pp.normalize_total + pp.normalize_clr pp.regress_out pp.scale pp.sample diff --git a/docs/references.bib b/docs/references.bib index aa2261fe7b..0b1453f989 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -111,6 +111,16 @@ @article{Blondel2008 pages = {P10008}, } +@article{Booeshaghi2022, + author = {Booeshaghi, A. Sina and Hallgrímsdóttir, Ingileif B. and Gálvez-Merchán, Ángel and Pachter, Lior}, + title = {Normalization for sampled count data}, + year = {2026}, + url = {https://doi.org/10.1101/2022.05.06.490859}, + doi = {10.1101/2022.05.06.490859}, + publisher = {Cold Spring Harbor Laboratory}, + journal = {bioRxiv}, +} + @article{Burczynski2006, author = {Burczynski, Michael E. and Peterson, Ron L. and Twine, Natalie C. and Zuberek, Krystyna A. and Brodeur, Brendan J. and Casciotti, Lori and Maganti, Vasu and Reddy, Padma S. and Strahs, Andrew and Immermann, Fred and Spinelli, Walter and Schwertschlag, Ulrich and Slager, Anna M. and Cotreau, Monette M. and Dorner, Andrew J.}, title = {Molecular Classification of Crohn’s Disease and Ulcerative Colitis Patients Using Transcriptional Profiles in Peripheral Blood Mononuclear Cells}, diff --git a/docs/release-notes/4160.feat.md b/docs/release-notes/4160.feat.md new file mode 100644 index 0000000000..0f6c30457b --- /dev/null +++ b/docs/release-notes/4160.feat.md @@ -0,0 +1 @@ +Add {func}`scanpy.pp.normalize_clr` for PFlog shifted centered log-ratio normalization, a variance-stabilizing, depth-invariant and rank-preserving count transform {cite:p}`Booeshaghi2022`. By default, it estimates negative-binomial overdispersion and uses the Anscombe pseudocount. The complete centered result is stored in `X` or the selected layer. Because CLR centering generally maps zeros to non-zero values, the output is dense. {smaller}`R Baber` diff --git a/src/scanpy/preprocessing/__init__.py b/src/scanpy/preprocessing/__init__.py index 10641e730d..2544c4aa64 100644 --- a/src/scanpy/preprocessing/__init__.py +++ b/src/scanpy/preprocessing/__init__.py @@ -9,7 +9,7 @@ from ._harmony import harmony_integrate from ._hashsolo import hashsolo 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 @@ -37,6 +37,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 d0ce3ac175..63987e042e 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 @@ -304,3 +305,187 @@ def normalize_total( # noqa: PLR0912 elif not inplace: return dat return None + + +def _estimate_overdispersion(x: np.ndarray | CSBase | DaskArray) -> float: + r"""Estimate the negative-binomial overdispersion :math:`α` from raw counts. + + Fits :math:`\mathrm{Var}_g = μ_g + α \cdot μ_g^2` across genes, where + :math:`μ_g` and :math:`\mathrm{Var}_g` are the per-gene mean and (population) + variance over cells. The model is linear in :math:`α`, so the ordinary + least-squares solution is closed form + + .. math:: + α = \frac{\sum_g (\mathrm{Var}_g - μ_g) \, μ_g^2}{\sum_g μ_g^4}, + + which is exactly the minimizer a non-linear `curve_fit` would converge to, + but without the dependency. :func:`~fast_array_utils.stats.mean_var` is + dispatched for dense, sparse and dask input alike, so the only dask-specific + step is computing the two final scalar sums. + """ + mu, var = mean_var(x, axis=0, correction=0) + mu2 = mu**2 + numerator = np.sum((var - mu) * mu2) + denominator = np.sum(mu2 * mu2) + if isinstance(x, DaskArray): + import dask + + numerator, denominator = dask.compute(numerator, denominator) + if denominator == 0.0: + msg = ( + "Cannot estimate overdispersion: every gene has zero mean. " + "Pass a positive `alpha` explicitly." + ) + raise ValueError(msg) + alpha = float(numerator / denominator) + if not alpha > 0: + msg = ( + f"Estimated overdispersion is non-positive (alpha = {alpha}); " + "pass a positive `alpha` explicitly." + ) + raise ValueError(msg) + return alpha + + +def _log1p_sparse_block(x: np.ndarray | CSBase) -> np.ndarray | CSBase: + """Apply log1p to a dense or sparse block while preserving sparse zeros.""" + if isinstance(x, CSBase): + x = x.copy() + x.data = np.log1p(x.data) + return x + return np.log1p(x) + + +def _normalize_clr_helper( + x: np.ndarray | CSBase | DaskArray, + *, + alpha: float | None, +) -> tuple[np.ndarray | DaskArray, np.ndarray | DaskArray]: + """Compute the dense PFlog / shifted-CLR matrix and cell depths.""" + # 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() + + if alpha is None: + alpha = _estimate_overdispersion(x) + elif not alpha > 0: + msg = ( + f"`alpha` must be positive to compute PFlog, got {alpha}. " + "The data may be underdispersed." + ) + raise ValueError(msg) + + x = x * (4.0 * float(alpha)) + + if isinstance(x, DaskArray): + log_values = x.map_blocks( + _log1p_sparse_block, dtype=np.float64, meta=x._meta.astype(np.float64) + ) + else: + log_values = _log1p_sparse_block(x) + + row_center = stats.sum(log_values, axis=1) / x.shape[1] + if not isinstance(row_center, DaskArray): + row_center = np.asarray(row_center).ravel() + if isinstance(log_values, CSBase): + log_values = log_values.toarray() + return log_values - row_center[:, None], cell_depths + + +def normalize_clr( + adata: AnnData, + *, + alpha: float | None = None, + layer: str | None = None, + inplace: bool = True, + copy: bool = False, +) -> AnnData | dict[str, np.ndarray | DaskArray] | None: + r"""Normalize counts with the shifted centered log-ratio (PFlog) transform. + + If `alpha` is not provided, it is estimated from the input matrix. PFlog is + then computed as + + .. math:: + T(x)_i = \log(1 + 4 α x_i) + - \frac{1}{D} \sum_{j=1}^D \log(1 + 4 α x_j), + + which is equivalent to centering + :math:`\log(x_i + 1 / (4 α))` because the constant :math:`\log(4 α)` + cancels during CLR centering. + + .. note:: + CLR centering generally maps zeros to non-zero values, so the resulting + matrix is dense even when the input is sparse. + + Parameters + ---------- + adata + The annotated data matrix of shape `n_obs` × `n_vars`. + Rows correspond to cells and columns to genes. + alpha + Negative-binomial overdispersion of the dataset (``var = μ + α·μ²``). + If `None`, it is estimated from the input matrix. A positive numeric + value uses that value directly. PFlog applies + ``log1p(4 * alpha * x)`` before CLR centering. + 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 matrix or updates `adata`, depending + on `inplace`. + + Example + ------- + >>> import numpy as np + >>> from anndata import AnnData + >>> import scanpy as sc + >>> adata = AnnData(np.array([[1, 2, 30], [4, 50, 6]], dtype="float32")) + >>> sc.pp.normalize_clr(adata, alpha=0.5) + >>> np.allclose(adata.X.sum(axis=1), 0, atol=1e-5) + 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_arr(adata, layer=layer) + if isinstance(x, CSCBase): + x = x.tocsr() + if not inplace: + x = x.copy() + if issubclass(x.dtype.type, int | np.integer): + x = x.astype(np.float64) + + start = logg.info("normalizing counts per cell via PFlog") + + x, cell_depths = _normalize_clr_helper(x, 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) + if inplace: + _set_obs_rep(adata, 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 1058b737e9..5b276ca898 100644 --- a/tests/test_normalization.py +++ b/tests/test_normalization.py @@ -11,15 +11,21 @@ from scipy import sparse import scanpy as sc +from scanpy._compat import CSBase, DaskArray from scanpy.preprocessing._normalization import _compute_nnz_median from testing.scanpy._helpers import ( _check_check_values_warnings, check_rep_mutation, check_rep_results, ) +from testing.scanpy._pytest.marks import needs # TODO: Add support for sparse-in-dask -from testing.scanpy._pytest.params import ARRAY_TYPES, ARRAY_TYPES_DENSE +from testing.scanpy._pytest.params import ( + ARRAY_TYPES, + ARRAY_TYPES_DENSE, + ARRAY_TYPES_MEM, +) if TYPE_CHECKING: from collections.abc import Callable @@ -350,3 +356,205 @@ def test_normalize_total_target_sum_ignores_zero_count_cells(array_type): # median of the non-zero row sums (10, 20, 30) is 20, not 15 np.testing.assert_allclose(stats.sum(adata.X, axis=1)[1:], 20.0) assert_equal(conv.to_dense(adata.X), conv.to_dense(expected.X)) + + +# ------------------------------------------------------------------------------ +# normalize_clr (shifted CLR / PFlog) +# ------------------------------------------------------------------------------ + +# A small count matrix with no empty cells, used for the value/equivalence tests. +X_clr = np.array( + [[5, 0, 3, 2], [1, 1, 0, 4], [0, 7, 2, 1], [3, 3, 3, 3]], dtype="float32" +) + + +def _estimate_alpha_reference(x) -> float: + """Calculate OLS overdispersion for reference.""" + x = np.asarray(to_ndarray(x), dtype=np.float64) + mu = x.mean(axis=0) + var = (x**2).mean(axis=0) - mu**2 + mu2 = mu**2 + return float(np.sum((var - mu) * mu2) / np.sum(mu2 * mu2)) + + +def _clr_reference(x, *, alpha=None) -> np.ndarray: + """Calculate PFlog densely for reference.""" + x = np.asarray(to_ndarray(x), dtype=np.float64) + if alpha is None: + alpha = _estimate_alpha_reference(x) + log_u = np.log1p(4.0 * alpha * x) + return log_u - log_u.mean(axis=1, keepdims=True) + + +def _materialize(x): + if hasattr(x, "compute"): + x = x.compute() + return to_ndarray(x) + + +@pytest.mark.parametrize("array_type", ARRAY_TYPES_MEM) +@pytest.mark.parametrize("dtype", ["float32", "int64"]) +def test_normalize_clr_values(array_type, dtype): + """Check values against the reference and zero-sum cells.""" + adata = AnnData(array_type(X_clr).astype(dtype)) + sc.pp.normalize_clr(adata) + result = _materialize(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) +def test_normalize_clr_explicit_alpha(array_type): + alpha = 0.5 + adata = AnnData(array_type(X_clr).astype("float32")) + sc.pp.normalize_clr(adata, alpha=alpha) + np.testing.assert_allclose( + _materialize(adata.X), + _clr_reference(X_clr, alpha=alpha), + rtol=1e-5, + atol=1e-5, + ) + + +@pytest.mark.parametrize("array_type", ARRAY_TYPES_MEM) +def test_normalize_clr_estimated_alpha_matches_explicit(array_type): + """Check that the default estimate matches explicit `alpha`.""" + estimated = _estimate_alpha_reference(X_clr) + assert estimated > 0 + + estimated_adata = AnnData(array_type(X_clr).astype("float32")) + sc.pp.normalize_clr(estimated_adata) + + explicit = AnnData(array_type(X_clr).astype("float32")) + sc.pp.normalize_clr(explicit, alpha=estimated) + np.testing.assert_allclose( + _materialize(estimated_adata.X), + _materialize(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): + """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_estimated_alpha_zero_mean_raises(): + """Alpha cannot be estimated when every gene mean is zero.""" + adata = AnnData(np.zeros((3, 4), dtype="float32")) + with pytest.raises(ValueError, match="Cannot estimate overdispersion"): + sc.pp.normalize_clr(adata) + + +@pytest.mark.parametrize("array_type", ARRAY_TYPES_MEM) +def test_normalize_clr_zero_cell(array_type): + """Keep an empty cell finite and all-zero.""" + x = X_clr.copy() + x[1] = 0 # make the second cell empty + adata = AnnData(array_type(x)) + with pytest.warns(UserWarning, match="Some cells have zero counts"): + sc.pp.normalize_clr(adata) + result = _materialize(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( + _materialize(out["X"]), _clr_reference(X_clr), rtol=1e-5, atol=1e-5 + ) + assert set(out) == {"X"} + # 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(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` selects the input layer and leaves `X` untouched.""" + adata = AnnData( + sparse.csr_matrix(X_clr), # noqa: TID251 + layers={"counts": sparse.csr_matrix(X_clr)}, # noqa: TID251 + ) + x_before = to_ndarray(adata.X).copy() + sc.pp.normalize_clr(adata, layer="counts") + + np.testing.assert_array_equal(to_ndarray(adata.X), x_before) + np.testing.assert_allclose( + adata.layers["counts"], + _clr_reference(X_clr), + rtol=1e-5, + atol=1e-5, + ) + assert isinstance(adata.layers["counts"], np.ndarray) + + +def test_normalize_clr_densifies_sparse_input(): + adata = AnnData(sparse.csr_matrix(X_clr)) # noqa: TID251 + sc.pp.normalize_clr(adata) + + assert isinstance(adata.X, np.ndarray) + np.testing.assert_allclose(adata.X, _clr_reference(X_clr), rtol=1e-5) + + +@needs.dask +@pytest.mark.parametrize("alpha", [None, 0.5], ids=["estimated", "explicit"]) +@pytest.mark.parametrize( + "sparse_blocks", [False, True], ids=["dense_dask", "sparse_dask"] +) +def test_normalize_clr_dask(sparse_blocks, alpha): + 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, alpha=alpha) + + result = _materialize(adata.X) + np.testing.assert_allclose( + result, _clr_reference(X_clr, alpha=alpha), rtol=1e-5, atol=1e-5 + ) + assert isinstance(adata.X, DaskArray) + assert isinstance(adata.X._meta, np.ndarray) + + +def test_normalize_clr_view(): + adata = AnnData(X_clr.copy()) + v = adata[:, :] + with pytest.warns(UserWarning, match=r"Received a view"): + sc.pp.normalize_clr(v) + assert not v.is_view diff --git a/tests/test_package_structure.py b/tests/test_package_structure.py index 6b2b3a69b8..a8d6b07e4f 100644 --- a/tests/test_package_structure.py +++ b/tests/test_package_structure.py @@ -102,6 +102,9 @@ class ExpectedSig(TypedDict): copy_sigs["sc.pp.scale"]["first_name"] = "data" copy_sigs["sc.pp.sqrt"]["first_name"] = "data" # other partial exceptions +copy_sigs["sc.pp.normalize_clr"]["return_ann"] = ( + "AnnData | dict[str, np.ndarray | DaskArray] | None" +) copy_sigs["sc.pp.normalize_total"]["return_ann"] = copy_sigs[ "sc.experimental.pp.normalize_pearson_residuals" ]["return_ann"] = "AnnData | dict[str, np.ndarray] | None"