Add feature-complete Python port of SelectSim (v0.3.0)#1
Conversation
- Add utils.py with helper functions (add, pairwise_indices, matrix_to_pairwise_vector, create_pair_template) - Add template.py with generate_s and template_obj_gen functions - Update pyproject.toml with new dependencies (joblib, tqdm, pytest, pytest-cov, pyreadr)
- Add weights.py with generate_w_mean_tmb and generate_w_block functions - Implements TMB-based sample weighting using fold-change penalty
- Add null_model.py with null_model_parallel and retrieve_outliers functions - Implements rejection sampling algorithm preserving gene frequencies - Supports parallel execution via joblib - Includes outlier detection and removal (top 10%)
- Add stats.py with all statistical computation functions - Includes am_stats, al_stats, pairwise overlap functions - Implements effect_size, estimate_fdr2, binary_yule - Adds interaction_table for generating final results with 23 columns - Supports CO/ME classification and FDR computation
- Implement full 13-step analysis pipeline in selectsim.py - Fix parameter order bug in AlterationLandscape constructor call - Add module-level docstrings to alteration_landscape.py - Update __init__.py with all public exports - Set package version to 0.1.0
- Add conftest.py with test fixtures (simple_test_data, luad_run_data) - Add test_alteration_landscape.py for AlterationLandscape class - Add test_template.py for template generation functions - Add test_weights.py for weight matrix generation - Add test_null_model.py for null model simulation - Add test_stats.py for statistics functions - Add test_selectsim.py for integration tests - Tests use pytest framework with synthetic data
- Add comprehensive documentation in README.md - Include installation instructions - Add quick start example - Document key parameters and output format - Add API reference section - Include development instructions
- Fix exp_r_es_norm indexing from [select, :] to [:, select] to correctly select pairs (columns) not permutations (rows) - Widen test tolerance for outlier fraction (10% ± 10%)
… for null model Ports the remaining pieces of the R SelectSim package that had no Python equivalent: gam_utils.r (MAF filtering + GAM construction) as selectsim.gam, and selectX_plot.R (ridge plots, obs/exp scatter) as selectsim.plotting (matplotlib/seaborn). Adds selectsim.io for Parquet/Zarr helpers, and an optional backend="jax" / store="zarr" path on null_model_parallel (default numpy/in-memory behavior unchanged) to vectorize the per-gene selection loop that was the main hot path. Also exports the previously-unexported estimate_p_val/estimate_pairwise_p and fixes a malformed docstring in alteration_landscape.py that broke Sphinx builds.
Exports the R package's example datasets (luad_run_data, luad_result, luad_maf, oncokb gene lists, variant_catalogue) from their .rda form to Parquet (scripts/export_r_data.R + csv_to_parquet.py), bundled under tests/data/ so tests no longer need a sibling R checkout or pyreadr. Rewrites the LUAD fixtures in conftest.py to load the bundled Parquet files, and replaces the two previously-skipped, empty-bodied parity tests in test_selectsim.py with real assertions: deterministic columns (overlap, freq, support) match R to ~1e-12, and permutation-dependent columns (nES) correlate at 0.9998 with R's reference output.
…scaffolding Full API reference via autosummary (one short page per function instead of per-module automodule dumps), two executable Jupyter tutorials mirroring the R vignettes (introduction, data_processing), and a switch from sphinx-rtd-theme to sphinx-book-theme (matching the sibling CellCharter project's docs), with sphinx-autodoc-typehints and sphinx-copybutton. Adds standard OSS scaffolding that was entirely missing: LICENSE (MIT), CITATION.cff, NEWS.md, logo/favicon assets, and .github/ (CI workflow for tests, a docs-build-and-deploy workflow that's ready but inert until the repo is public, and issue templates).
Rewrites pyproject.toml from Poetry's legacy [tool.poetry.*] tables to standard PEP 621 ([project], [project.optional-dependencies]) plus PEP 735 [dependency-groups] for dev-only tooling, switching the build backend from poetry-core to hatchling. Fixes a latent bug found in the process: the old docs extra referenced "Sphinx" (wrong case) instead of the dependency table key "sphinx", which made `poetry check` fail silently. poetry.lock is replaced by uv.lock. `uv sync` now creates an isolated project-local .venv/ instead of Poetry writing this package's dependencies directly into whatever environment happened to be active. Verified 126/126 tests pass on Python 3.10, 3.11, and 3.12.
Both vertical segments (red = actual overlap, blue = mean background) were scaled to the KDE density *at that exact x-position*, matching R's own `ymin + density * scale * iscale`, except R evaluates that expression at the ridge's peak-density row and reuses it for both segments regardless of where x falls, while the Python port re-evaluated density at each marker's own x via np.interp. Since the actual overlap is typically in the tail of the null distribution for a significant pair, its local density is near zero, making the red line invisible for exactly the pairs it matters most for. Both ridge_plot_ed and ridge_plot_ed_compare now draw both segments to the ridge's peak height, matching R's actual behavior.
filter_maf_column(inclusive=False, fixed=False) with multiple values negated each value's containment mask independently and unioned the results, so a row only got dropped if it matched every excluded value simultaneously rather than any one of them (e.g. a row containing "A" but not "B" survived a values=["A","B"] exclusion, when it should have been dropped for containing "A"). Found by code review, live-reproduced, fixed by building the "contains any excluded value" mask first and negating once. Also adds a dedupe=True parameter (propagated through the filter_maf_* wrapper functions via **kwargs) so callers chaining several filter calls can opt out of the now-redundant per-step drop_duplicates() scan after the first one -- pure boolean-mask filtering can't introduce a duplicate that wasn't already in the input, so dedup only ever needs to happen once. Co-Authored-By: Claude <noreply@anthropic.com>
store="zarr" fully materialized every permutation in memory before ever writing to disk, contradicting its own docstring and defeating the reason to use it -- now streams in bounded batches (default 200 permutations), confirmed ~5x lower peak memory on a 2000-permutation/100-gene/300-sample run, with bit-identical output to store="memory" for the same seed (verified for both the numpy and jax backends, including that chunking the jax backend preserves the exact same PRNGKey stream as a single-shot call). The jax backend's jit-compiled kernel was rebuilt as a fresh closure on every null_model_parallel(backend="jax") call, so JAX's compilation cache could never hit across calls with identical shapes -- moved to a module-level cache built once, with templates/gene_freq passed as explicit (non-closed-over) arguments so the cached function works across datasets. n_cores now warns instead of silently no-op'ing when backend="jax" (jax parallelizes via vmap, not multiprocessing, so the argument has no effect). retrieve_outliers: n_sim was accepted but never used (len(null) always used instead); now actually sub-samples. Its 90th-percentile outlier threshold was also off-by-one vs R's 1-indexed dev2[round(0.90*length)] formula, making Python's cutoff one rank stricter than R's; aligned to match. ridge_plot_ed crashed with "max() iterable argument is empty" whenever every displayed pair's null-model background was empty (e.g. zero retained permutations); fixed with the same observed-value/literal fallback ridge_plot_ed_compare already used. Also factored the peak-height marker geometry shared by both ridge functions into one helper (_ridge_row_geometry) so a future fix to that logic only needs to be made once. Co-Authored-By: Claude <noreply@anthropic.com>
…NEWS Bumps version to 0.2.0 (pyproject.toml, selectsim.__version__, CITATION.cff, docs release) -- this round changed observable behavior (store="zarr" now genuinely streams, n_cores warns with backend="jax", retrieve_outliers' outlier flagging shifts slightly), which under 0.x semver belongs in a minor bump rather than a patch release. Trims README.md down to match the R package's README structure and length (title, method figure, citation, install, quick start, doc links, citation, license, contact) -- it had grown a full API reference and parameter/output sections that now duplicate the Sphinx docs, which didn't exist when the original README was written. Adds a NEWS.md note that this round of work (the bug fixes in the previous commit, plus the earlier module ports, documentation, and uv migration) was developed with Claude Code assistance, directed and verified by Arvind Iyer throughout -- kept in NEWS.md rather than README to keep the README itself minimal, matching the R package's style. Co-Authored-By: Claude <noreply@anthropic.com>
New selectsim.plotting.oncoprint_pair(gene1, gene2, obj, simulation_index=0): a two-gene "mini oncoprint" with a per-sample TMB bar, stacking the real (observed) co-mutation pattern directly above the pattern from one retained null-model permutation, illustrating the impact of TMB on co-mutation. Requested and iteratively refined against a reference figure from the accompanying manuscript; verified the annotated counts match the real package pipeline exactly (am_weight_pairwise_alteration_overlap) before that annotation was later simplified away per feedback. Design, settled after several rounds of feedback: - Samples grouped by mutation status (co-mutated / primary-gene-only / neither / secondary-gene-only), each group keeping its original column order -- no TMB or other sort is applied anywhere (verified in tests by spying on numpy.argsort/numpy.sort). - TMB row: a single fill_between step-polygon rather than one bar() patch per sample, which anti-aliases into visible seams at ~500 samples. - Gene rows: an imshow raster with a transparent "unmutated" color and a thin black sample-track baseline drawn underneath, showing through wherever a gene isn't altered (standard oncoprint convention) -- mutated cells are opaque and cover the line. - Highlight: translucent fill on the TMB row, clean (unfilled) outline on the gene rows so the green/white track colors don't get tinted. - Deliberately minimal: no annotation text/arrows (that's left to the user, e.g. via Illustrator, as in the original figure) -- just the highlight and a plain co-mutation count. Co-Authored-By: Claude <noreply@anthropic.com>
New selectsim.plotting.oncoprint(genes, obj, simulation_index=0, title=None):
generalizes oncoprint_pair from two genes to an arbitrary list, in the same
visual style -- a per-sample TMB bar above one binary track per gene, with
the "Observed" panel and one retained null-model permutation side by side.
- Genes ordered top-to-bottom by decreasing observed mutation frequency
(shared between panels, since the null model preserves each gene's
frequency exactly); samples ordered by a "memo sort" generalizing
oncoprint_pair's 2-gene grouping to N genes (mutation pattern across
genes, most-mutated-gene first, no TMB or other sort -- verified in
tests by spying on numpy.argsort/numpy.sort).
- One matplotlib axes per gene (not a single combined heatmap image), so
gene rows get real vertical spacing via hspace, matching oncoprint_pair.
- Gene names and the "TMB" axis label are only shown once, on the
Observed (left) panel -- both are shared/identical across panels, and
repeating them on the adjacent Simulation panel just overlapped it.
- Panel title shows which specific permutation is displayed (e.g.
"Simulation (permutation #42)"), in both oncoprint and oncoprint_pair,
since a single instance out of many retained permutations shouldn't be
captioned as if it were "the" null result.
- title defaults to f"Oncoprint (n={n_samples})"; pass something dataset-
specific (e.g. "LUAD oncoprint (n=502)") since the function has no way
to know the cohort's name on its own.
- "TMB" axis label no longer inherits theme_publication's bold axis-label
weight, in both oncoprint and oncoprint_pair.
Co-Authored-By: Claude <noreply@anthropic.com>
…ax backend The numpy null-model backend (_simulation_fixed_ones) looped over genes in Python to pick each row's top-k residual columns; replaced with the same vectorized double-argsort rank trick the jax backend already used, eliminating the per-gene Python loop entirely. Verified bit-identical to the previous implementation (full R-parity suite still passes). store="zarr" runs defaulted the on-disk array's chunk_permut to 1 while writing in batches of up to 200 permutations, forcing each batch write to be split into up to 200 separate compress/write calls. Chunking is now aligned to the write-batch size; benchmarked at ~2.5x faster zarr writes in isolation. Benchmarked the jax backend against the (now vectorized) numpy backend on real LUAD-sized data, CPU-only (no GPU available): jax was consistently 5-10x *slower*, dominated by per-call JIT dispatch overhead that this workload's array sizes never amortize. Removed backend="jax" from null_model_parallel, the selectsim[fast] extra, and jax mentions in README/docs, rather than keep a backend that's strictly worse here. If GPU acceleration is needed later, benchmark against the vectorized numpy backend on real GPU hardware before reintroducing it. Also fixed estimate_p_val/estimate_pairwise_p (used when estimate_pairwise=True): gene_names.index(...) did an O(n_genes) list scan per gene per pair; estimate_p_val now takes a prebuilt gene_index: Dict[str, int] built once by its caller instead. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018g7ADBXHN2yWDkf2NPyNm3
Upper-triangular gene x gene heatmap of interaction type (CO/ME/not significant) and strength (|nES| as color intensity) across all tested pairs at once -- the same visual idiom as maftools::somaticInteractions / cBioPortal's co-occurrence panel. Also fixes it never being reachable as ss.significance_heatmap: it was implemented, tested (against selectsim.plotting directly), and listed in plotting.py's own __all__, but never imported/re-exported from the top-level selectsim package. Pulls the CO/ME/NS color convention (previously only inline in obs_exp_scatter) out into a shared SIG_COLORS constant so both plots stay in sync. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018g7ADBXHN2yWDkf2NPyNm3
selectX gained store/store_path/memory_budget_gb parameters mirroring null_model_parallel's, plus a new "auto" option: estimates the null model's memory footprint from n_genes * n_samples * n_permut (with an empirically-grounded 5x safety factor covering the null-model list itself, in-flight residual/template matrices, and joblib's inter-process copies) and picks "memory" or "zarr" against a detected (os.sysconf) or caller-supplied RAM budget. Lets one selectX() call scale safely across machines of different sizes without the caller needing to know cohort dimensions in advance. Outlier removal on a zarr-backed null model now streams a filtered copy on disk (via null_model._filter_zarr_null_store) instead of materializing retained permutations into a Python list, keeping peak memory O(1 permutation) end to end. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018g7ADBXHN2yWDkf2NPyNm3
…to 0.3.0 The introduction tutorial only demonstrated obs_exp_scatter and ridge_plot_ed. Added sections (executed against the bundled LUAD data) for significance_heatmap, oncoprint, oncoprint_pair, and ridge_plot_ed_compare (the last using an illustrative second selectX run since the bundled data is a single cohort -- real usage would compare two cohorts/timepoints). Bumps version to 0.3.0 for this round's jax-backend removal (a breaking change to null_model_parallel's public API), estimate_p_val's parameter rename, and the new features. NEWS.md has the full changelog. CITATION.cff's version/date-released are left as-is -- those describe an actual tagged release, not yet cut. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018g7ADBXHN2yWDkf2NPyNm3
ℹ️ Recent review info⚙️ Run configurationConfiguration used: defaults Review profile: CHILL Plan: Free Run ID: ⛔ Files ignored due to path filters (18)
📒 Files selected for processing (42)
📝 WalkthroughWalkthroughThe release migrates packaging to Hatchling, adds the SelectSim analysis pipeline with MAF/GAM processing, null-model storage, statistics, plotting, documentation, CI workflows, benchmark scripts, project metadata, and broad automated validation. ChangesSelectSim release
Estimated code review effort: 5 (Critical) | ~120 minutes Note 🎁 Summarized by CodeRabbit FreeYour organization is on the Free plan. CodeRabbit will generate a high-level summary and a walkthrough for each pull request. For a comprehensive line-by-line review, please upgrade your subscription to CodeRabbit Pro by visiting https://app.coderabbit.ai/login. Comment |
Summary
selectX()pipeline,AlterationLandscape, template/weightgeneration, permutation null model, and all overlap/effect-size/, and FDR statistics.
selectsim.gam) for building Gene Alteration Matrices from raw MAF files.obs_exp_scatter,ridge_plot_ed(_compare),significance_heatmap,and two new visualizations not in the R package,
oncoprint/oncoprint_pair.fixed a Zarr chunk-size mismatch (~2.5x faster writes), and removed the JAX backend
after benchmarking it 5-10x slower than NumPy on real data (CPU-only, no GPU to
amortize its dispatch overhead).
store=" auto": picks an in-memory vs. on-disk (zarr) null model automatically basedOn estimated size vs. available RAM, so large cohorts don't need manual tuning.
TCGA LUAD example data, and an R-parity test suite validating output against a real
SelectSim::selectX()R run.See
NEWS.mdfor the detailed, dated changelog.Test plan
pytest), including R-parity numerical validationmake htmlindocs/)Summary by CodeRabbit
New Features
Documentation
Chores