diff --git a/CLAUDE.md b/CLAUDE.md index 47aa6ab..b1004f6 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -73,12 +73,12 @@ checked empirically (inspect the per-gene CPM distribution). | File | Role | | ---------------------- | ------------------------------------------------------------------------------------------------------- | | `src/pdex/__init__.py` | `pdex()` entry point and full pipeline logic | -| `src/pdex/_math.py` | Numba JIT-compiled `fold_change()`, `percent_change()`, and `mwu()`/`mwu_one_vs_rest()` wrappers over `numba-mwu`; `pseudobulk()` dispatcher; `cpm_bulk()` pooled-CPM view for the filter; the `bulk_matrix_pre_transform_mean()`/`pseudobulk_from_pre_mean()`/`cpm_from_gene_means()` trio powering `"all"` mode's global-sum optimization | +| `src/pdex/_math.py` | Numba JIT-compiled `log2_fold_change()`, `percent_change()`, and `mwu()`/`mwu_one_vs_rest()` wrappers over `numba-mwu`; `pseudobulk()` dispatcher; `cpm_bulk()` pooled-CPM view for the filter; the `bulk_matrix_pre_transform_mean()`/`pseudobulk_from_pre_mean()`/`cpm_from_gene_means()` trio powering `"all"` mode's global-sum optimization | | `src/pdex/_utils.py` | `set_numba_threadpool()` — sets Numba thread count before JIT warmup; `_available_cpus()` — affinity-aware CPU count (respects cgroup/SLURM limits); `_detect_is_log1p()` heuristic | ### Performance Design -- Numba JIT compilation accelerates per-cell/per-gene math (`fold_change`, `percent_change`, `_log1p_col_mean`, `_expm1_vec`) +- Numba JIT compilation accelerates per-cell/per-gene math (`log2_fold_change`, `percent_change`, `_log1p_col_mean`, `_expm1_vec`) - `numba-mwu` (external dep, `>=0.2.0`) provides Numba-accelerated Mann-Whitney U kernels for **both** the pairwise case (`mannwhitneyu_columns`/`mannwhitneyu_sparse`, used by `"ref"` and `"on_target"` modes) and the one-vs-rest case (`mannwhitneyu_one_vs_rest`/`_sparse`, used by `"all"` mode) — the one-vs-rest kernels originated in pdex and were upstreamed into `numba-mwu` since the optimization is domain-agnostic (see that package's CLAUDE.md for the algorithm). - Sparse CSR matrices are handled by reusing pre-computed non-targeting column indices to avoid redundant dense conversion - Parallelism is controlled via `threads` passed to `set_numba_threadpool()` @@ -99,7 +99,6 @@ The returned Polars DataFrame (or pandas DataFrame when `as_pandas=True`) has co | `ref_mean` | float | Pseudobulk mean for the reference, always in natural (count) space | | `target_membership` | int | Number of cells in the target group | | `ref_membership` | int | Number of cells in the reference | -| `fold_change` | float | **Deprecated** alias for `log2_fold_change` (identical values). Retained for one release; emits a `FutureWarning` on every `pdex(...)` call and will be removed in pdex 0.3.0. | | `log2_fold_change` | float | log2((target_mean + epsilon) / (ref_mean + epsilon)) — computed from pseudobulk means. `epsilon` defaults to `1e-9` (finite-guard), so by default there are no `±inf`/`NaN`: one-sided zeros become large-but-finite and `0/0` is `0.0`. With `epsilon == 0`, `0/0` is still defined as `0.0` (not `NaN`) but one-sided zeros yield `±inf`. | | `percent_change` | float | (target_mean - ref_mean) / (ref_mean + epsilon) — computed from pseudobulk means. With the default `epsilon=1e-9` there are no non-finite values; with `epsilon == 0`, `0/0` is `0.0` (not `NaN`) and a zero reference with nonzero target yields `+inf`. | | `p_value` | float | Mann-Whitney U p-value (per-cell vectors) | diff --git a/src/pdex/__init__.py b/src/pdex/__init__.py index 290cf4f..e8c0d8e 100644 --- a/src/pdex/__init__.py +++ b/src/pdex/__init__.py @@ -256,7 +256,7 @@ def pdex( epsilon: Pseudocount added to the **native (count-space) means** — the denominator (and, for ``log2_fold_change``, the numerator) — before computing - ``fold_change`` and ``percent_change``. It is never applied to the CPM view + ``log2_fold_change`` and ``percent_change``. It is never applied to the CPM view used by ``cpm_filter``. When ``epsilon > 0``, extreme values from near-zero reference means (scRNA-seq sparsity artifact) are dampened toward zero, and one-sided zeros become large-but-finite instead of ``±inf``. Has no effect on @@ -307,7 +307,7 @@ def pdex( pl.DataFrame | pd.DataFrame One row per (group, feature) pair with columns: ``target``, ``feature``, ``target_mean``, ``ref_mean``, ``target_membership``, ``ref_membership``, - ``fold_change``, ``log2_fold_change``, ``percent_change``, ``p_value``, + ``log2_fold_change``, ``percent_change``, ``p_value``, ``statistic``, ``fdr``. ``target_mean`` and ``ref_mean`` are always in **natural (count) space**. @@ -322,12 +322,8 @@ def pdex( and ``percent_change`` define this as ``0.0`` (no change) rather than ``NaN``. One-sided zeros still produce ``±inf``. - ``fold_change`` is a **deprecated** alias for ``log2_fold_change`` - (identical values). It is retained for one release to ease migration - and will be removed in pdex 0.3.0. New code should read - ``log2_fold_change`` directly. A :class:`FutureWarning` is emitted - on every ``pdex(...)`` call. The MWU ``p_value`` and - ``statistic`` are computed directly on the per-cell expression vectors. + The MWU ``p_value`` and ``statistic`` are computed directly on the + per-cell expression vectors. For ``mode="ref"``, the reference group itself is excluded from the output. @@ -345,14 +341,6 @@ def pdex( if epsilon < 0: raise ValueError(f"epsilon must be non-negative, got {epsilon}") - warnings.warn( - "The `fold_change` column in pdex output is deprecated and will be " - "removed in pdex 0.3.0. Use `log2_fold_change` instead — it contains " - "the same values (`log2(target_mean / ref_mean)`).", - FutureWarning, - stacklevel=2, - ) - # Set the global threadpool for numba set_numba_threadpool(threads) @@ -502,7 +490,6 @@ def _assemble_group_frame( "ref_mean": ref_mean, "target_membership": np.full(feature.shape[0], target_membership), "ref_membership": np.full(feature.shape[0], ref_membership), - "fold_change": lfc, "log2_fold_change": lfc, "percent_change": pc, "p_value": pvalue, @@ -790,7 +777,6 @@ def _pdex_on_target( "ref_mean": ref_mean, "target_membership": group_mask.size, "ref_membership": ref_membership, - "fold_change": lfc, "log2_fold_change": lfc, "percent_change": pc, "p_value": p_value, @@ -809,7 +795,6 @@ def _pdex_on_target( "ref_mean": pl.Float64, "target_membership": pl.Int64, "ref_membership": pl.Int64, - "fold_change": pl.Float64, "log2_fold_change": pl.Float64, "percent_change": pl.Float64, "p_value": pl.Float64, diff --git a/tests/test_pdex.py b/tests/test_pdex.py index 177d974..a5f38b5 100644 --- a/tests/test_pdex.py +++ b/tests/test_pdex.py @@ -17,7 +17,6 @@ "ref_mean", "target_membership", "ref_membership", - "fold_change", "log2_fold_change", "percent_change", "p_value", @@ -87,7 +86,7 @@ def test_fold_change_sign(self, small_adata): for group_name in ["A", "B"]: group_rows = result.filter(pl.col("target") == group_name) # Mean fold change should be positive since we boosted these groups - mean_fc = group_rows["fold_change"].mean() + mean_fc = group_rows["log2_fold_change"].mean() assert mean_fc > 0, f"Expected positive fold change for group {group_name}" # type: ignore def test_statistics_against_scipy(self, small_adata): @@ -177,7 +176,7 @@ def test_sparse_dense_agreement(self, small_adata, small_adata_sparse): assert dense_result.shape == sparse_result.shape - for col in ["p_value", "statistic", "fold_change", "percent_change"]: + for col in ["p_value", "statistic", "log2_fold_change", "percent_change"]: np.testing.assert_allclose( dense_result[col].to_numpy(), sparse_result[col].to_numpy(), @@ -237,7 +236,7 @@ def test_fold_change_sign(self, small_adata): """Group B was boosted the most, so its fold change vs rest should be positive.""" result = pdex(small_adata, groupby="guide", mode="all", is_log1p=False) group_b_rows = result.filter(pl.col("target") == "B") - mean_fc = group_b_rows["fold_change"].mean() + mean_fc = group_b_rows["log2_fold_change"].mean() assert mean_fc > 0 # type: ignore def test_statistics_against_scipy(self, small_adata): @@ -277,7 +276,7 @@ def test_sparse_dense_agreement(self, small_adata, small_adata_sparse): assert dense_result.shape == sparse_result.shape - for col in ["p_value", "statistic", "fold_change", "percent_change"]: + for col in ["p_value", "statistic", "log2_fold_change", "percent_change"]: np.testing.assert_allclose( dense_result[col].to_numpy(), sparse_result[col].to_numpy(), @@ -351,7 +350,7 @@ def test_sparse_dense_agreement(self, multi_group_adata, multi_group_adata_spars for col in [ "p_value", "statistic", - "fold_change", + "log2_fold_change", "percent_change", "target_mean", "ref_mean", @@ -572,7 +571,7 @@ def test_sparse_dense_agreement(self, on_target_adata, on_target_adata_sparse): ) assert dense_result.shape == sparse_result.shape - for col in ["p_value", "statistic", "fold_change", "percent_change"]: + for col in ["p_value", "statistic", "log2_fold_change", "percent_change"]: np.testing.assert_allclose( dense_result[col].to_numpy(), sparse_result[col].to_numpy(), @@ -744,7 +743,7 @@ def test_both_log1p_paths_agree(self, small_adata, small_adata_log1p): """pdex on raw counts with is_log1p=False and on log1p counts with is_log1p=True should yield identical results across all output columns. - Pseudobulk means back-transform to the same count space, so fold_change and + Pseudobulk means back-transform to the same count space, so log2_fold_change and percent_change must match. The MWU statistic and p_value operate on the raw cell-level values (which differ between the two inputs), so they are NOT expected to match — only the pseudobulk-derived columns are tested here. @@ -763,7 +762,7 @@ def test_both_log1p_paths_agree(self, small_adata, small_adata_log1p): is_log1p=True, geometric_mean=True, ) - for col in ["target_mean", "ref_mean", "fold_change", "percent_change"]: + for col in ["target_mean", "ref_mean", "log2_fold_change", "percent_change"]: np.testing.assert_allclose( raw_result[col].to_numpy(), log_result[col].to_numpy(), @@ -797,7 +796,7 @@ def test_ref_mode_backed_matches_inmemory(self, small_adata, small_adata_backed) inmem = pdex(small_adata, groupby="guide", mode="ref", is_log1p=False) backed = pdex(small_adata_backed, groupby="guide", mode="ref", is_log1p=False) assert inmem.shape == backed.shape - for col in ["p_value", "statistic", "fold_change", "percent_change"]: + for col in ["p_value", "statistic", "log2_fold_change", "percent_change"]: np.testing.assert_allclose( inmem[col].to_numpy(), backed[col].to_numpy(), @@ -809,7 +808,7 @@ def test_all_mode_backed_matches_inmemory(self, small_adata, small_adata_backed) inmem = pdex(small_adata, groupby="guide", mode="all", is_log1p=False) backed = pdex(small_adata_backed, groupby="guide", mode="all", is_log1p=False) assert inmem.shape == backed.shape - for col in ["p_value", "statistic", "fold_change", "percent_change"]: + for col in ["p_value", "statistic", "log2_fold_change", "percent_change"]: np.testing.assert_allclose( inmem[col].to_numpy(), backed[col].to_numpy(), @@ -849,7 +848,7 @@ def test_zero_in_both_is_zero_not_nan(self, small_adata, mode): assert (gene0["target_mean"].to_numpy() == 0).all() assert (gene0["ref_mean"].to_numpy() == 0).all() - for col in ["log2_fold_change", "fold_change", "percent_change"]: + for col in ["log2_fold_change", "percent_change"]: values = gene0[col].to_numpy() assert not np.isnan(values).any(), f"{col} contains NaN" np.testing.assert_array_equal(values, 0.0) @@ -870,7 +869,7 @@ def test_on_target_zero_in_both_is_zero_not_nan(self, on_target_adata): row = result.filter(pl.col("target") == "A") assert row["target_mean"].to_numpy()[0] == 0 assert row["ref_mean"].to_numpy()[0] == 0 - for col in ["log2_fold_change", "fold_change", "percent_change"]: + for col in ["log2_fold_change", "percent_change"]: value = row[col].to_numpy()[0] assert not np.isnan(value), f"{col} is NaN" assert value == 0.0 @@ -1129,7 +1128,7 @@ def test_backed_sparse_matches_inmemory_sparse( for col in [ "target_mean", "ref_mean", - "fold_change", + "log2_fold_change", "percent_change", "p_value", "statistic",