Skip to content

refactor: completely deprecated bare fold_change column #80

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
noamteyssier merged 1 commit into from
Jul 14, 2026
Merged

refactor: completely deprecated bare fold_change column #80

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 2 additions & 3 deletions CLAUDE.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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()`
Expand All @@ -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) |
Expand Down
23 changes: 4 additions & 19 deletions src/pdex/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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**.
Expand All @@ -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.

Expand All @@ -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)

Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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,
Expand All @@ -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,
Expand Down
27 changes: 13 additions & 14 deletions tests/test_pdex.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,6 @@
"ref_mean",
"target_membership",
"ref_membership",
"fold_change",
"log2_fold_change",
"percent_change",
"p_value",
Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -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(),
Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -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(),
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -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(),
Expand Down Expand Up @@ -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.
Expand All @@ -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(),
Expand Down Expand Up @@ -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(),
Expand All @@ -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(),
Expand Down Expand Up @@ -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)
Expand All @@ -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
Expand Down Expand Up @@ -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",
Expand Down
Loading