From df3b8fb633eed9a179943e37d0be73af3ca45f4b Mon Sep 17 00:00:00 2001 From: Leon Hafner Date: Fri, 26 Jun 2026 00:11:27 +0000 Subject: [PATCH 1/3] add CPM filter --- CLAUDE.md | 30 ++++- pyproject.toml | 2 +- src/pdex/__init__.py | 295 ++++++++++++++++++++++++++++++++++++------- src/pdex/_math.py | 24 ++++ tests/conftest.py | 40 ++++++ tests/test_math.py | 49 ++++++- tests/test_pdex.py | 225 ++++++++++++++++++++++++++++++++- 7 files changed, 605 insertions(+), 60 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index b6e8b0a..89ac38c 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -36,13 +36,13 @@ uv run ty check ### Core Pipeline (`src/pdex/__init__.py`) -The main entry point is `pdex(adata, groupby, mode, threads, is_log1p, geometric_mean, as_pandas, epsilon, **kwargs)`, which: +The main entry point is `pdex(adata, groupby, mode, threads, is_log1p, geometric_mean, as_pandas, epsilon, cpm_filter, **kwargs)`, which: 1. Validates the `groupby` column in `adata.obs` 2. Extracts unique groups (filters NaN and empty strings) 3. Identifies a reference group (defaults to `"non-targeting"` in `"ref"` and `"on_target"` modes) 4. For each non-reference group, slices the expression matrix, computes pseudobulk (mean), fold change, percent change, and Mann-Whitney U statistic vs the reference -5. Applies per-group FDR correction (scipy) and returns a Polars DataFrame (or pandas if `as_pandas=True`) +5. Optionally drops genes below the `cpm_filter` floor (see below), then applies FDR correction over the surviving genes (scipy) and returns a Polars DataFrame (or pandas if `as_pandas=True`) Three modes: @@ -52,12 +52,28 @@ Three modes: Unexpected `**kwargs` for any mode trigger a `UserWarning`. +### CPM floor filter (`cpm_filter`) + +`cpm_filter` (default `None` = off) is an opt-in **two-view** floor filter. A native +(unnormalized) view drives the reported means, LFC, MWU, and output; a separate **CPM +view** drives only the keep/drop decision, so the output is never normalized. Per group +per gene the pooled (bulk) CPM is `Σcounts_gene / Σcounts_all * 1e6` (computed on counts — +`expm1` is applied first when `is_log1p`). A `(target, gene)` row is **dropped** when the +gene's CPM is `<= T` in **both** the target and the reference (kept iff `target_cpm > T` +**or** `ref_cpm > T`; strict `>`, negative `T` keeps everything). The CPM ratio is +scale-invariant, so `T` means the same regardless of input normalization. The drop is +independent of the MWU result, and **FDR is corrected over the surviving genes only**. +Applies to all three modes (in `on_target`, a group whose target gene is a floor gene is +dropped). Emits a `UserWarning` if the data contains negative values. `T = 5` is a +reasonable starting point, but the optimal threshold is dataset-dependent and should be +checked empirically (inspect the per-gene CPM distribution). + ### Key Files | 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()` wrappers; `pseudobulk()` dispatcher | +| `src/pdex/_math.py` | Numba JIT-compiled `fold_change()`, `percent_change()`, and `mwu()` wrappers; `pseudobulk()` dispatcher; `cpm_bulk()` pooled-CPM view for the filter | | `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 @@ -80,14 +96,16 @@ The returned Polars DataFrame (or pandas DataFrame when `as_pandas=True`) has co | `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. Features unexpressed in both groups (`target_mean == ref_mean == 0`, only with `epsilon == 0`) give `0/0`, defined as `0.0` (not `NaN`); one-sided zeros still yield `±inf`. | -| `percent_change` | float | (target_mean - ref_mean) / (ref_mean + epsilon) — computed from pseudobulk means. Features unexpressed in both groups (`target_mean == ref_mean == 0`, only with `epsilon == 0`) give `0/0`, defined as `0.0` (not `NaN`); a zero reference with nonzero target still yields `+inf`. | +| `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) | | `statistic` | float | Mann-Whitney U statistic | -| `fdr` | float | FDR-corrected p-value, applied per-group across genes. For `on_target` mode, applied across all groups. | +| `fdr` | float | FDR-corrected p-value, applied per-group across genes. For `on_target` mode, applied across all groups. When `cpm_filter` is set, the correction is over the **surviving** genes only. | `target_mean` and `ref_mean` are always in natural (count) space regardless of `is_log1p` or `geometric_mean`. FDR is corrected within each group (across genes) for `ref` and `all` modes. For `on_target` mode, it is applied across all resulting p-values. +`epsilon` defaults to `1e-9`; pass `epsilon=0.0` to recover legacy `±inf` for one-sided zeros. +When `cpm_filter` is set, genes failing the floor are **dropped**, so the output has fewer than `n_targets × n_genes` rows (a fully-filtered comparison yields a height-0 frame with the full schema). ### Public API (`__all__`) diff --git a/pyproject.toml b/pyproject.toml index e83f869..9824f54 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "pdex" -version = "0.2.4" +version = "0.2.5" description = "Parallel differential expression for single-cell perturbation sequencing" readme = "README.md" authors = [{ name = "noam teyssier", email = "noam.teyssier@arcinstitute.org" }] diff --git a/src/pdex/__init__.py b/src/pdex/__init__.py index b01ba2a..f8fa9a7 100644 --- a/src/pdex/__init__.py +++ b/src/pdex/__init__.py @@ -12,7 +12,14 @@ from scipy.stats import false_discovery_control from tqdm import tqdm -from pdex._math import log2_fold_change, mwu, percent_change, pseudobulk +from pdex._math import ( + bulk_matrix_arithmetic, + cpm_bulk, + log2_fold_change, + mwu, + percent_change, + pseudobulk, +) from ._utils import _detect_is_log1p, set_numba_threadpool @@ -143,6 +150,42 @@ def _isolate_matrix( return np.asarray(result) +def _x_has_negative(x: np.ndarray | csr_matrix | None) -> bool: + """Best-effort check for negative values in the expression matrix. + + Negative values break the CPM computation used by ``cpm_filter`` (counts + cannot be negative). In-memory dense/sparse matrices are checked in full; + backed/other inputs fall back to a bounded sample of the leading rows. + """ + if x is None: + return False + if isinstance(x, csr_matrix): + return bool(x.data.size and (x.data < 0).any()) + if isinstance(x, np.ndarray): + return bool(x.size and (x < 0).any()) + arr = np.asarray(x[: min(1000, x.shape[0])]) + return bool(arr.size and (arr < 0).any()) + + +def _per_cell_library_sizes(adata: ad.AnnData, is_log1p: bool) -> np.ndarray: + """Per-cell total expression over all genes, in natural (count) space. + + Needed by ``on_target`` mode's CPM filter: the single-gene slices don't carry + the library size, so it is precomputed once over the full matrix (``expm1`` is + applied first when ``is_log1p``). Returns a 1-D float64 array of length n_obs. + """ + x = _isolate_matrix(adata, np.arange(adata.n_obs)) + if is_log1p: + if isinstance(x, csr_matrix): + x = x.copy() + np.expm1(x.data, out=x.data) + else: + x = np.expm1(np.asarray(x, dtype=np.float64)) + if isinstance(x, csr_matrix): + return np.asarray(x.sum(axis=1)).ravel().astype(np.float64) + return np.asarray(x, dtype=np.float64).sum(axis=1) + + def pdex( adata: ad.AnnData, groupby: str, @@ -151,7 +194,8 @@ def pdex( is_log1p: bool | None = None, geometric_mean: bool = True, as_pandas: bool = False, - epsilon: float = 0.0, + epsilon: float = 1e-9, + cpm_filter: float | None = None, **kwargs, ) -> pl.DataFrame | pd.DataFrame: """Run parallel differential expression analysis on single-cell data. @@ -203,23 +247,44 @@ def pdex( If ``True``, return a :class:`pandas.DataFrame` instead of a :class:`polars.DataFrame`. Requires ``pyarrow``. epsilon: - Pseudocount added to the denominator (and, for ``log2_fold_change``, the - numerator) before computing ``fold_change`` and ``percent_change``. When - ``epsilon > 0``, extreme values from near-zero reference means (scRNA-seq - sparsity artifact) are dampened toward zero. Has no effect on the - Mann-Whitney U p-value or FDR. Default ``0.0`` preserves existing behaviour; - regardless of ``epsilon``, features unexpressed in both groups report - ``0.0`` (no change) rather than ``NaN`` (see Returns). - - **Recommended usage:** For scRNA-seq CRISPRi/CRISPRa screens where many - genes are unexpressed in the reference group, start with ``epsilon=0.5``. - This provides modest dampening without substantially compressing fold changes - for well-expressed genes. For complete suppression of the sparsity artifact, - combine with a ``min_mean_expression`` pre-filter on the reference group — - ``epsilon`` alone cannot eliminate low p-values arising from per-cell - distributional shifts in near-zero genes. + 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 + 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 + the Mann-Whitney U p-value or FDR. Regardless of ``epsilon``, features + unexpressed in both groups report ``0.0`` (no change) rather than ``NaN`` + (see Returns). + + Default ``1e-9`` acts as a tiny finite-guard so the output contains no + ``±inf`` by default. Pass ``epsilon=0.0`` to recover the legacy behaviour + where one-sided zeros yield ``±inf``. For stronger dampening of the sparsity artifact in scRNA-seq + CRISPRi/CRISPRa screens, larger values (e.g. ``0.5``) trade fold-change + fidelity for floor suppression; prefer combining the tiny default with + ``cpm_filter`` to remove the noise floor outright. Must be non-negative. Raises :class:`ValueError` if negative. + cpm_filter: + Optional counts-per-million (CPM) floor filter. ``None`` (default) disables + it. When set to a threshold ``T``, a ``(target, feature)`` row is **dropped** + from the output when the gene's pooled (bulk) CPM is ``<= T`` in **both** the + target group and the reference (kept if ``target_cpm > T`` **or** + ``ref_cpm > T``). The pooled CPM is ``Σcounts_gene / Σcounts_all * 1e6`` per + group, computed on a separate internal CPM view (counts are recovered via + ``expm1`` when ``is_log1p``); the reported ``target_mean``/``ref_mean`` stay + in native count space — the output is never normalised. Because the CPM is a + ratio it is scale-invariant, so ``T`` means the same regardless of how the + input was normalised. The drop is independent of the Mann-Whitney U result, + and **FDR is corrected over the surviving genes only** (the filter changes + the multiple-testing universe). A negative ``T`` keeps everything. Emits a + :class:`UserWarning` if the data contains negative values (CPM assumes + non-negative expression). + + ``T = 5`` is a reasonable starting point, but the optimal threshold is + dataset-dependent (it tracks the separation between the noise floor and + genuinely expressed genes); inspect the per-gene CPM distribution of your + data and tune ``T`` empirically rather than relying on a fixed default. **kwargs: Mode-specific keyword arguments: @@ -297,6 +362,14 @@ def pdex( log.info("is_log1p=%s, geometric_mean=%s", is_log1p, geometric_mean) + if cpm_filter is not None and _x_has_negative(adata.X): # ty: ignore[invalid-argument-type] + msg = ( + "cpm_filter is set but adata.X contains negative values; CPM assumes " + "non-negative expression, so the filter results may be meaningless." + ) + log.warning(msg) + warnings.warn(msg, UserWarning, stacklevel=2) + if mode == "ref": reference = kwargs.pop("reference", DEFAULT_REFERENCE) if kwargs: @@ -313,6 +386,7 @@ def pdex( geometric_mean=geometric_mean, is_log1p=is_log1p, epsilon=epsilon, + cpm_filter=cpm_filter, ) elif mode == "all": if kwargs: @@ -327,6 +401,7 @@ def pdex( geometric_mean=geometric_mean, is_log1p=is_log1p, epsilon=epsilon, + cpm_filter=cpm_filter, ) elif mode == "on_target": gene_col = kwargs.pop("gene_col", None) @@ -348,6 +423,7 @@ def pdex( geometric_mean=geometric_mean, is_log1p=is_log1p, epsilon=epsilon, + cpm_filter=cpm_filter, ) else: raise ValueError(f"Invalid mode: {mode}") @@ -357,6 +433,78 @@ def pdex( return result +def _cpm_keep_mask( + target_matrix: np.ndarray | csr_matrix, + ref_cpm: np.ndarray, + is_log1p: bool, + cpm_filter: float, +) -> np.ndarray: + """Boolean keep mask: keep a gene iff target OR reference pooled CPM exceeds T. + + ``ref_cpm`` is precomputed by the caller (it is constant within a comparison). + """ + target_cpm = cpm_bulk(target_matrix, is_log1p) + return (target_cpm > cpm_filter) | (ref_cpm > cpm_filter) + + +def _assemble_group_frame( + *, + target: str, + feature_names: np.ndarray | pd.Index, + target_bulk: np.ndarray, + ref_bulk: np.ndarray, + target_membership: int, + ref_membership: int, + lfc: np.ndarray, + pc: np.ndarray, + pvalue: np.ndarray, + statistic: np.ndarray, + keep: np.ndarray | None, +) -> pl.DataFrame: + """Build one group's result frame, applying ``keep`` and correcting FDR over survivors. + + When ``keep`` is ``None`` no filtering is applied (legacy behaviour). The FDR is + always computed over the rows that remain, so it reflects the post-filter + multiple-testing universe. An all-``False`` mask yields a height-0 frame with the + full schema (safe for :func:`polars.concat`). + """ + feature = np.asarray(feature_names) + target_mean = np.asarray(target_bulk).ravel() + ref_mean = np.asarray(ref_bulk).ravel() + lfc = np.asarray(lfc).ravel() + pc = np.asarray(pc).ravel() + pvalue = np.asarray(pvalue).ravel() + statistic = np.asarray(statistic).ravel() + + if keep is not None: + feature = feature[keep] + target_mean = target_mean[keep] + ref_mean = ref_mean[keep] + lfc = lfc[keep] + pc = pc[keep] + pvalue = pvalue[keep] + statistic = statistic[keep] + + fdr = false_discovery_control(pvalue) if pvalue.size else np.empty(0, dtype=float) + + return pl.DataFrame( + { + "target": np.full(feature.shape[0], target), + "feature": feature, + "target_mean": target_mean, + "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, + "statistic": statistic, + "fdr": fdr, + } + ) + + def _pdex_ref( adata: ad.AnnData, groupby: str, @@ -364,6 +512,7 @@ def _pdex_ref( geometric_mean: bool = True, is_log1p: bool = False, epsilon: float = 0.0, + cpm_filter: float | None = None, ) -> pl.DataFrame: unique_groups, unique_group_indices = _unique_groups(adata.obs, groupby) log.info("Found %d groups (excluding reference)", len(unique_groups) - 1) @@ -376,6 +525,9 @@ def _pdex_ref( ref_bulk = pseudobulk(ref_matrix, geometric_mean=geometric_mean, is_log1p=is_log1p) ref_membership = ref_mask.size + # Reference pooled CPM is constant across target groups (CPM view, filter only) + ref_cpm = cpm_bulk(ref_matrix, is_log1p) if cpm_filter is not None else None + # Either sparse_column_index or ref_matrix ref_data = ( sparse_column_index(ref_matrix) @@ -405,24 +557,26 @@ def _pdex_ref( mwu_statistic = mwu_result.statistic mwu_pvalue = np.asarray(mwu_result.pvalue).clip(0, 1) - mwu_fdr = false_discovery_control(mwu_pvalue) + + if cpm_filter is None: + keep = None + else: + assert ref_cpm is not None # set whenever cpm_filter is not None + keep = _cpm_keep_mask(group_matrix, ref_cpm, is_log1p, cpm_filter) results.append( - pl.DataFrame( - { - "target": group_name, - "feature": feature_names, - "target_mean": np.asarray(group_bulk).ravel(), - "ref_mean": np.asarray(ref_bulk).ravel(), - "target_membership": group_mask.size, - "ref_membership": ref_membership, - "fold_change": lfc, - "log2_fold_change": lfc, - "percent_change": pc, - "p_value": mwu_pvalue, - "statistic": mwu_statistic, - "fdr": mwu_fdr, - } + _assemble_group_frame( + target=group_name, + feature_names=feature_names, + target_bulk=group_bulk, + ref_bulk=ref_bulk, + target_membership=group_mask.size, + ref_membership=ref_membership, + lfc=lfc, + pc=pc, + pvalue=mwu_pvalue, + statistic=mwu_statistic, + keep=keep, ) ) return pl.concat(results) @@ -434,6 +588,7 @@ def _pdex_all( geometric_mean: bool = True, is_log1p: bool = False, epsilon: float = 0.0, + cpm_filter: float | None = None, ) -> pl.DataFrame: unique_groups, unique_group_indices = _unique_groups(adata.obs, groupby) log.info("Found %d groups for 1-vs-rest comparison", len(unique_groups)) @@ -468,24 +623,26 @@ def _pdex_all( mwu_statistic = mwu_result.statistic mwu_pvalue = np.asarray(mwu_result.pvalue).clip(0, 1) - mwu_fdr = false_discovery_control(mwu_pvalue) + + if cpm_filter is None: + keep = None + else: + rest_cpm = cpm_bulk(rest_matrix, is_log1p) + keep = _cpm_keep_mask(group_matrix, rest_cpm, is_log1p, cpm_filter) results.append( - pl.DataFrame( - { - "target": group_name, - "feature": feature_names, - "target_mean": np.asarray(group_bulk).ravel(), - "ref_mean": np.asarray(rest_bulk).ravel(), - "target_membership": group_mask.size, - "ref_membership": rest_mask.size, - "fold_change": lfc, - "log2_fold_change": lfc, - "percent_change": pc, - "p_value": mwu_pvalue, - "statistic": mwu_statistic, - "fdr": mwu_fdr, - } + _assemble_group_frame( + target=group_name, + feature_names=feature_names, + target_bulk=group_bulk, + ref_bulk=rest_bulk, + target_membership=group_mask.size, + ref_membership=rest_mask.size, + lfc=lfc, + pc=pc, + pvalue=mwu_pvalue, + statistic=mwu_statistic, + keep=keep, ) ) @@ -500,6 +657,7 @@ def _pdex_on_target( geometric_mean: bool = True, is_log1p: bool = False, epsilon: float = 0.0, + cpm_filter: float | None = None, ) -> pl.DataFrame: unique_groups, unique_group_indices = _unique_groups(adata.obs, groupby) ref_index = _identify_reference_index(unique_groups, reference) @@ -520,6 +678,11 @@ def _pdex_on_target( len(group_gene_map), ) + # Per-cell library sizes for the CPM filter (single-gene slices lack them) + lib_cell = ( + _per_cell_library_sizes(adata, is_log1p) if cpm_filter is not None else None + ) + rows = [] for group_idx in tqdm( range(len(unique_groups)), @@ -544,6 +707,20 @@ def _pdex_on_target( group_col = np.asarray(group_col).reshape(-1, 1) ref_col = np.asarray(ref_col).reshape(-1, 1) + # CPM filter: drop the row iff the target gene is <= T in both sides. + # Uses arithmetic-mean pooled CPM (consistent with cpm_bulk), independent + # of the geometric/arithmetic choice for the reported means. + if cpm_filter is not None: + assert lib_cell is not None # set whenever cpm_filter is not None + target_arith = float(bulk_matrix_arithmetic(group_col, is_log1p)[0]) + ref_arith = float(bulk_matrix_arithmetic(ref_col, is_log1p)[0]) + t_lib = float(lib_cell[group_mask].mean()) if group_mask.size else 0.0 + r_lib = float(lib_cell[ref_mask].mean()) if ref_mask.size else 0.0 + target_cpm = target_arith / t_lib * 1e6 if t_lib > 0 else 0.0 + ref_cpm = ref_arith / r_lib * 1e6 if r_lib > 0 else 0.0 + if not (target_cpm > cpm_filter or ref_cpm > cpm_filter): + continue + target_mean = float( pseudobulk(group_col, geometric_mean=geometric_mean, is_log1p=is_log1p)[0] ) @@ -578,6 +755,26 @@ def _pdex_on_target( } ) + if not rows: + # No surviving rows (e.g. every target gene filtered out, or no group + # mapped to a gene): return a height-0 frame with the full schema. + return pl.DataFrame( + schema={ + "target": pl.Utf8, + "feature": pl.Utf8, + "target_mean": pl.Float64, + "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, + "statistic": pl.Float64, + "fdr": pl.Float64, + } + ) + df = pl.DataFrame(rows) fdr = false_discovery_control(df["p_value"].to_numpy()) return df.with_columns(pl.Series("fdr", fdr)) diff --git a/src/pdex/_math.py b/src/pdex/_math.py index 42fd6a3..c67e321 100644 --- a/src/pdex/_math.py +++ b/src/pdex/_math.py @@ -105,6 +105,30 @@ def bulk_matrix_geometric( return _expm1_vec(log_mean) +def cpm_bulk(matrix: np.ndarray | csr_matrix, is_log1p: bool) -> np.ndarray: + """Per-gene pooled (bulk) counts-per-million for one group's cell-by-gene matrix. + + Used only for the ``cpm_filter`` keep/drop decision — never for the reported + means. The arithmetic mean is taken in natural (count) space (back-transformed + via ``expm1`` when ``is_log1p``, reusing :func:`bulk_matrix_arithmetic`), then + + cpm[gene] = gene_mean[gene] / sum(gene_means) * 1e6 + + which equals the pooled ``Σcounts_gene / Σcounts_all * 1e6`` (the per-cell count + cancels). This is **scale-invariant**: a uniform rescaling of the input cancels + in the ratio, so a threshold in CPM means the same regardless of how the input + was normalised. + + A group whose total is zero (empty or all-zero slice) gets a denominator of 1.0 + so every gene's CPM is ``0.0`` (and is dropped by any positive threshold) rather + than ``inf``/``nan``. Returns a flat float64 array of length ``n_genes``. + """ + gene_means = bulk_matrix_arithmetic(matrix, is_log1p) + total = float(gene_means.sum()) + denom = total if total != 0.0 else 1.0 + return gene_means / denom * 1e6 + + @nb.njit(parallel=True) def log2_fold_change(x: np.ndarray, y: np.ndarray, epsilon: float = 0.0) -> np.ndarray: """Calculates the log2-fold change between two arrays. diff --git a/tests/conftest.py b/tests/conftest.py index 1472c30..cb225a3 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -64,6 +64,46 @@ def on_target_adata_sparse(on_target_adata): return adata +@pytest.fixture +def cpm_floor_adata(rng): + """AnnData purpose-built for the cpm_filter: 3 groups, 10 cells each, 5 genes. + + - gene_0..2: well expressed in every group (high CPM, always kept). + - gene_3: one-sided — zero in the non-targeting reference, expressed in A and B + (kept by the OR rule; produces a one-sided-zero LFC). + - gene_4: pure floor — zero in every group (dropped whenever both sides <= T). + """ + n_per = 10 + groups = ["non-targeting", "A", "B"] + obs_groups = np.repeat(groups, n_per) + n_cells = len(obs_groups) + n_genes = 5 + + X = rng.poisson(lam=5, size=(n_cells, n_genes)).astype(np.float64) + X[n_per : 2 * n_per, :3] += 3 # boost genes 0..2 in A + X[2 * n_per :, :3] += 6 # boost genes 0..2 in B + + # gene_3: zero in reference, expressed in A and B (one-sided) + X[:n_per, 3] = 0.0 + # gene_4: zero everywhere (pure floor) + X[:, 4] = 0.0 + + obs = pd.DataFrame( + {"guide": obs_groups}, + index=np.array([f"cell_{i}" for i in range(n_cells)]), + ) + var = pd.DataFrame(index=np.array([f"gene_{i}" for i in range(n_genes)])) + return ad.AnnData(X=X, obs=obs, var=var) + + +@pytest.fixture +def cpm_floor_adata_sparse(cpm_floor_adata): + """cpm_floor_adata with sparse CSR X matrix.""" + adata = cpm_floor_adata.copy() + adata.X = csr_matrix(adata.X) + return adata + + @pytest.fixture def small_adata_log1p(small_adata): """small_adata with X replaced by log1p-transformed values.""" diff --git a/tests/test_math.py b/tests/test_math.py index 79344ea..171ecf6 100644 --- a/tests/test_math.py +++ b/tests/test_math.py @@ -1,8 +1,14 @@ """Tests for pdex._math (log2_fold_change, percent_change, bulk_matrix_geometric).""" import numpy as np +from scipy.sparse import csr_matrix -from pdex._math import bulk_matrix_geometric, log2_fold_change, percent_change +from pdex._math import ( + bulk_matrix_geometric, + cpm_bulk, + log2_fold_change, + percent_change, +) class TestFoldChange: @@ -186,3 +192,44 @@ def test_zeros_handled(self): result = bulk_matrix_geometric(x, is_log1p=False) expected = np.expm1(np.log1p(x).mean(axis=0)) np.testing.assert_allclose(result, expected) + + +class TestCpmBulk: + """Tests for cpm_bulk (pooled CPM view used by the cpm_filter).""" + + def test_known_values_dense(self): + """cpm[g] = Σcounts_g / Σcounts_all * 1e6.""" + x = np.array([[1.0, 3.0], [3.0, 5.0], [5.0, 1.0]]) + # gene sums: 9, 9; total 18 -> 0.5 each + result = cpm_bulk(x, is_log1p=False) + np.testing.assert_allclose(result, [500_000.0, 500_000.0]) + + def test_sums_to_one_million(self): + x = np.array([[2.0, 4.0, 0.0], [6.0, 1.0, 3.0]]) + np.testing.assert_allclose(cpm_bulk(x, is_log1p=False).sum(), 1e6) + + def test_sparse_matches_dense(self): + x = np.array([[1.0, 3.0, 0.0], [3.0, 5.0, 2.0], [5.0, 1.0, 0.0]]) + dense = cpm_bulk(x, is_log1p=False) + sparse = cpm_bulk(csr_matrix(x), is_log1p=False) + np.testing.assert_allclose(dense, sparse) + + def test_log1p_agrees_with_counts(self): + """is_log1p=True on log1p(counts) matches is_log1p=False on counts.""" + counts = np.array([[2.0, 4.0], [6.0, 8.0], [10.0, 2.0]]) + from_counts = cpm_bulk(counts, is_log1p=False) + from_log = cpm_bulk(np.log1p(counts), is_log1p=True) + np.testing.assert_allclose(from_counts, from_log, rtol=1e-10) + + def test_all_zero_is_zero_not_nan(self): + """An all-zero group yields all-zero CPM (denominator guard), never nan/inf.""" + x = np.zeros((3, 4)) + result = cpm_bulk(x, is_log1p=False) + assert not np.isnan(result).any() + assert not np.isinf(result).any() + np.testing.assert_array_equal(result, np.zeros(4)) + + def test_scale_invariant(self): + """Uniformly rescaling counts does not change the CPM (ratio cancels).""" + x = np.array([[1.0, 3.0, 7.0], [3.0, 5.0, 2.0], [5.0, 1.0, 4.0]]) + np.testing.assert_allclose(cpm_bulk(x, False), cpm_bulk(100.0 * x, False)) diff --git a/tests/test_pdex.py b/tests/test_pdex.py index 08bc14d..2a74e54 100644 --- a/tests/test_pdex.py +++ b/tests/test_pdex.py @@ -143,11 +143,11 @@ def test_epsilon_accepted(self, small_adata): result = pdex(small_adata, groupby="guide", is_log1p=False, epsilon=0.5) assert isinstance(result, pl.DataFrame) - def test_epsilon_zero_matches_default(self, small_adata): - """epsilon=0.0 produces identical results to omitting the parameter.""" + def test_default_epsilon_is_tiny_finite_guard(self, small_adata): + """Omitting epsilon uses the 1e-9 default (not 0.0).""" default_result = pdex(small_adata, groupby="guide", is_log1p=False) explicit_result = pdex( - small_adata, groupby="guide", is_log1p=False, epsilon=0.0 + small_adata, groupby="guide", is_log1p=False, epsilon=1e-9 ) assert isinstance(default_result, pl.DataFrame) assert isinstance(explicit_result, pl.DataFrame) @@ -736,3 +736,222 @@ def test_one_sided_zero_still_infinite(self, small_adata): # log2(0 / ref) -> -inf; percent_change (0 - ref) / ref -> -1.0 assert np.isneginf(gene0["log2_fold_change"].to_numpy()).all() np.testing.assert_allclose(gene0["percent_change"].to_numpy(), -1.0) + + +def _pairs(df) -> set: + """Set of (target, feature) tuples in a result frame.""" + return set(zip(df["target"].to_list(), df["feature"].to_list())) + + +class TestCpmFilter: + """Tests for the cpm_filter (bulk-CPM floor) parameter.""" + + def test_both_sides_below_threshold_dropped(self, cpm_floor_adata): + """gene_4 (zero in every group) is dropped from the output.""" + result = pdex( + cpm_floor_adata, groupby="guide", mode="ref", is_log1p=False, cpm_filter=5 + ) + assert "gene_4" not in result["feature"].to_list() + + def test_one_side_above_threshold_kept(self, cpm_floor_adata): + """gene_3 (zero in ref, expressed in target) survives via the OR rule.""" + result = pdex( + cpm_floor_adata, groupby="guide", mode="ref", is_log1p=False, cpm_filter=5 + ) + assert "gene_3" in result["feature"].to_list() + + def test_negative_threshold_keeps_everything(self, cpm_floor_adata): + """A negative threshold keeps all genes (CPM >= 0 > T).""" + unfiltered = pdex(cpm_floor_adata, groupby="guide", mode="ref", is_log1p=False) + result = pdex( + cpm_floor_adata, groupby="guide", mode="ref", is_log1p=False, cpm_filter=-1 + ) + assert _pairs(result) == _pairs(unfiltered) + + def test_zero_threshold_strict_drops_only_exact_zero(self, cpm_floor_adata): + """T=0 drops genes whose pooled CPM is exactly 0 (strict >), keeps the rest.""" + result = pdex( + cpm_floor_adata, groupby="guide", mode="ref", is_log1p=False, cpm_filter=0.0 + ) + features = set(result["feature"].to_list()) + assert "gene_4" not in features # cpm 0, not > 0 -> dropped + assert {"gene_0", "gene_1", "gene_2", "gene_3"} <= features + + def test_none_matches_unfiltered(self, cpm_floor_adata): + """cpm_filter=None is identical to omitting it.""" + omitted = pdex(cpm_floor_adata, groupby="guide", mode="ref", is_log1p=False) + explicit = pdex( + cpm_floor_adata, + groupby="guide", + mode="ref", + is_log1p=False, + cpm_filter=None, + ) + assert isinstance(omitted, pl.DataFrame) + assert isinstance(explicit, pl.DataFrame) + assert omitted.equals(explicit) + + def test_filtered_is_subset_with_unchanged_values(self, cpm_floor_adata): + """Surviving rows keep their exact means/p-values; only the row set shrinks.""" + unfiltered = pdex(cpm_floor_adata, groupby="guide", mode="ref", is_log1p=False) + filtered = pdex( + cpm_floor_adata, groupby="guide", mode="ref", is_log1p=False, cpm_filter=5 + ) + assert _pairs(filtered) < _pairs(unfiltered) # strict subset + # Surviving rows keep their exact target_mean / p_value (look up by key) + full = { + (t, f): (tm, pv) + for t, f, tm, pv in zip( + unfiltered["target"].to_list(), + unfiltered["feature"].to_list(), + unfiltered["target_mean"].to_list(), + unfiltered["p_value"].to_list(), + ) + } + for t, f, tm, pv in zip( + filtered["target"].to_list(), + filtered["feature"].to_list(), + filtered["target_mean"].to_list(), + filtered["p_value"].to_list(), + ): + np.testing.assert_allclose(tm, full[(t, f)][0]) + np.testing.assert_allclose(pv, full[(t, f)][1]) + + def test_scale_invariant_kept_set(self, cpm_floor_adata): + """Uniformly rescaling counts does not change which rows survive.""" + scaled = cpm_floor_adata.copy() + scaled.X = scaled.X * 100.0 + base = pdex( + cpm_floor_adata, groupby="guide", mode="ref", is_log1p=False, cpm_filter=5 + ) + rescaled = pdex( + scaled, groupby="guide", mode="ref", is_log1p=False, cpm_filter=5 + ) + assert _pairs(base) == _pairs(rescaled) + + def test_no_inf_with_default_epsilon(self, cpm_floor_adata): + """Filter + default epsilon (1e-9) leaves no inf/nan in the LFC columns.""" + result = pdex( + cpm_floor_adata, groupby="guide", mode="ref", is_log1p=False, cpm_filter=5 + ) + for col in ("log2_fold_change", "percent_change"): + vals = result[col].to_numpy() + assert not np.isinf(vals).any(), col + assert not np.isnan(vals).any(), col + + def test_epsilon_zero_still_allows_one_sided_inf(self, cpm_floor_adata): + """With epsilon=0, a surviving one-sided zero (gene_3) keeps its +inf LFC.""" + result = pdex( + cpm_floor_adata, + groupby="guide", + mode="ref", + is_log1p=False, + cpm_filter=5, + epsilon=0.0, + ) + gene3 = result.filter(pl.col("feature") == "gene_3") + # gene_3: ref_mean 0, target_mean > 0 -> log2(t/0) = +inf + assert np.isposinf(gene3["log2_fold_change"].to_numpy()).all() + + def test_fdr_over_survivors(self, cpm_floor_adata): + """FDR is recomputed over surviving genes, not the full gene set.""" + filtered = pdex( + cpm_floor_adata, groupby="guide", mode="ref", is_log1p=False, cpm_filter=5 + ) + a = filtered.filter(pl.col("target") == "A") + # FDR matches BH over exactly the surviving p-values + recomputed = stats.false_discovery_control(a["p_value"].to_numpy()) + np.testing.assert_allclose(a["fdr"].to_numpy(), recomputed) + + # ... and differs from BH over the full (unfiltered) p-value set + unfiltered = pdex(cpm_floor_adata, groupby="guide", mode="ref", is_log1p=False) + a_full = unfiltered.filter(pl.col("target") == "A") + full_fdr = dict( + zip( + a_full["feature"].to_list(), + stats.false_discovery_control(a_full["p_value"].to_numpy()), + ) + ) + surviving_full = np.array([full_fdr[f] for f in a["feature"].to_list()]) + assert not np.allclose(a["fdr"].to_numpy(), surviving_full) + + def test_sparse_matches_dense(self, cpm_floor_adata, cpm_floor_adata_sparse): + """Sparse input yields the same kept set and values as dense.""" + dense = pdex( + cpm_floor_adata, groupby="guide", mode="ref", is_log1p=False, cpm_filter=5 + ) + sparse = pdex( + cpm_floor_adata_sparse, + groupby="guide", + mode="ref", + is_log1p=False, + cpm_filter=5, + ) + assert _pairs(dense) == _pairs(sparse) + d = dense.sort(["target", "feature"]) + s = sparse.sort(["target", "feature"]) + np.testing.assert_allclose( + d["target_mean"].to_numpy(), s["target_mean"].to_numpy() + ) + + def test_log1p_kept_set_matches_raw(self, cpm_floor_adata): + """CPM is computed on counts, so log1p input gives the same kept set.""" + log_adata = cpm_floor_adata.copy() + log_adata.X = np.log1p(log_adata.X) + raw = pdex( + cpm_floor_adata, groupby="guide", mode="ref", is_log1p=False, cpm_filter=5 + ) + logged = pdex( + log_adata, groupby="guide", mode="ref", is_log1p=True, cpm_filter=5 + ) + assert _pairs(raw) == _pairs(logged) + + def test_all_mode_drops_floor(self, cpm_floor_adata): + """In 1-vs-rest mode, the floor gene is dropped for every group.""" + result = pdex( + cpm_floor_adata, groupby="guide", mode="all", is_log1p=False, cpm_filter=5 + ) + assert "gene_4" not in result["feature"].to_list() + # every group still present (they retain expressed genes) + assert set(result["target"].to_list()) == {"non-targeting", "A", "B"} + + def test_on_target_drops_floor_group(self, cpm_floor_adata): + """on_target: a group whose target gene is a floor gene is dropped.""" + adata = cpm_floor_adata.copy() + adata.obs["target_gene"] = ( + adata.obs["guide"].map( + {"non-targeting": "gene_0", "A": "gene_3", "B": "gene_4"} + ) + ).astype(object) + result = pdex( + adata, + groupby="guide", + mode="on_target", + gene_col="target_gene", + is_log1p=False, + cpm_filter=5, + ) + # A targets gene_3 (expressed) -> kept; B targets gene_4 (floor) -> dropped + assert result["target"].to_list() == ["A"] + assert result["feature"].to_list() == ["gene_3"] + # FDR over the single surviving row + assert (result["fdr"] >= 0).all() and (result["fdr"] <= 1).all() + + def test_all_dropped_returns_empty_with_schema(self, cpm_floor_adata): + """A threshold above every gene's CPM yields a height-0 frame, full schema.""" + result = pdex( + cpm_floor_adata, + groupby="guide", + mode="ref", + is_log1p=False, + cpm_filter=1e12, + ) + assert result.height == 0 + assert set(result.columns) == EXPECTED_COLUMNS + + def test_negative_values_warn(self, cpm_floor_adata): + """cpm_filter on data with negative values emits a UserWarning.""" + adata = cpm_floor_adata.copy() + adata.X[0, 0] = -1.0 + with pytest.warns(UserWarning, match="negative values"): + pdex(adata, groupby="guide", mode="ref", is_log1p=False, cpm_filter=5) From cfe150aada8e6ec228cbd47df829ac1505977518 Mon Sep 17 00:00:00 2001 From: Abhinav Adduri Date: Fri, 26 Jun 2026 18:11:56 +0000 Subject: [PATCH 2/3] fix: handle backed sparse cpm filter preflight --- src/pdex/__init__.py | 15 +++++++----- tests/test_pdex.py | 57 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 66 insertions(+), 6 deletions(-) diff --git a/src/pdex/__init__.py b/src/pdex/__init__.py index f8fa9a7..e9a7e56 100644 --- a/src/pdex/__init__.py +++ b/src/pdex/__init__.py @@ -1,6 +1,6 @@ import logging import warnings -from typing import Literal +from typing import Any, Literal, cast import anndata as ad import numpy as np @@ -159,11 +159,14 @@ def _x_has_negative(x: np.ndarray | csr_matrix | None) -> bool: """ if x is None: return False - if isinstance(x, csr_matrix): - return bool(x.data.size and (x.data < 0).any()) - if isinstance(x, np.ndarray): - return bool(x.size and (x < 0).any()) - arr = np.asarray(x[: min(1000, x.shape[0])]) + sample = x + if not isinstance(sample, (np.ndarray, csr_matrix)): + sample_obj = cast(Any, sample) + sample = sample_obj[: min(1000, sample_obj.shape[0])] + if issparse(sample): + sample = csr_matrix(sample) + return bool(sample.data.size and (sample.data < 0).any()) + arr = np.asarray(sample) return bool(arr.size and (arr < 0).any()) diff --git a/tests/test_pdex.py b/tests/test_pdex.py index 2a74e54..c0eb5f7 100644 --- a/tests/test_pdex.py +++ b/tests/test_pdex.py @@ -1,9 +1,11 @@ """Integration tests for pdex() and _pdex_ref().""" +import anndata as ad import numpy as np import polars as pl import pytest from scipy import stats +from scipy.sparse import csr_matrix from pdex import DEFAULT_REFERENCE, pdex @@ -894,6 +896,48 @@ def test_sparse_matches_dense(self, cpm_floor_adata, cpm_floor_adata_sparse): d["target_mean"].to_numpy(), s["target_mean"].to_numpy() ) + def test_backed_sparse_matches_inmemory_sparse( + self, cpm_floor_adata_sparse, tmp_path + ): + """Backed sparse input supports cpm_filter and matches in-memory sparse.""" + path = tmp_path / "cpm_floor_sparse.h5ad" + cpm_floor_adata_sparse.write_h5ad(path) + backed = ad.read_h5ad(path, backed="r") + + inmem = pdex( + cpm_floor_adata_sparse, + groupby="guide", + mode="ref", + is_log1p=False, + cpm_filter=5, + ) + backed_result = pdex( + backed, + groupby="guide", + mode="ref", + is_log1p=False, + cpm_filter=5, + ) + + assert _pairs(inmem) == _pairs(backed_result) + i = inmem.sort(["target", "feature"]) + b = backed_result.sort(["target", "feature"]) + for col in [ + "target_mean", + "ref_mean", + "fold_change", + "percent_change", + "p_value", + "statistic", + "fdr", + ]: + np.testing.assert_allclose( + i[col].to_numpy(), + b[col].to_numpy(), + rtol=1e-6, + err_msg=f"Mismatch in column {col}", + ) + def test_log1p_kept_set_matches_raw(self, cpm_floor_adata): """CPM is computed on counts, so log1p input gives the same kept set.""" log_adata = cpm_floor_adata.copy() @@ -955,3 +999,16 @@ def test_negative_values_warn(self, cpm_floor_adata): adata.X[0, 0] = -1.0 with pytest.warns(UserWarning, match="negative values"): pdex(adata, groupby="guide", mode="ref", is_log1p=False, cpm_filter=5) + + def test_backed_sparse_negative_values_warn(self, cpm_floor_adata, tmp_path): + """The negative-value preflight handles backed sparse arrays.""" + adata = cpm_floor_adata.copy() + adata.X[0, 0] = -0.5 + adata.X = csr_matrix(adata.X) + path = tmp_path / "negative_sparse.h5ad" + adata.write_h5ad(path) + backed = ad.read_h5ad(path, backed="r") + + with pytest.warns(UserWarning, match="negative values"): + with pytest.raises(ValueError, match="Sparse MWU requires non-negative"): + pdex(backed, groupby="guide", mode="ref", is_log1p=False, cpm_filter=5) From 7986044ce25d1a894bd02d5d64ca5677943379ba Mon Sep 17 00:00:00 2001 From: Abhinav Adduri Date: Fri, 26 Jun 2026 19:22:20 +0000 Subject: [PATCH 3/3] docs: document cpm filter and epsilon semantics --- README.md | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index 5d891a6..3f54929 100644 --- a/README.md +++ b/README.md @@ -80,9 +80,17 @@ results = pdex( | `is_log1p` | `bool \| None` | `None` | Whether data is log1p-transformed. Auto-detected if `None` | | `geometric_mean` | `bool` | `True` | Use geometric mean for pseudobulk (vs arithmetic) | | `as_pandas` | `bool` | `False` | Return a pandas DataFrame instead of Polars | +| `epsilon` | `float` | `1e-9` | Pseudocount used for `log2_fold_change` and `percent_change`; pass `0.0` for legacy one-sided `±inf` values | +| `cpm_filter` | `float \| None` | `None` | Optional pooled-CPM floor filter; drops rows where both target and reference CPM are at or below the threshold | | `reference` | `str` | `"non-targeting"` | Reference group name (modes: `ref`, `on_target`) | | `gene_col` | `str` | — | Column mapping groups to target genes (mode: `on_target`) | +### CPM filter + +`cpm_filter` is an opt-in floor filter. When set to a threshold `T`, a `(target, feature)` row is dropped only when the gene's pooled CPM is `<= T` in both the target group and the reference group. Rows are kept when either side has CPM `> T`. + +The CPM view is used only for filtering: reported means, fold changes, MWU statistics, and p-values are still computed from the original expression scale. When `is_log1p=True`, counts are recovered with `expm1` before CPM is computed. FDR is corrected over the surviving genes only. + ## Output Returns a Polars DataFrame (or pandas if `as_pandas=True`) with one row per (group, gene) pair: @@ -96,8 +104,8 @@ Returns a Polars DataFrame (or pandas if `as_pandas=True`) with one row per (gro | `target_membership` | Number of cells in the target group | | `ref_membership` | Number of cells in the reference | | `fold_change` | **Deprecated alias** for `log2_fold_change` (identical values). Will be removed in pdex 0.3.0. | -| `log2_fold_change` | log2(target_mean / ref_mean). Genes unexpressed in both groups (0/0) report `0.0`, not `NaN`. | -| `percent_change` | (target_mean - ref_mean) / ref_mean. Genes unexpressed in both groups (0/0) report `0.0`, not `NaN`. | +| `log2_fold_change` | log2((target_mean + epsilon) / (ref_mean + epsilon)). With default `epsilon=1e-9`, one-sided zeros are large finite values; with `epsilon=0.0`, one-sided zeros yield `±inf`. Genes unexpressed in both groups (0/0) report `0.0`, not `NaN`. | +| `percent_change` | (target_mean - ref_mean) / (ref_mean + epsilon). With default `epsilon=1e-9`, zero-reference cases are finite; with `epsilon=0.0`, a zero reference with nonzero target yields `+inf`. Genes unexpressed in both groups (0/0) report `0.0`, not `NaN`. | | `p_value` | Mann-Whitney U p-value | | `statistic` | Mann-Whitney U statistic | -| `fdr` | FDR-corrected p-value (per-group, across genes). For `on_target` mode, this is applied across all groups. | +| `fdr` | FDR-corrected p-value (per-group, across genes). For `on_target` mode, this is applied across all groups. When `cpm_filter` is set, FDR is corrected over surviving genes only. |