feat(preprocessing): add optimal shifted CLR (PFlog1pPF) normalization

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

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},

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`

docs/release-notes/4160.feat.md

@@ -0,0 +1 @@
+Add {func}`scanpy.pp.normalize_clr` for PFlog shifted centered log-ratio normalization, a variance-stabilizing, depth-invariant and rank-preserving count transform {cite:p}`Booeshaghi2022`. The default representation stores sparse shifted-log values together with per-cell centering offsets, and {func}`scanpy.pp.pca` can use this representation directly for memory-efficient PCA without materializing the dense centered matrix. {smaller}`R Baber`

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",

src/scanpy/preprocessing/_normalization.py

@@ -7,13 +7,17 @@
 import numpy as np
 from fast_array_utils import stats
 from fast_array_utils.numba import njit
+from fast_array_utils.stats import mean_var
+from scipy import sparse
 
 from .. import logging as logg
 from .._compat import CSBase, CSCBase, CSRBase, DaskArray, warn
 from .._utils import axis_mul_or_truediv, dematrix, view_to_actual
 from ..get import _get_obs_rep, _set_obs_rep
 
 if TYPE_CHECKING:
+    from typing import Literal
+
     from anndata import AnnData
 
 
@@ -304,3 +308,281 @@ def normalize_total(  # noqa: PLR0912
     elif not inplace:
         return dat
     return None
+
+
+def _estimate_overdispersion(x: np.ndarray | CSBase | DaskArray) -> float:
+    r"""Estimate the negative-binomial overdispersion :math:`α` from raw counts.
+
+    Fits :math:`\mathrm{Var}_g = μ_g + α \cdot μ_g^2` across genes, where
+    :math:`μ_g` and :math:`\mathrm{Var}_g` are the per-gene mean and (population)
+    variance over cells. The model is linear in :math:`α`, so the ordinary
+    least-squares solution is closed form
+
+    .. math::
+        α = \frac{\sum_g (\mathrm{Var}_g - μ_g) \, μ_g^2}{\sum_g μ_g^4},
+
+    which is exactly the minimizer a non-linear `curve_fit` would converge to,
+    but without the dependency. :func:`~fast_array_utils.stats.mean_var` is
+    dispatched for dense, sparse and dask input alike, so the only dask-specific
+    step is computing the two final scalar sums.
+    """
+    mu, var = mean_var(x, axis=0, correction=0)
+    mu2 = mu**2
+    numerator = np.sum((var - mu) * mu2)
+    denominator = np.sum(mu2 * mu2)
+    if isinstance(x, DaskArray):
+        import dask
+
+        numerator, denominator = dask.compute(numerator, denominator)
+    if denominator == 0.0:
+        msg = (
+            "Cannot estimate overdispersion: every gene has zero mean. "
+            "Pass a positive `alpha` explicitly."
+        )
+        raise ValueError(msg)
+    alpha = float(numerator / denominator)
+    if not alpha > 0:
+        msg = (
+            f"Estimated overdispersion is non-positive (alpha = {alpha}); "
+            "pass a positive `alpha` explicitly."
+        )
+        raise ValueError(msg)
+    return alpha
+
+
+def _log1p_sparse_block(x: np.ndarray | CSBase) -> np.ndarray | CSBase:
+    """Apply log1p to a dense or sparse block while preserving sparse zeros."""
+    if isinstance(x, CSBase):
+        x = x.copy()
+        x.data = np.log1p(x.data)
+        return x
+    return np.log1p(x)
+
+
+def _densify_shifted_clr(
+    log_values: np.ndarray | CSBase | DaskArray, row_center: np.ndarray | DaskArray
+) -> np.ndarray | DaskArray:
+    dense_log = log_values.toarray() if isinstance(log_values, CSBase) else log_values
+    return dense_log - row_center[:, None]
+
+
+def _normalize_clr_helper(  # noqa: PLR0912
+    x: np.ndarray | CSBase | DaskArray,
+    *,
+    target: Literal["auto", "mean", "median"] | float,
+    alpha: float | None,
+) -> tuple[
+    np.ndarray | CSBase | DaskArray,
+    np.ndarray | DaskArray,
+    np.ndarray | DaskArray,
+    dict[str, object],
+]:
+    """Compute sparse PFlog / shifted-CLR values and row centers."""
+    # Keep the depths lazy for dask; `.ravel()` would otherwise materialize them.
+    cell_depths = stats.sum(x, axis=1)
+    if not isinstance(x, DaskArray):
+        cell_depths = np.asarray(cell_depths).ravel()
+
+    mean_depth = (
+        float(cell_depths.mean().compute())
+        if isinstance(cell_depths, DaskArray)
+        else float(cell_depths.mean())
+    )
+    if alpha is not None:
+        if not alpha > 0:
+            msg = (
+                f"`alpha` must be positive to compute PFlog, got {alpha}. "
+                "The data may be underdispersed."
+            )
+            raise ValueError(msg)
+        scale: float | np.ndarray | DaskArray = 4.0 * float(alpha)
+        target_sum = 4.0 * float(alpha) * mean_depth
+        resolved_target = "alpha"
+    elif target == "auto":
+        alpha = _estimate_overdispersion(x)
+        scale = 4.0 * alpha
+        target_sum = 4.0 * alpha * mean_depth
+        resolved_target = "auto"
+    else:
+        alpha = None
+        if target == "mean":
+            target_sum = mean_depth
+        elif target == "median":
+            if isinstance(cell_depths, DaskArray):
+                msg = (
+                    "`target='median'` is not supported for dask input; pass "
+                    "`target='auto'`, `target='mean'`, a positive numeric target, "
+                    "or explicit `alpha`."
+                )
+                raise NotImplementedError(msg)
+            target_sum = float(np.median(cell_depths))
+        elif isinstance(target, bool) or not isinstance(target, int | float):
+            msg = "`target` must be 'auto', 'mean', 'median', or a positive number."
+            raise ValueError(msg)
+        else:
+            target_sum = float(target)
+        if not target_sum > 0:
+            msg = f"`target` must resolve to a positive depth, got {target_sum}."
+            raise ValueError(msg)
+        resolved_target = str(target)
+        x = axis_mul_or_truediv(
+            x, cell_depths / target_sum, op=truediv, axis=0, allow_divide_by_zero=False
+        )
+
+    if alpha is not None:
+        x = x * scale
+
+    if isinstance(x, DaskArray):
+        log_values = x.map_blocks(
+            _log1p_sparse_block, dtype=np.float64, meta=x._meta.astype(np.float64)
+        )
+    else:
+        log_values = _log1p_sparse_block(x)
+
+    row_center = stats.sum(log_values, axis=1) / x.shape[1]
+    if not isinstance(row_center, DaskArray):
+        row_center = np.asarray(row_center).ravel()
+    report = dict(
+        target=resolved_target,
+        alpha=None if alpha is None else float(alpha),
+        k=float(target_sum),
+        mean_depth=mean_depth,
+    )
+    return log_values, row_center, cell_depths, report
+
+
+def normalize_clr(
+    adata: AnnData,
+    *,
+    target: Literal["auto", "mean", "median"] | float = "auto",
+    alpha: float | None = None,
+    key_added: str = "pflog",
+    densify: bool = False,
+    layer: str | None = None,
+    inplace: bool = True,
+    copy: bool = False,
+) -> AnnData | dict[str, np.ndarray | CSBase | DaskArray | dict[str, object]] | None:
+    r"""Normalize counts with the shifted centered log-ratio (PFlog) transform.
+
+    With `target="auto"` (or explicit `alpha`), computes PFlog as
+
+    .. math::
+        T(x)_i = \log(1 + 4 α x_i)
+            - \frac{1}{D} \sum_{j=1}^D \log(1 + 4 α x_j),
+
+    which is equivalent to centering
+    :math:`\log(x_i + 1 / (4 α))` because the constant :math:`\log(4 α)`
+    cancels during CLR centering. For memory efficiency, the sparse log values
+    are stored separately from the per-cell centering vector.
+
+    .. note::
+        `target="auto"` is the recommended PFlog default. Fixed targets
+        (``"mean"``, ``"median"``, or a positive number) apply shifted CLR
+        after rescaling each cell to the same total count.
+
+    Parameters
+    ----------
+    adata
+        The annotated data matrix of shape `n_obs` × `n_vars`.
+        Rows correspond to cells and columns to genes.
+    target
+        ``"auto"`` estimates the negative-binomial overdispersion and applies
+        PFlog. ``"mean"``, ``"median"``, or a positive numeric value set a
+        fixed total count used before applying the shifted CLR transform.
+    alpha
+        Negative-binomial overdispersion of the dataset (``var = μ + α·μ²``).
+        When given, it overrides `target` and applies PFlog with sparse values
+        ``log1p(4 * alpha * x)``.
+    key_added
+        Layer key used to store the sparse log values. The row-centering vector
+        is stored in ``adata.obs[f"{key_added}_center"]``.
+    densify
+        If `True`, also materialize the dense centered matrix into `X` or the
+        selected `layer`. This is intended for small data or generic downstream
+        tools that cannot consume the sparse-plus-row-center representation.
+    layer
+        Layer to normalize instead of `X`.
+    inplace
+        Whether to update `adata` or return a dictionary with the normalized
+        matrix.
+    copy
+        Whether to modify a copied input object. Not compatible with
+        `inplace=False`.
+
+    Returns
+    -------
+    Returns a dictionary with normalized values and metadata or updates `adata`,
+    depending on `inplace`.
+
+    Example
+    -------
+    >>> import numpy as np
+    >>> from anndata import AnnData
+    >>> import scanpy as sc
+    >>> adata = AnnData(np.array([[1, 2, 30], [4, 50, 6]], dtype="float32"))
+    >>> sc.pp.normalize_clr(adata, alpha=0.5)
+    >>> "pflog" in adata.layers and "pflog_center" in adata.obs
+    True
+    """
+    if copy:
+        if not inplace:
+            msg = "`copy=True` cannot be used with `inplace=False`."
+            raise ValueError(msg)
+        adata = adata.copy()
+
+    view_to_actual(adata)
+
+    x = _get_obs_rep(adata, layer=layer)
+    if isinstance(x, CSCBase):
+        x = x.tocsr()
+    if not inplace:
+        x = x.copy()
+    if issubclass(x.dtype.type, int | np.integer):
+        x = x.astype(np.float64)
+
+    start = logg.info("normalizing counts per cell via PFlog")
+
+    log_values, row_center, cell_depths, report = _normalize_clr_helper(
+        x, target=target, alpha=alpha
+    )
+
+    if not isinstance(cell_depths, DaskArray) and not np.all(cell_depths > 0):
+        warn("Some cells have zero counts", UserWarning)
+
+    row_center_obs = (
+        np.asarray(row_center.compute()).ravel()
+        if isinstance(row_center, DaskArray)
+        else np.asarray(row_center).ravel()
+    )
+    if not isinstance(log_values, DaskArray | CSBase):
+        log_values = sparse.csr_matrix(log_values)  # noqa: TID251
+
+    dense_x = _densify_shifted_clr(log_values, row_center) if densify else None
+    row_center_key = f"{key_added}_center"
+    metadata = dict(
+        encoding_type="shifted_clr",
+        row_center_key=row_center_key,
+        params=report | {"layer": layer, "densify": densify},
+    )
+    dat = dict(
+        X=dense_x if densify else log_values,
+        row_center=row_center_obs,
+        metadata=metadata,
+    )
+    if inplace:
+        adata.layers[key_added] = log_values
+        adata.obs[row_center_key] = row_center_obs
+        adata.uns[key_added] = metadata
+        if densify:
+            _set_obs_rep(adata, dense_x, layer=layer)
+
+    logg.info(
+        "    finished ({time_passed})",
+        time=start,
+    )
+
+    if copy:
+        return adata
+    elif not inplace:
+        return dat
+    return None