From ba89e992fed90600e0c55fb0ee6fea57cc71e762 Mon Sep 17 00:00:00 2001 From: Erik Date: Tue, 22 Sep 2026 15:14:31 +0200 Subject: [PATCH] fix: shift signed CLR only for differential_abundance Wilcoxon Keep mean difference on the original matrix, like PixelatorR RunDAA. Co-authored-by: Cursor --- CHANGELOG.md | 3 + .../pna/analysis/_differential_abundance.py | 97 +++++++++++-------- .../analysis/test_differential_abundance.py | 11 ++- 3 files changed, 69 insertions(+), 42 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 166f75457..0ff637ab4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -24,6 +24,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - `pixelator.pna.analysis.filter_proximity_scores` to filter a proximity score table by marker abundance (Python analog of pixelatorR `FilterProximityScores`). ### Changed +- `differential_abundance` runs Wilcoxon on a per-marker shift to non-negative + values when the selected matrix contains negatives, and still reports mean + difference on the original matrix. - Updated pixelgen-pixelator-core to 0.2.0 improving peak memory usage in the graph step by ~20%. - `density_scatter_plot` now lives in `pixelator.plot` (previously `pixelator.mpx.plot`). - `uei_count` is now optional on PNA edgelists in `sample_calling` and the graph component diff --git a/src/pixelator/pna/analysis/_differential_abundance.py b/src/pixelator/pna/analysis/_differential_abundance.py index cf68efe0b..752b0a924 100644 --- a/src/pixelator/pna/analysis/_differential_abundance.py +++ b/src/pixelator/pna/analysis/_differential_abundance.py @@ -5,7 +5,6 @@ from __future__ import annotations -import warnings from collections.abc import Iterable, Sequence from typing import Literal @@ -36,6 +35,7 @@ } _SCANPY_KEY = "_pixelator_rank_genes_groups" +_NONNEG_TEST_LAYER = "_pixelator_daa_nonneg" _RESULT_COLUMNS = [ "marker", "p", @@ -90,14 +90,28 @@ def differential_abundance( the same as ``RunDAA``. A further **global FDR** across separately computed cell-type tables (PAT 05 notebook glue) is **not** applied here. - Scanpy's Wilcoxon historically assumes non-negative values. Pixelator's - default CLR (``clr_transformation(..., non_negative=True)``) is - non-negative and is safe to pass via ``layer``. Signed CLR (negatives) does - not crash the Wilcoxon call in current scanpy, but percent-expressed - (values ``> 0``) is uninformative and log-fold changes would be invalid; - this helper therefore reports mean difference on the original matrix. - For signed CLR, pass counts in ``X``, a non-negative CLR layer (the - pixelator default), or shift each marker so its minimum is 0. + For each target vs reference (and each ``group_vars`` stratum) the + procedure is: + + 1. Restrict to those cells and to ``features`` if given. + 2. On that original matrix, compute ``difference`` as + ``mean(target) - mean(reference)`` and ``pct_1`` / ``pct_2`` as the + fraction of cells with value ``> 0``. This is the effect size + ``RunDAA`` reports with Seurat ``mean.fxn = rowMeans`` and + ``fc.name = "difference"``, not scanpy's log-fold change. + 3. If any value is negative (signed CLR), build a **test-only** copy: + for each marker, subtract its minimum when that minimum is negative, + so every marker is non-negative. The same additive shift is applied + to every cell, so ranks between groups are unchanged. + 4. Call ``scanpy.tl.rank_genes_groups(..., method="wilcoxon")`` on that + non-negative matrix (the original matrix if it was already + non-negative). Keep the raw p-values; discard scanpy's log-fold + changes and its per-call p-adjustment. + + Pixelator's default CLR (``clr_transformation(..., non_negative=True)``) + is already non-negative, so step 3 is a no-op. Percent-expressed on + signed CLR remains the fraction of cells with value ``> 0`` on the + original matrix, which is not a count-based detection rate. Args: adata: AnnData of components × markers. ``contrast_column`` and any @@ -137,7 +151,8 @@ def differential_abundance( * ``p_adj`` — adjusted p-value (see above). Non-finite raw p-values are left as NaN and excluded from the adjustment, like R ``p.adjust``. - * ``difference`` — mean(target) − mean(reference) on the selected matrix + * ``difference`` — mean(target) − mean(reference) on the original + selected matrix (not on the non-negative test copy) * ``pct_1`` — fraction of target cells with value ``> 0`` (scanpy ``pts`` / Seurat ``pct.1``) * ``pct_2`` — fraction of reference cells with value ``> 0`` (Seurat @@ -344,39 +359,34 @@ def _run_one_scanpy_wilcoxon( if features is not None: work = work[:, features].copy() matrix = _values_matrix(work, layer) + effects = _mean_difference_and_pct(work, contrast_column, target, reference, layer) has_negatives = _has_negatives(matrix) - if not warned_negatives and has_negatives: - logger.warning( - "Selected matrix contains negative values. scanpy Wilcoxon p-values " - "are still computed, but percent-expressed (value > 0) is not " - "meaningful on signed CLR. Mean difference is reported on the " - "original values. Use counts, or non-negative CLR " - "(clr_transformation(..., non_negative=True))." - ) - warned_negatives = True - - # Scanpy still computes log-fold changes internally; they are unused here - # (we report mean difference) and are NaN on signed CLR. - with warnings.catch_warnings(): - if has_negatives: - warnings.filterwarnings( - "ignore", - message="invalid value encountered in log2", - category=RuntimeWarning, + scanpy_layer: str | None + if has_negatives: + if not warned_negatives: + logger.info( + "Selected matrix contains negative values. Wilcoxon is run on a " + "per-marker shift to non-negative values; mean difference is " + "still computed on the original matrix." ) - sc.tl.rank_genes_groups( - work, - groupby=contrast_column, - groups=[target], - reference=reference, - method="wilcoxon", - use_raw=False, - layer=layer, - n_genes=work.n_vars, - key_added=_SCANPY_KEY, - ) + warned_negatives = True + work.layers[_NONNEG_TEST_LAYER] = _shift_markers_to_nonnegative(matrix) + scanpy_layer = _NONNEG_TEST_LAYER + else: + scanpy_layer = layer + + sc.tl.rank_genes_groups( + work, + groupby=contrast_column, + groups=[target], + reference=reference, + method="wilcoxon", + use_raw=False, + layer=scanpy_layer, + n_genes=work.n_vars, + key_added=_SCANPY_KEY, + ) ranked = sc.get.rank_genes_groups_df(work, group=target, key=_SCANPY_KEY) - effects = _mean_difference_and_pct(work, contrast_column, target, reference, layer) result = pd.DataFrame( { "marker": ranked["names"].astype(str), @@ -439,6 +449,13 @@ def _has_negatives(matrix) -> bool: return bool(minimum < 0) +def _shift_markers_to_nonnegative(matrix) -> np.ndarray: + """Shift each marker so its minimum is at least 0, for the Wilcoxon test only.""" + dense = matrix.toarray() if issparse(matrix) else np.asarray(matrix) + mins = np.nanmin(dense, axis=0) + return dense - np.minimum(mins, 0) + + def _mean_difference_and_pct( adata: AnnData, contrast_column: str, diff --git a/tests/pna/analysis/test_differential_abundance.py b/tests/pna/analysis/test_differential_abundance.py index 74fc216de..6ac0e2d14 100644 --- a/tests/pna/analysis/test_differential_abundance.py +++ b/tests/pna/analysis/test_differential_abundance.py @@ -184,8 +184,11 @@ def test_differential_abundance_numpy_obsm_with_features(): assert len(result) == 2 -def test_differential_abundance_clr_like_negatives_do_not_crash(): +def test_differential_abundance_signed_clr_keeps_original_difference_sign(): adata = _tiny_adata(with_negatives=True) + original = np.asarray(adata.X).copy() + assert original.min() < 0 + result = differential_abundance( adata, contrast_column="condition", @@ -194,6 +197,10 @@ def test_differential_abundance_clr_like_negatives_do_not_crash(): ) assert list(result.columns) == EXPECTED_COLUMNS - assert result["difference"].notna().all() + np.testing.assert_allclose(np.asarray(adata.X), original) planted = result.set_index("marker").loc["M0"] + treated = original[adata.obs["condition"].eq("treated").to_numpy(), 0].mean() + control = original[adata.obs["condition"].eq("control").to_numpy(), 0].mean() + assert planted["difference"] == pytest.approx(treated - control) assert planted["difference"] > 0 + assert np.isfinite(planted["p"])