From f93fe2432701d76b263c324099cc5d6df91cb98d Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Mon, 21 Sep 2026 10:23:02 -0400 Subject: [PATCH 01/83] Plan the pysceptre backend for the high-MOI power analysis Writes down how the power simulation moves off R/sceptre and onto pysceptre: the per-target replicate-stacking that makes one call cover a whole replicate chunk, the pseudo-target keys that keep each replicate's CRT draws independent, and the two export gaps in pysceptre that block it. Records the consequences that are easy to get wrong later -- the per-gene null model becomes the faithful refit for free, so the comparison baseline is power_null_fit and not power_as_is -- and gates the work on a benchmark that can stop it. Co-Authored-By: Claude Opus 5 (1M context) --- docs/pysceptre-backend.md | 395 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 395 insertions(+) create mode 100644 docs/pysceptre-backend.md diff --git a/docs/pysceptre-backend.md b/docs/pysceptre-backend.md new file mode 100644 index 0000000..164393b --- /dev/null +++ b/docs/pysceptre-backend.md @@ -0,0 +1,395 @@ +--- +title: Plan - pysceptre backend +nav_order: 8 +--- + +# Plan: move the high-MOI power analysis onto pysceptre + +Replace the R/sceptre engine inside the power simulation with +[pysceptre](https://github.com/broadinstitute/pysceptre), so that the whole pipeline is one Python +package, `pixi.toml` drops R entirely, and the per-simulation cost falls by the margin pysceptre +already demonstrates on discovery analysis. + +This document is a plan, not a record of work done. Nothing below has been implemented. + +## 0. Branches, and why + +| Repo | Branch | Rule | +|---|---|---| +| `WattEG` | **`feat/pysceptre-backend`** (created for this work) | `main` stays the R implementation. `WattEG-paper` was written against it and every number in that paper has to remain reproducible from `main` without archaeology | +| `pysceptre` | **`feature/watteg-simulation-support`** in `../pysceptre` (to create; matching its existing `feature/` prefix) | `pysceptre-paper` freezes pysceptre `main` for the same reason. The changes in §5 land there and merge only once the equivalence checks in §8 pass | +| `WattEG` | `legacy` | untouched (the Snakemake implementation) | + +The R path is not deleted when the Python path lands. It is the reference the Python path is +measured against, and it is what the paper describes; retiring it is a separate decision, taken +after §8 reports. + +## 1. What actually moves + +The pipeline is five steps plus two support steps. Only one of them is expensive, and only one of +them needs R. + +| Today | Fate under the Python backend | +|---|---| +| `src/prepare_sim_input.R` | **ported to Python**, reading a pysceptre export instead of the `.rds` (§4). Everything it computes — poscounts size factors, normalised means, dispersions, the discovery threshold — is computable from the exported matrix | +| `src/split_pairs.R` | ported, trivially | +| `src/fit_null_models.R` | **deleted** (§3.1) | +| `src/merge_null_models.R` | **deleted** (§3.1) | +| `src/run_power_simulation.R` + `lib/simulate.R` + `lib/pert_input.R` | **ported to Python**; the NB draw, guide-to-guide variability, centering and seeding are WattEG's method and stay in WattEG (§6) | +| `src/consolidate_replicates.R` | ported (pyarrow) | +| `src/compute_power.R` + `lib/stats.R` | ported (Wilson interval) | +| `src/summarize_power.R`, `src/fit_power_curve.R` | ported | +| `patches/`, `lib/apply_patch.R`, `src/check_sceptre_api.R`, `src/install_sceptre.R`, `src/install_ondisc.R`, `src/audit_dependencies.R` | **deleted.** All six exist only because the pipeline reaches into unexported sceptre S4 slots and patches its CRT path | +| `lib/sceptre_io.R` | **deleted from the pipeline**; its odm-materialisation logic moves to the one-off export (§4) | +| `src/make_test_data.R` | ported, or replaced by a synthetic fixture generated in Python | + +Everything in `workflow/slurm_executor/` and `workflow/compare_*.R` is scaffolding around the R +steps and follows them. + +## 2. The core mapping, and the reason it is fast + +Today the unit of work is **one `run_discovery_analysis()` call per (target, replicate)**: on +`day0_grna20` that is 3,026 targets x 100 simulations = 302,600 R calls per effect size, each one +carrying the full 567,690-cell bookkeeping for a median of 9 gene pairs. The measured cost model is +`1.140s + 0.5561s x pairs` per (target, simulation), **48 % of it in the per-target term**. + +pysceptre's `run_discovery_analysis` takes plain arrays — `response_matrix` (n_genes x n_cells), +`covariate_matrix`, `grna_target_cells` (dict target -> 0-based cell indices), and a `pairs` frame. +That signature lets one call cover **a whole replicate chunk of one target**: + +``` +for target in split: + rows = [f"{gene}@{rep}" for rep in reps for gene in target_genes] # pseudo-genes + counts = vstack(simulate(target, gene, rep) for rep in reps) # (genes*reps, n_cells) + pairs = DataFrame(response_id=rows, grna_target=[f"{target}@{rep}" ...]) # pseudo-targets + run_discovery_analysis(counts, rows, covariates, {f"{target}@{rep}": cells}, pairs, ...) +``` + +Five consequences, each verified against the pysceptre source rather than assumed: + +**2.1 Per-gene null fits come for free and are per-replicate.** `fit_all_genes` fits every row of +the response matrix independently (`_GENE_BATCH_WIDTH = 1`, `pipeline/discovery.py:89` — pinned at 1 +deliberately so a fit cannot depend on its batch companions). Stacking replicates as rows therefore +fits each `(gene, replicate)` null model **on that replicate's own simulated counts**. That is the +`cleared` configuration in [Status]({{ site.baseurl }}{% link status.md %}) — the faithful reference +that `null_fit` was built to approximate — obtained at no extra cost. §3.1. + +**2.2 Replicates must not share a CRT draw.** pysceptre seeds each target's resampling stream from +the *target's name* (`target_seed_sequence`), which is what makes its results independent of chunking +and target order. One shared key for all replicates would hand every replicate the **same** synthetic +treated-cell index sets, correlating the replicates: R redraws them per call, and the Wilson interval +in `compute_power` assumes independent Bernoulli trials. Hence the `target@rep` pseudo-target keys +above — independent streams by construction, and chunk-layout invariance inherited for free. + +The key must carry the **effect size** as well as the replicate — `target@es@rep` — or a single +`seed` would hand every effect size in a sweep the same synthetic index sets. R's `derive_seed` +includes `effect_size` for exactly this reason. (The alternative, deriving the per-call `seed` from +the effect size, works too but makes the invariance harder to see; prefer the key.) + +The cost is that the target's binomial GLM (perturbation status on covariates) is refit once per +replicate although the input is identical — precisely the redundancy the sceptre patch in `patches/` +removes on the R side, worth 999 -> 634 CPU-h there. pysceptre batches those fits across a target +chunk, so the redundancy is far cheaper than R's. **It is accepted, not fixed** — removing it would +mean changing the pysceptre package, and §5.3 says why that bar is not met. + +**2.3 Memory is a storage question, not a fitting question.** Both row-reading paths — `_gene_fit_job` +(`discovery.py:257`, which at `_GENE_BATCH_WIDTH = 1` fills a one-row float64 `Y`) and `_gene_job` +(`discovery.py:804`) — go through `_get_row` (`discovery.py:190`), which densifies **one row at a +time** and casts to float there. So the stacked matrix can be stored as `int16` or sparse without +any change to pysceptre. For the median target: 9 genes x 100 replicates x +567,690 cells is 1.0 GB as dense `int16`, against 8.2 GB as float64. The max target (36 pairs) is +4.1 GB, which is why `reps_per_chunk` stays a parameter and becomes the memory knob it always +implicitly was. Overflow guard: counts above `int16` range must promote, not wrap. + +**2.4 Call the inner entry point, not the public one.** `pipeline.api.run_discovery_analysis` +sizes `B2`/`B3` from `len(pairs)` (R's `run_qc` rule). With pseudo-pairs the pair count is inflated +by the replicate count, which under `resampling_approximation = "no_approximation"` would inflate +the `B3` budget ~100x. Call `pipeline.discovery.run_discovery_ntcells_complement` directly and pass +`B1`/`B2`/`B3` from the export's `metadata.json` — which is what R does, since the template carries +the real analysis's `@B1/@B2/@B3` slots and the simulation never re-derives them. + +**2.5 Where the time will actually go — and why "30 s" is not the answer.** A full +permutation-mode pysceptre discovery analysis on a real dataset runs in ~30 s. That number does not +transfer to this pipeline, and reading it as though it does would set the wrong expectations and +optimise the wrong thing. One discovery analysis fits **237 gene nulls** and tests ~35,000 pairs +once. One effect size of a 100-replicate power sweep fits **34,886 x 100 ≈ 3.5 million** per-(gene, +replicate) nulls, draws 3.5 million negative-binomial count vectors of 567,690 values, and runs 3.5 +million resampling tests. It is roughly four orders of magnitude more gene-fitting work than the +analysis whose runtime is 30 s. + +So the useful reading of that 30 s is a **per-unit rate**, and the three terms it decomposes into +are what the phase-3 benchmark has to report separately: + +1. drawing the counts (`rnbinom` was 7 % of R's time; in Python it will be a larger share, because + everything around it got faster); +2. the per-(gene, replicate) Poisson IRLS null fit — 59 % of R's time was `glm.fit`, and §2.1 makes + this term unavoidable rather than hoistable; +3. the per-pair resampling test against that target's CRT draws. + +Whichever dominates is where any later optimisation goes. The redundant binomial fits of §2.2 are a +fourth term and, on this reasoning, the smallest of the four — which is the quantitative case for +§5.3 staying unbuilt. + +## 3. Statistical decisions this forces, stated up front + +### 3.1 The null model becomes the faithful refit + +`fit_null_models.R` / `merge_null_models.R` and the whole `--null-precomputations` bundle exist for +one reason: R refitting the per-gene null inside every call costs **4.3x**, so the fits were hoisted +into their own Nextflow processes, fitted once per `(gene, simulation)` on a null simulation, and +injected. Status records the verdict: `null_fit` is **exactly equivalent** to the faithful +`cleared` refit (0 flips in 265 calls), and `as_is` — the inherited real-data cache — understates +power by +0.0063 mean. + +Under §2.1 the faithful refit is what happens anyway. So: delete both processes, delete the seed- +matching guard between bundle and simulation, delete the `@response_precomputations` trap in +`slim_sceptre_object`. **Expect the Python results to match `power_null_fit/`, not `power_as_is/`** +— that is the comparison baseline in §8, and picking the wrong one would manufacture a +0.006 +discrepancy out of a known, already-settled difference. + +### 3.2 Dispersions still come from real data + +`row_data$dispersion` is `1/theta` from sceptre's `@response_precomputations`, fitted on the **real** +counts — it sets the noise the simulation exists to reproduce, and it must not become a property of +the simulated counts. pysceptre has its own validated theta estimator (`glm/nb_theta.py`), so the +Python `prepare_sim_input` computes theta from the real matrix rather than reading a cache. +**Acceptance gate:** per-gene `theta` from `nb_theta.py` against `@response_precomputations$theta` +on the fixture object, reported as a distribution, before anything downstream is trusted. +`build_dispersion_vector`'s hard error on missing/non-finite dispersions carries over. + +### 3.3 Seeding contract is preserved + +Today: `set.seed(derive_seed(seed, target, rep, effect_size))` before each replicate, and a separate +`rep = 0` key for per-target setup, which is what makes results invariant to split layout and +replicate chunking. Python equivalent: `np.random.SeedSequence` spawned from the same four-part key +(hashed, not R's integer arithmetic — the streams differ from R's either way). The invariance test +(1x4 against 2x2 chunking, and two different split layouts) ports directly and becomes a unit test +rather than an sbatch script. + +### 3.4 The `log_2_fold_change < 0` condition + +pysceptre returns `fold_change`, not `log_2_fold_change`; `compute_power`'s one-sided condition +becomes `fold_change < 1`, which is the same predicate. `pct_change_es` and its CI are a bonus the R +path never had. + +## 4. The input boundary — the one real design decision + +`prepare_sim_input.R` reads a `sceptre_object` `.rds`. Nothing in Python reads that. Two options: + +| | Keep `prepare_sim_input.R` in R | **Recommended: `.h5mu` in, R export run once, outside the pipeline** | +|---|---|---| +| Env | still needs `r-base` + pinned sceptre + ondisc + the patch machinery | pure Python; `pixi.toml` drops R | +| User cost | none | one `Rscript export_sceptre_dataset.R` per dataset, in a container we provide | +| Pipeline surface | unchanged | samplesheet column becomes `dataset` (`.h5mu`) instead of `sceptre_object` | + +Take the second. The stated goal is one Python package, and a single R step in the pipeline keeps the +entire R toolchain — pin, patch, `check-api`, two source installs — alive to serve it. The export is +a per-dataset, one-off, already-written script (`pysceptre/scripts/export_sceptre_dataset.R` + +`make_h5mu.py`) that handles odm-backed and in-memory matrices alike, so `lib/sceptre_io.R`'s +materialisation logic has an owner. + +**Keep a `--sceptre-object` path on the R side for one release** as an escape hatch and so §8 can run +both engines from the same object. + +What the export already carries and WattEG needs: the count matrix (`--all-genes`, required — +poscounts size factors are a whole-gene reduction), the covariate matrix restricted to +`cells_in_use`, QC-passing pairs, `discovery_result` (from which the nominal threshold is derived +exactly as `discovery_threshold()` does today), and in `metadata.json` the `side_code`, +`resampling_approximation`, `run_permutations`, `control_group_complement`, `B1/B2/B3`, +`multiple_testing_alpha` and both `n_nonzero_*` thresholds. That covers `analysis_mode` and every +pysceptre argument. + +What it does **not** carry is enough cells — see §5.2, which is a blocking gap, not a detail. + +## 5. What pysceptre must change — the export only + +**The constraint that shapes this whole section: the pysceptre *package* does not change.** Both +gaps below are in `scripts/`, which pysceptre's own `CLAUDE.md` marks as *not shipped in the wheel* — +they are dataset-export tooling, not the statistical engine. Nothing in `src/pysceptre/` is touched, +so nothing this plan does can move a pysceptre result, and the validation burden stays on WattEG +where it belongs. + +Verified gaps, not speculation. + +**5.1 Individual targeting-gRNA assignments — blocking.** The export writes +`grna_assignments$grna_group_idxs`, which is the **union of each target's gRNAs**, plus individual +*non-targeting* gRNAs (`scripts/sceptre_export_lib.R:112-139`). WattEG's guide-to-guide variability +(`create_guide_pert_status`, `create_effect_size_matrix`, `guide_sd = 0.13`) needs **per-gRNA** +membership for targeting guides, and the `grna_id -> grna_target` map. Add both to the export: a +third block of assignment rows with `unit_kind = "targeting_grna"`, and `grna_target_data_frame` +written out whole. + +While the export is open: it writes only `response_id` and `grna_target` for the QC-passing pairs +(`qc_passing_pairs`), but the R simulation output carries `n_nonzero_trt`, `n_nonzero_cntrl` and +`pass_qc` from `@discovery_pairs_with_info` — real-data diagnostics, constant across replicates, and +the first thing anyone looks at when a pair's power is surprising. Write that frame whole, or drop +those three columns from the byte-compatibility promise in §6. Prefer writing it. + +All additive; no existing reader changes. + +**5.2 All cells, not just `cells_in_use` — blocking.** `prepare_sim_input.R:263-270` calls +`compute_expression_stats()` on `get_response_matrix(so)`, the **whole** matrix: 586,309 columns on +`day0_grna20`, against 567,690 in `cells_in_use`. Poscounts size factors are a per-cell reduction +over a per-gene geometric mean taken across *all* cells, so computing them on the QC-passing subset +gives different size factors and different normalised means — which is the input the simulation +draws from. The export writes `cells_in_use` only (`sceptre_export_lib.R:37,61,85`). Add an +`--all-cells` mode that writes the full matrix plus a `cells_in_use` index vector. + +Until that lands, phase 2's column-by-column gate **will** fail, and it would be easy to +misattribute the failure to theta (§3.2). It also constrains Stage A: matrices dumped from Python +have to be indexed the way `template@cells_in_use` expects before R can test them. + +**5.3 A shared target fit across aliased targets — considered and NOT planned.** The `target@es@rep` +keys of §2.2 make pysceptre refit each target's binomial GLM once per replicate although the input +is identical. On the R side removing that redundancy was worth 999 -> 634 CPU-h, which is why it +gets a mention at all. Here it does not: pysceptre batches those fits across a target chunk, and a +full permutation-mode discovery analysis on a real dataset runs in **~30 s**, so the engine is not +plausibly the bottleneck in a simulation whose per-replicate cost is dominated by drawing counts and +fitting per-gene nulls (§2.5). + +Adopting it would mean an API change inside `src/pysceptre/` — an optional `target_fit_key` letting +several target keys share one fit while keeping their own CRT stream. That is a change to how +pysceptre works, so the bar is not "it would be faster": it is the phase-3 benchmark showing the +redundant fits are a **large** share of simulation wall clock. Absent that number, this stays +unbuilt, and the plan assumes it never gets built. + +**5.4 Nothing else.** The engine is used as published. If a change to `discovery.py` turns out to be +needed, that is a signal the mapping in §2 is wrong, not that pysceptre needs a WattEG-shaped hole in +it. + +## 6. Where the code lives + +The simulation model — NB draw from `mean x size_factor x effect_size`, per-guide effect sizes, +re-centering, the seeding scheme — is **WattEG's method**, not part of sceptre, and does not go into +pysceptre. New package in this repo: + +``` +src/watteg/ + sim_input.py # the container, ported from lib/sim_input.R + expression.py # poscounts size factors, normalised means, theta + perturbation.py # pert_input, guide status, effect-size matrix (lib/pert_input.R, lib/simulate.R) + simulate.py # draw_counts + engine.py # the pysceptre call: pseudo-gene/pseudo-target assembly + power.py # Wilson interval, power, MDES (lib/stats.R, compute_power.R) + seeds.py # derive_seed / SeedSequence + cli/ # one entry point per pipeline step +``` + +`pyproject.toml` with console scripts, so the Nextflow modules call `watteg-prepare-sim-input` etc. +rather than `Rscript src/...`. Output file names, columns and TSV/Parquet layouts stay **byte- +compatible** with the R path wherever they can — `consolidate_replicates`, `compute_power` and +`summarize_power` outputs are what the paper's figures read. + +## 7. The DAG afterwards + +``` +samplesheet -> PREPARE_SIM_INPUT -> SPLIT_PAIRS -> POWER_SIMULATION (split x effect size x rep chunk) + -> CONSOLIDATE_REPLICATES -> COMPUTE_POWER -> SUMMARIZE_POWER +``` + +Eight processes become six; `FIT_NULL_MODELS` and `MERGE_NULL_MODELS` go, and with them +`reps_per_null_chunk`, `test_max_null_reps`, the divisibility check on them, and one join in +`main.nf`. `reps_per_chunk` stays and becomes load-bearing for memory (§2.3). + +## 8. Validation — what would make this believable + +R and Python cannot agree draw for draw: different RNGs, different CRT streams, and §3.1 changes the +null model relative to what the R sweeps ran. Validate in stages, against the reference outputs +already in `WattEG-paper` rather than re-running R. + +**Which reference is which** — checked, because getting it backwards manufactures a discrepancy out +of a settled difference. `power_sweep/` holds the **six-effect-size sweep** (`power_es0.05` through +`power_es0.5`) and its `prepared/` carries `null_precomputations.rds`, so it ran the `null_fit` +configuration that §3.1 reproduces. `power_sweep_null/` holds `power_es0.0.tsv` only: it is the +**es = 0 null arm**, not the `null_fit` configuration. Stage B reads `power_sweep/`; Stage C reads +`power_sweep_null/`. + +**Stage 0 — measure the noise floor first.** The 0/265 flips `null_fit` achieved against `cleared` +were possible only because both ran the *same* R RNG stream: identical CRT index sets, differing +only in the null coefficients. R sceptre against pysceptre is two **independent** resampling streams +(4,999 draws plus a skew-normal fit each), so near-threshold p-values will cross the threshold +exactly as two R runs with different seeds do. Re-run the R panel with a second seed and record its +own flip rate and Δpower spread. That is the bar. Every criterion below is expressed against it +rather than against a number picked in advance. + +**Stage A — isolate the engine.** Simulate counts in Python, dump the matrices (indexed to match +`template@cells_in_use`, §5.2), and test the *same* matrices through both engines — R sceptre with +`@response_precomputations` cleared, and pysceptre. With the simulation held fixed, any difference is +the engine and its resampling stream. +- Report: max |Δp|, median |Δp|, Spearman, and the **threshold-flip count** at the discovery + threshold — the statistic `threshold_check.R` already produces. +- Acceptance: flip rate and |Δp| spread **inside the stage-0 floor** on the same 3-target x 53-pair + panel, with no directional bias in the flips. R's 7–0 against `as_is` is what bias looks like; a + 4–3 split is not. + +**Stage B — end to end, independent runs.** Full 100 simulations at effect size 0.15, Python against +`power_sweep/.../power_es0.15.tsv`. +- Report: per-pair Δpower distribution, the fraction exceeding each pair's Wilson half-width, the + mean shift, and the count crossing the 0.8 line in each direction. +- Acceptance: mean shift consistent with zero — two independent 100-draw estimates of the same + binomial p differ by about `sqrt(2p(1-p)/100)`, so ≈0.07 per pair at p = 0.5 and ≈0 in the mean + over 34,886 pairs — and a **symmetric** count of 0.8-line crossings. The 12.6 % of pairs status.md + calls ambiguous will move in both directions; that is expected, not a failure. What would not be + expected is a mean shift, or crossings running one way, which is precisely the signature `as_is` + showed (11,649 up against 3,631 down). + +**Stage C — the null arm.** Effect size 0 against `power_sweep_null/`: the empirical rejection rate +should sit at the nominal threshold in both. This is the pipeline's own calibration check, and the +one stage with an absolute bar rather than a relative one. + +**Stage D — invariance.** 1x4 against 2x2 replicate chunking, and two split layouts, byte-identical +(§3.3). Unit test, not a cluster job. + +## 9. Out of scope, said explicitly + +- **Low-MOI screens.** pysceptre covers the complement-control-group + CRT high-MOI path only. A + low-MOI object must fail at `prepare_sim_input` with a clear message naming the R path, not + silently produce numbers from the wrong control group. The check reads `control_group_complement` + and `run_permutations` out of the export metadata. +- **`--n-control-cells` / `--cell-batches`.** Measured to cost 21–60 % of power and off by default; + recommend not porting, and deleting the flags rather than carrying dead paths. Your call — say so + if they should survive. +- **`run_permutations = TRUE` screens.** pysceptre supports permutations, but its draws are sized by + the largest target in the run, which interacts badly with per-target calls. Refuse for now. + +## 10. Infrastructure + +- `pixi.toml`: drop `r-base`, `r-optparse`, `r-matrix`, `r-rcpp`, `r-dplyr`, `r-data.table`, + `r-purrr`, `r-crayon`, `r-parallelly`, `r-withr`, `r-nanoparquet`, `SCEPTRE_REF`/`SCEPTRE_SHA`, + `ONDISC_REF`/`ONDISC_SHA`, and the `setup` / `check-api` tasks. Add `python`, `numpy`, `scipy`, + `pandas`, `pyarrow`, `numba`, `mudata`. Keep `nextflow`. +- **pysceptre is still a pin.** It is on neither conda-forge nor bioconda, so it enters as a pixi + `[pypi-dependencies]` git dependency pinned by commit — the same shape as `SCEPTRE_SHA`, minus the + patch and the API check. `pixi.toml`'s comment block explaining why sceptre is pinned gets + rewritten, not deleted, and `check_sceptre_api.R`'s job — assert the pin still matches what we + call — passes to pysceptre's own test suite plus this repo's. +- A new container image for the `gcb` profile. The R image is not reusable. +- **`conf/*.config` needs recalibrating from scratch.** The last ten commits on `main` tuned + `POWER_SIMULATION`'s memory and machine type around R's 2.27 GB median / 3.27 GB max. Python's + footprint is dominated by the stacked count matrix (§2.3) and is a different function of + `reps_per_chunk` and pairs-per-target. Do not carry the closures over; re-measure, then rewrite + them. +- `.githooks/pre-commit` rejects camelCase **R** identifiers; add the Python equivalents (ruff, + matching pysceptre's `ruff.toml`) rather than leaving Python unlinted. + +## 11. Phases, with a gate that can stop the work + +1. **Export gap** (pysceptre branch): §5.1, plus a round-trip test that the individual targeting-gRNA + unions reproduce `grna_group_idxs` exactly. *Gate: they do.* +2. **`prepare_sim_input` in Python** against the fixture: size factors, normalised means, theta, + threshold, pairs — each compared to the R output column by column. *Gate: §3.2.* +3. **Benchmark before committing to the shape.** 3 targets x 100 simulations on the fixture, timing + (a) one call per (target, replicate), (b) replicates stacked per §2, at 10/25/50/100 replicates + per chunk, with peak RSS — and **broken down into the four terms of §2.5**, not reported as a + single wall clock. *Gate: this sets `reps_per_chunk`, gives the first honest estimate of CPU-h + per effect size against R's 634, and is the only evidence that could reopen §5.3. If it is not + comfortably ahead of R, stop and re-plan rather than porting the remaining four steps.* +4. **`run_power_simulation` in Python** + Stage A validation. +5. **The four cheap steps** (`split_pairs`, `consolidate_replicates`, `compute_power`, + `summarize_power`, `fit_power_curve`) — mechanical, byte-compatible outputs. +6. **Nextflow rewiring**, new container, resource recalibration (§10). +7. **Stage B/C/D validation** at full scale on one effect size. +8. **Docs**: `usage.md`, `methods.md` and `status.md` rewritten for the Python path; the R path + documented as the reference implementation it has become. + +Phases 1–4 are the work; 5–6 are mechanical; 7 is the one that decides whether `main` moves. From 308b40e7564147e7f57c1ab6a58b91bceafd934c Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Mon, 21 Sep 2026 11:00:36 -0400 Subject: [PATCH 02/83] Record what phase 1 found Both export gaps are closed in pysceptre. Two things the work turned up are worth keeping in the plan rather than in a commit message nobody reads again. The gRNA -> target map is many-to-many: 1,673 of day0's guides sit inside overlapping candidate elements, so perturbation.py must read grna_target_data_frame and never the per-unit annotation. Collapsing it would leave 216 of 3,071 targets simulating with an incomplete guide set, silently. And --all-cells now has a measured price: 2.2% median on raw gene means, 0.46% on size factors, 0.36% on the normalised means the simulation draws from. That bounds a question phase 2 should raise once it reproduces R -- whether QC-failed cells belong in that geometric mean at all -- and says the answer is worth about the third decimal of power. Co-Authored-By: Claude Opus 5 (1M context) --- docs/pysceptre-backend.md | 57 +++++++++++++++++++++++++++++++++++++-- 1 file changed, 55 insertions(+), 2 deletions(-) diff --git a/docs/pysceptre-backend.md b/docs/pysceptre-backend.md index 164393b..d5395ed 100644 --- a/docs/pysceptre-backend.md +++ b/docs/pysceptre-backend.md @@ -203,6 +203,11 @@ What it does **not** carry is enough cells — see §5.2, which is a blocking ga ## 5. What pysceptre must change — the export only +**Status: done, on `feature/watteg-simulation-support` in `../pysceptre`** (commits `7e29efa`, +`5117651`). Both gaps are closed and verified against the real day0 object; §5.1 and §5.2 below are +kept as the record of what they were and of what closing them turned up. §5.3 remains unbuilt, as +planned. + **The constraint that shapes this whole section: the pysceptre *package* does not change.** Both gaps below are in `scripts/`, which pysceptre's own `CLAUDE.md` marks as *not shipped in the wheel* — they are dataset-export tooling, not the statistical engine. Nothing in `src/pysceptre/` is touched, @@ -219,6 +224,20 @@ membership for targeting guides, and the `grna_id -> grna_target` map. Add both third block of assignment rows with `unit_kind = "targeting_grna"`, and `grna_target_data_frame` written out whole. +**What closing it turned up, and what it means for §6.** The gRNA -> target map is +**many-to-many**: 1,673 of day0's 43,736 guides sit inside two or three *overlapping* candidate +elements and so belong to two or three targets (45,463 design rows against 43,736 distinct ids). +R handles this without comment — `grna_map$grna_id[grna_map$grna_target == target]` selects by +target, so a shared guide is simply returned for both. Anything that collapses the map to one +target per guide — a `dict`, a `match()`, the per-unit `var` annotation — drops those guides from +every target but one, which on day0 would leave **216 of 3,071 targets simulating with an +incomplete guide set**, silently. + +So `perturbation.py` must read guides-per-target from `grna_target_data_frame`, **never** from the +gRNA assay's `var` annotation, which has one row per unit and therefore records `""` for +a shared guide. The export asserts the union round-trip exhaustively over every target at write +time, which is what caught this; all 3,071 day0 targets reproduce exactly. + While the export is open: it writes only `response_id` and `grna_target` for the QC-passing pairs (`qc_passing_pairs`), but the R simulation output carries `n_nonzero_trt`, `n_nonzero_cntrl` and `pass_qc` from `@discovery_pairs_with_info` — real-data diagnostics, constant across replicates, and @@ -239,6 +258,36 @@ Until that lands, phase 2's column-by-column gate **will** fail, and it would be misattribute the failure to theta (§3.2). It also constrains Stage A: matrices dumped from Python have to be indexed the way `template@cells_in_use` expects before R can test them. +**Closed by `--all-cells`, and this is what it buys.** Measured across the 18,619 cells QC removes +on day0, computed both ways: + +| Quantity | Median shift | Max | +|---|---:|---:| +| Raw gene mean | 2.2 % | 6.2 % | +| poscounts size factor (in-use cells) | 0.46 % | 1.3 % | +| Size-factor-normalised gene mean | 0.36 % | 3.0 % | + +The file keeps **one cell space** — under `--all-cells` the matrix columns, covariate rows and +every gRNA unit are absolute positions together — and `load_export` subsets back to `cells_in_use` +by default, so an analysis reads either file identically and only the simulation passes +`all_cells=True`. Verified on day0: the default read of the `--all-cells` export is identical to +the plain one across all 92,622,239 nonzeros, the covariates, all 46,789 units and the pair table. + +One thing the round-trip assertion forced into the open: **gRNA membership is post-QC in both +spaces.** `@grna_assignments` is built after QC while `@initial_grna_assignment_list` is the +pre-QC input, so a target's guides between them cover cells the target does not — 518 against 493 +on day0's first target. The guides are restricted to `cells_in_use`, which keeps the union +invariant true in every file and costs nothing, since those cells have no covariates and no test +sees them. `--all-cells` therefore adds cells to the **expression side only**. + +**A question for phase 2, not for the port.** Whether cells QC removed *should* enter the per-gene +geometric mean that sets the size factors is a scientific question, and the honest answer is that +R's implementation includes them because it reads the whole matrix, not because anyone chose it. +The port reproduces R first — that is what the phase-2 gate is for — and the table above bounds +what the choice is worth: a 0.36 % shift in the gene mean the simulation draws from, against +effect sizes of 5–50 %, can move power at the third decimal at most. Worth raising once the Python +path reproduces the R one, and not before. + **5.3 A shared target fit across aliased targets — considered and NOT planned.** The `target@es@rep` keys of §2.2 make pysceptre refit each target's binomial GLM once per replicate although the input is identical. On the R side removing that redundancy was worth 999 -> 634 CPU-h, which is why it @@ -374,8 +423,12 @@ one stage with an absolute bar rather than a relative one. ## 11. Phases, with a gate that can stop the work -1. **Export gap** (pysceptre branch): §5.1, plus a round-trip test that the individual targeting-gRNA - unions reproduce `grna_group_idxs` exactly. *Gate: they do.* +1. ~~**Export gap** (pysceptre branch): §5.1 and §5.2, plus a round-trip test that the individual + targeting-gRNA unions reproduce `grna_group_idxs` exactly.~~ **Done** — `7e29efa`, `5117651` on + `feature/watteg-simulation-support`. *Gate passed:* all 3,071 day0 targets reproduce exactly, + asserted at export time rather than in a test that can be skipped; the default read of an + `--all-cells` export is identical to a plain one on the real screen; 230 tests green, with the + export-format contract covered by 10 new ones that need neither R nor a real dataset. 2. **`prepare_sim_input` in Python** against the fixture: size factors, normalised means, theta, threshold, pairs — each compared to the R output column by column. *Gate: §3.2.* 3. **Benchmark before committing to the shape.** 3 targets x 100 simulations on the fixture, timing From c530cc2188ba1686117258c3fbe20ffb0a3b6219 Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Mon, 21 Sep 2026 11:04:04 -0400 Subject: [PATCH 03/83] Call the size of the size-factor question an estimate, not a bound The 0.36% shift in the normalised gene mean says the direction, not the magnitude: as_is moved mean power by +0.0063 from coefficient differences far larger than that, and nobody has run this one. Phase 2 measures it. Co-Authored-By: Claude Opus 5 (1M context) --- docs/pysceptre-backend.md | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/docs/pysceptre-backend.md b/docs/pysceptre-backend.md index d5395ed..3d8f2e8 100644 --- a/docs/pysceptre-backend.md +++ b/docs/pysceptre-backend.md @@ -283,10 +283,11 @@ sees them. `--all-cells` therefore adds cells to the **expression side only**. **A question for phase 2, not for the port.** Whether cells QC removed *should* enter the per-gene geometric mean that sets the size factors is a scientific question, and the honest answer is that R's implementation includes them because it reads the whole matrix, not because anyone chose it. -The port reproduces R first — that is what the phase-2 gate is for — and the table above bounds -what the choice is worth: a 0.36 % shift in the gene mean the simulation draws from, against -effect sizes of 5–50 %, can move power at the third decimal at most. Worth raising once the Python -path reproduces the R one, and not before. +The port reproduces R first — that is what the phase-2 gate is for — and the table above is the +order-of-magnitude argument for deferring it: a 0.36 % shift in the gene mean the simulation draws +from, against effect sizes of 5–50 %. That is an estimate, not a measurement — `as_is` shifted +mean power by +0.0063 from coefficient differences far larger than this, so the direction is right +and the size is not established. Measure it once the Python path reproduces the R one. **5.3 A shared target fit across aliased targets — considered and NOT planned.** The `target@es@rep` keys of §2.2 make pysceptre refit each target's binomial GLM once per replicate although the input From 9fe30fc77f10ffa9e884ad842a0986605c0d85e0 Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Mon, 21 Sep 2026 11:07:06 -0400 Subject: [PATCH 04/83] Note that the real-data regression passed too test_day0_regression, 6/6 against a re-export of day0 with the new export. The field-level comparison already showed the engine's whole input surface was unchanged; this runs the engine over it and agrees. Co-Authored-By: Claude Opus 5 (1M context) --- docs/pysceptre-backend.md | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/docs/pysceptre-backend.md b/docs/pysceptre-backend.md index 3d8f2e8..f8d73cc 100644 --- a/docs/pysceptre-backend.md +++ b/docs/pysceptre-backend.md @@ -429,7 +429,9 @@ one stage with an absolute bar rather than a relative one. `feature/watteg-simulation-support`. *Gate passed:* all 3,071 day0 targets reproduce exactly, asserted at export time rather than in a test that can be skipped; the default read of an `--all-cells` export is identical to a plain one on the real screen; 230 tests green, with the - export-format contract covered by 10 new ones that need neither R nor a real dataset. + export-format contract covered by 10 new ones that need neither R nor a real dataset; and + `test_day0_regression` passes 6/6 against a re-export of day0 (4 min, 34,886 pairs), so the + export changes move nothing the engine reads. 2. **`prepare_sim_input` in Python** against the fixture: size factors, normalised means, theta, threshold, pairs — each compared to the R output column by column. *Gate: §3.2.* 3. **Benchmark before committing to the shape.** 3 targets x 100 simulations on the fixture, timing From deee55ffcb8d8a5255cfce165af6f1059901fb10 Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Mon, 21 Sep 2026 12:14:15 -0400 Subject: [PATCH 05/83] Point the plan at the squashed commit The export work was squash-merged into pysceptre's 0.1.1rc as d96d48f, so the three commit hashes the plan cited now resolve only on the feature branch. d96d48f is the reference that will still mean something later. Co-Authored-By: Claude Opus 5 (1M context) --- docs/pysceptre-backend.md | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/docs/pysceptre-backend.md b/docs/pysceptre-backend.md index f8d73cc..33c9fa7 100644 --- a/docs/pysceptre-backend.md +++ b/docs/pysceptre-backend.md @@ -17,7 +17,7 @@ This document is a plan, not a record of work done. Nothing below has been imple | Repo | Branch | Rule | |---|---|---| | `WattEG` | **`feat/pysceptre-backend`** (created for this work) | `main` stays the R implementation. `WattEG-paper` was written against it and every number in that paper has to remain reproducible from `main` without archaeology | -| `pysceptre` | **`feature/watteg-simulation-support`** in `../pysceptre` (to create; matching its existing `feature/` prefix) | `pysceptre-paper` freezes pysceptre `main` for the same reason. The changes in §5 land there and merge only once the equivalence checks in §8 pass | +| `pysceptre` | **`0.1.1rc`** in `../pysceptre` | `pysceptre-paper` freezes pysceptre `main` for the same reason. The §5 export work was done on `feature/watteg-simulation-support` and squash-merged into `0.1.1rc` as **`d96d48f`**, alongside the analytical-power work; `main` is unchanged | | `WattEG` | `legacy` | untouched (the Snakemake implementation) | The R path is not deleted when the Python path lands. It is the reference the Python path is @@ -203,8 +203,11 @@ What it does **not** carry is enough cells — see §5.2, which is a blocking ga ## 5. What pysceptre must change — the export only -**Status: done, on `feature/watteg-simulation-support` in `../pysceptre`** (commits `7e29efa`, -`5117651`). Both gaps are closed and verified against the real day0 object; §5.1 and §5.2 below are +**Status: done, in `../pysceptre` on `0.1.1rc`** — squash-merged as **`d96d48f`** ("ADD analytical +per-pair power, per-gRNA exports, and a minimal container"), which carries the export work together +with changes of its own. The three commits it squashes (`7e29efa`, `5117651`, `37f9d81`) survive +only on `feature/watteg-simulation-support`, so `d96d48f` is the reference that will keep +resolving. Both gaps are closed and verified against the real day0 object; §5.1 and §5.2 below are kept as the record of what they were and of what closing them turned up. §5.3 remains unbuilt, as planned. @@ -425,8 +428,8 @@ one stage with an absolute bar rather than a relative one. ## 11. Phases, with a gate that can stop the work 1. ~~**Export gap** (pysceptre branch): §5.1 and §5.2, plus a round-trip test that the individual - targeting-gRNA unions reproduce `grna_group_idxs` exactly.~~ **Done** — `7e29efa`, `5117651` on - `feature/watteg-simulation-support`. *Gate passed:* all 3,071 day0 targets reproduce exactly, + targeting-gRNA unions reproduce `grna_group_idxs` exactly.~~ **Done** — squashed into + `d96d48f` on `0.1.1rc`. *Gate passed:* all 3,071 day0 targets reproduce exactly, asserted at export time rather than in a test that can be skipped; the default read of an `--all-cells` export is identical to a plain one on the real screen; 230 tests green, with the export-format contract covered by 10 new ones that need neither R nor a real dataset; and From 3ddbfc47aca78b001bd25179684cf708957e804f Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Mon, 21 Sep 2026 12:33:49 -0400 Subject: [PATCH 06/83] Port the expression statistics to Python, and gate them against R First of the five steps to move off R. compute_expression_stats() is the only place in the pipeline that reads real expression values, so it is where a silent numerical difference would do the most damage and where the gate belongs. Against an R run on the same machine the port is bit-identical: the per-gene geometric means agree exactly and the per-cell medians to one ulp. Against the published day0 sim_input.rds both differ by at most 6.2e-12, and that residual is the reference's platform rather than the port -- R's sum() accumulates in LDOUBLE, 80 bits on the x86 cluster the sweep ran on and plain double on arm64. A local R run reproduces only 30 of 20,000 published size factors bit for bit and differs by the same 6.2e-12. So 1e-10 is what "reproduces R" can mean across platforms, and the comparison prints a bit-identical count so the distinction stays visible. The gate earned itself immediately. pysceptre's export stores counts as uint16 to keep the file small, and numpy's log of a 16-bit integer array returns float32 -- so every logarithm was being computed in single precision. It moved the size factors by 2.5e-7 and the normalised gene means by 1.2e-9: far too small to notice, far too large to be floating point, and invisible to any tolerance chosen by eye. Also scaffolds the package -- pyproject with a console script, ruff config matching pysceptre's, and pysceptre as an editable path dependency until the pixi environment is rebuilt. Co-Authored-By: Claude Opus 5 (1M context) --- .gitignore | 3 + pyproject.toml | 34 +++++++ ruff.toml | 6 ++ src/watteg/__init__.py | 3 + src/watteg/cli/__init__.py | 0 src/watteg/expression.py | 140 +++++++++++++++++++++++++++ workflow/compare_expression_stats.py | 124 ++++++++++++++++++++++++ workflow/dump_r_sim_input.R | 60 ++++++++++++ 8 files changed, 370 insertions(+) create mode 100644 pyproject.toml create mode 100644 ruff.toml create mode 100644 src/watteg/__init__.py create mode 100644 src/watteg/cli/__init__.py create mode 100644 src/watteg/expression.py create mode 100755 workflow/compare_expression_stats.py create mode 100755 workflow/dump_r_sim_input.R diff --git a/.gitignore b/.gitignore index b499dae..3720113 100644 --- a/.gitignore +++ b/.gitignore @@ -88,3 +88,6 @@ tests/data/ # Created by tests/config/run.sbatch -- sparse files whose SIZE is the fixture, 3.3 GB of nothing. tests/config/inputs/ +.venv/ +__pycache__/ +*.egg-info/ diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..f0596ea --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,34 @@ +[project] +name = "watteg" +version = "0.1.0.dev0" +description = "Power analysis for element-gene pairs in single-cell CRISPR screens" +readme = "README.md" +license = { file = "LICENSE" } +requires-python = ">=3.10" +dependencies = [ + "numpy>=1.24", + "scipy>=1.10", + "pandas>=2.0", + "h5py>=3.8", + "pyarrow>=14", + "pysceptre", +] + +[project.scripts] +watteg-prepare-sim-input = "watteg.cli.prepare_sim_input:main" + +[project.optional-dependencies] +dev = ["pytest>=7", "ruff>=0.5"] + +[build-system] +requires = ["hatchling"] +build-backend = "hatchling.build" + +[tool.hatch.build.targets.wheel] +packages = ["src/watteg"] + +# pysceptre is not on PyPI. Pinned by path while the port is in development; this becomes a +# git dependency pinned by commit when the pixi environment is rebuilt (see +# docs/pysceptre-backend.md section 10). +[tool.uv.sources] +pysceptre = { path = "../pysceptre", editable = true } diff --git a/ruff.toml b/ruff.toml new file mode 100644 index 0000000..6a6648a --- /dev/null +++ b/ruff.toml @@ -0,0 +1,6 @@ +# Matches pysceptre's, so the two halves of the port read the same. +line-length = 100 +target-version = "py310" + +[lint] +select = ["E", "F", "I", "UP", "B", "SIM"] diff --git a/src/watteg/__init__.py b/src/watteg/__init__.py new file mode 100644 index 0000000..29c0b6d --- /dev/null +++ b/src/watteg/__init__.py @@ -0,0 +1,3 @@ +"""Power analysis for element-gene pairs in single-cell CRISPR screens.""" + +__version__ = "0.1.0.dev0" diff --git a/src/watteg/cli/__init__.py b/src/watteg/cli/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/watteg/expression.py b/src/watteg/expression.py new file mode 100644 index 0000000..48879a5 --- /dev/null +++ b/src/watteg/expression.py @@ -0,0 +1,140 @@ +"""Per-gene and per-cell expression statistics: what the simulation draws from. + +A port of `compute_expression_stats()` in `src/prepare_sim_input.R`, and the only +place in the Python pipeline that reads real expression values. Everything +downstream works from the three numbers this produces -- a per-cell size factor, +a per-gene size-factor-normalised mean, and a per-gene raw mean -- plus the +dispersion in `dispersion.py`. + +**Computed over every cell, including the ones QC removed.** DESeq2 "poscounts" +size factors are a per-cell reduction against a per-gene geometric mean taken +across the whole matrix, so which cells are in the matrix changes the size +factors of the cells that stay: measured on day0, a median 0.46 % shift in size +factor and 0.36 % in normalised gene mean. The R implementation reads the +object's whole response matrix, so reproducing it means doing the same, which is +what `--all-cells` on pysceptre's export is for. Whether those cells *should* +count is a separate question -- see `docs/pysceptre-backend.md` section 5.2. + +DESeq2 itself is not used, in either language: `DESeqDataSetFromMatrix()` coerces +to dense, which at 292 x 586,309 is 1.4 GB and pointless. +""" + +from __future__ import annotations + +from dataclasses import dataclass + +import numpy as np +from scipy import sparse + + +@dataclass(frozen=True) +class ExpressionStats: + """Aligned to the matrix that produced them: one entry per gene, per cell.""" + + size_factors: np.ndarray # (n_cells,) + average_expression_all_cells: np.ndarray # (n_genes,) raw mean + normalized_mean: np.ndarray # (n_genes,) what the simulation draws from + density: float + + +def compute_expression_stats(counts: sparse.spmatrix) -> ExpressionStats: + """`counts` is (n_genes, n_cells), sparse, non-negative. + + Never densified. The matrix is worked on in CSC because every statistic here + is a per-cell reduction, and CSC is what makes a cell's nonzeros contiguous. + """ + csc = counts.tocsc(copy=False) + # Explicit zeros would take log(0) = -inf into the geometric mean and turn a + # whole gene unusable. R calls drop0() here for the same reason. + csc.eliminate_zeros() + n_genes, n_cells = csc.shape + if csc.data.size and csc.data.min() < 0: + raise ValueError("count matrix contains negative values") + + # EXPLICIT float64, and not a formality. pysceptre's export stores counts as + # `uint16` to keep the file small, and numpy's `log` of a 16-bit integer + # array returns **float32**: every logarithm below would be computed in + # single precision, silently. Measured against R on day0, that alone moved + # the size factors by 2.5e-7 and the normalised gene means by 1.2e-9 -- + # small enough to look like floating-point noise and large enough not to be. + data = csc.data.astype(np.float64, copy=False) + gene_of_nonzero = csc.indices + log_counts = np.log(data) + + # Geometric mean per gene, on the log scale, over the cells where the gene is + # nonzero but divided by *every* cell -- that is what makes this "poscounts" + # rather than an ordinary geometric mean. + log_geomeans = np.bincount(gene_of_nonzero, weights=log_counts, minlength=n_genes) / n_cells + row_totals = np.bincount(gene_of_nonzero, weights=data, minlength=n_genes) + log_geomeans[row_totals == 0] = -np.inf + usable_gene = np.isfinite(log_geomeans) + + # Each count against its gene's geometric mean. A gene with no usable mean + # contributes nothing to any cell's median rather than contributing a zero, + # which would drag every size factor towards 1. + ratio = log_counts - log_geomeans[gene_of_nonzero] + ratio[~usable_gene[gene_of_nonzero]] = np.nan + + size_factors = np.exp(_median_per_cell(ratio, csc.indptr, n_cells)) + + bad = ~np.isfinite(size_factors) | (size_factors <= 0) + if bad.any(): + raise ValueError( + f"{bad.sum()} of {n_cells} cells have no usable size factor (no nonzero counts in " + "any gene with a finite geometric mean). Filter these cells out before running the " + "power analysis." + ) + + # Raw mean, reported in the output but never simulated from. Kept distinct + # from the normalised mean below, which is what the simulation draws from. + average_expression_all_cells = row_totals / n_cells + + cell_of_nonzero = _cell_of_nonzero(csc.indptr, data.size) + normalized = data / size_factors[cell_of_nonzero] + normalized_mean = np.bincount(gene_of_nonzero, weights=normalized, minlength=n_genes) / n_cells + + return ExpressionStats( + size_factors=size_factors, + average_expression_all_cells=average_expression_all_cells, + normalized_mean=normalized_mean, + density=data.size / (float(n_genes) * n_cells), + ) + + +def _cell_of_nonzero(indptr: np.ndarray, n_nonzero: int) -> np.ndarray: + """Which cell each CSC nonzero belongs to.""" + return np.repeat(np.arange(indptr.size - 1), np.diff(indptr)) + + +def _median_per_cell(ratio: np.ndarray, indptr: np.ndarray, n_cells: int) -> np.ndarray: + """Median of each cell's ratios, NaN where a cell has none. + + One global sort rather than a loop over cells. The nonzeros are already + grouped by cell -- that is what CSC means -- so sorting by `(cell, ratio)` + leaves each cell's values contiguous *and* ordered, and every median is then + an index lookup. At 95.7M nonzeros a Python loop over 586,309 cells is + minutes; this is seconds. + + The even-length case averages the two middle values, which is what both R's + `median.default` and `np.median` do, so the two agree element for element. + """ + out = np.full(n_cells, np.nan) + if ratio.size == 0: + return out + + # NaN sorts last under lexsort, so a cell's usable values stay at the front + # of its slice and the median is taken over the count of usable ones. + cell_of_nonzero = _cell_of_nonzero(indptr, ratio.size) + order = np.lexsort((ratio, cell_of_nonzero)) + sorted_ratio = ratio[order] + + usable_per_cell = np.bincount(cell_of_nonzero[np.isfinite(ratio)], minlength=n_cells).astype( + np.int64 + ) + has_any = usable_per_cell > 0 + starts = indptr[:-1] + lower = starts[has_any] + (usable_per_cell[has_any] - 1) // 2 + upper = starts[has_any] + usable_per_cell[has_any] // 2 + out[has_any] = (sorted_ratio[lower] + sorted_ratio[upper]) / 2.0 + + return out diff --git a/workflow/compare_expression_stats.py b/workflow/compare_expression_stats.py new file mode 100755 index 0000000..5db0d1e --- /dev/null +++ b/workflow/compare_expression_stats.py @@ -0,0 +1,124 @@ +#!/usr/bin/env python3 +"""Phase 2's gate: the Python expression statistics against R's own output. + +R's reference is a `sim_input.rds` dumped by `workflow/dump_r_sim_input.R` -- +for day0, the one the published sweep actually ran on. Python reads the +`--all-cells` pysceptre export of the same sceptre object, so both sides see the +same counts over the same 586,309 cells and any difference is the arithmetic. + + workflow/compare_expression_stats.py + +Three quantities are compared, and `dispersion` is deliberately not one of them: +it comes from a model fit rather than from this arithmetic, and +`workflow/compare_dispersion.py` reports it as a distribution rather than as a +pass or a fail. + +**The tolerance is 1e-10, and the reason is the reference's platform rather +than the port.** Against an R run on the *same machine* this port is +bit-identical: the per-gene geometric means agree exactly and the per-cell +medians to one ulp. Against the published `sim_input.rds` both differ by up to +6.2e-12, because R's `sum()` accumulates in `LDOUBLE` -- 80-bit extended +precision on the x86 cluster the sweep ran on, plain `double` on arm64, where +`.Machine$sizeof.longdouble == 8`. Measured: a local R run reproduces only +30 of 20,000 published size factors bit for bit, and differs by the same 6.2e-12 +this port does. So the published outputs cannot be reproduced exactly by R +either, and 1e-10 is what "reproduces R" can mean across platforms. + +The `bit-identical` column is the part to watch: it should read 100 % against a +same-platform reference, and near 0 % against the published one. A tolerance set +loosely by eye would have hidden the bug this comparison caught -- the export +stores counts as `uint16`, and `np.log` of a 16-bit integer array returns +**float32**, which moved the size factors by 2.5e-7. +""" + +from __future__ import annotations + +import argparse +import sys +from pathlib import Path + +import numpy as np + +from watteg.expression import compute_expression_stats + +RTOL = 1e-10 + +FAILURES: list[str] = [] + + +def read_floats(path: Path) -> np.ndarray: + return np.array([float(x) for x in path.read_text().split()]) + + +def report(label: str, ours: np.ndarray, theirs: np.ndarray, rtol: float = RTOL) -> None: + ours = np.asarray(ours, dtype=float) + theirs = np.asarray(theirs, dtype=float) + if ours.shape != theirs.shape: + print(f"FAIL {label}: shape {ours.shape} vs R's {theirs.shape}") + FAILURES.append(label) + return + rel = np.abs(ours - theirs) / np.maximum(np.abs(theirs), np.finfo(float).tiny) + exact = int((ours == theirs).sum()) + ok = bool(np.allclose(ours, theirs, rtol=rtol, atol=0)) + print( + f"{'PASS' if ok else 'FAIL'} {label:30s} n={ours.size:>9,} " + f"max rel {rel.max():.2e} bit-identical {exact:,}/{ours.size:,}" + ) + if not ok: + FAILURES.append(label) + for i in np.argsort(rel)[-3:][::-1]: + print(f" [{i}] ours {ours[i]!r} R {theirs[i]!r} rel {rel[i]:.3e}") + + +def main() -> int: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("export", type=Path, help="an --all-cells pysceptre export") + parser.add_argument("reference", type=Path, help="output of workflow/dump_r_sim_input.R") + args = parser.parse_args() + + sys.path.insert(0, str(Path(__file__).resolve().parents[2] / "pysceptre" / "scripts")) + from sceptre_io import load_export + + export = load_export(args.export, all_cells=True) + counts = export.response_matrix + print(f"export: {counts.shape[0]} genes x {counts.shape[1]:,} cells, {counts.nnz:,} nonzeros") + if counts.shape[1] == export.metadata.get("n_cells_in_use"): + print( + " WARNING: this export has no QC-failed cells. R computes the statistics over every " + "cell in the object, so a comparison against it needs an --all-cells export." + ) + + stats = compute_expression_stats(counts) + print(f" density {stats.density:.1%}") + + ref_genes = (args.reference / "genes.txt").read_text().split() + print(f"reference: {len(ref_genes)} genes\n") + + gene_index = {g: i for i, g in enumerate(export.gene_ids)} + missing = [g for g in ref_genes if g not in gene_index] + if missing: + print(f"FAIL {len(missing)} of R's genes are absent from the export, e.g. {missing[:3]}") + FAILURES.append("gene alignment") + return 1 + rows = np.array([gene_index[g] for g in ref_genes]) + + report( + "row_data$mean", stats.normalized_mean[rows], read_floats(args.reference / "row_mean.txt") + ) + report( + "row_data$average_expression", + stats.average_expression_all_cells[rows], + read_floats(args.reference / "row_average_expression_all_cells.txt"), + ) + report( + "col_data$size_factors", + stats.size_factors, + read_floats(args.reference / "col_size_factors.txt"), + ) + + print("\n" + ("ALL PASS" if not FAILURES else f"{len(FAILURES)} FAILED: {FAILURES}")) + return 1 if FAILURES else 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/workflow/dump_r_sim_input.R b/workflow/dump_r_sim_input.R new file mode 100755 index 0000000..2a8f5c1 --- /dev/null +++ b/workflow/dump_r_sim_input.R @@ -0,0 +1,60 @@ +#!/usr/bin/env Rscript +# +# Dump an R sim_input.rds to text, at full precision, so the Python port can be compared to it. +# +# The counterpart is workflow/compare_expression_stats.py. Written at 17 significant digits and +# not with write.csv: write.csv gives 15, which is enough to make a bit-identical port look like +# it agrees only to 6e-12 -- measured, and misleading enough to have sent one investigation the +# wrong way. +# +# Usage: +# workflow/dump_r_sim_input.R +# +# The reference used for day0 is the sim_input the published sweep ran on: +# WattEG-paper/power_sweep/day0/day0/prepared/sim_input.rds + +suppressPackageStartupMessages(library(Matrix)) + +args <- commandArgs(trailingOnly = TRUE) +if (length(args) != 2) { + stop("usage: dump_r_sim_input.R ", call. = FALSE) +} +sim_path <- args[[1]] +out <- args[[2]] +dir.create(out, showWarnings = FALSE, recursive = TRUE) + +sim <- readRDS(sim_path) +cat("genes:", length(sim$genes), " cells:", length(sim$cells), "\n") + +# 17 significant digits round-trips a double exactly; formatC is used rather than format() so the +# width is not chosen per vector. +full <- function(x) formatC(as.numeric(x), digits = 17, format = "e") + +writeLines(sim$genes, file.path(out, "genes.txt")) +writeLines(sim$cells, file.path(out, "cells.txt")) + +for (column in colnames(sim$row_data)) { + writeLines(full(sim$row_data[[column]]), file.path(out, paste0("row_", column, ".txt"))) +} +for (column in colnames(sim$col_data)) { + values <- sim$col_data[[column]] + # Factors go out as their labels; the levels are recoverable from the labels and the order is + # not load-bearing on this side of the comparison. + writeLines(if (is.numeric(values)) full(values) else as.character(values), + file.path(out, paste0("col_", column, ".txt"))) +} + +for (nm in names(sim$perts)) { + m <- as(sim$perts[[nm]], "CsparseMatrix") + writeLines(rownames(m), file.path(out, paste0(nm, "_rows.txt"))) + # Triplets, 0-based, so the Python side can rebuild the matrix without knowing R's storage. + triplets <- summary(m) + utils::write.table( + data.frame(i = triplets$i - 1L, j = triplets$j - 1L), + file.path(out, paste0(nm, "_triplets.tsv")), + sep = "\t", row.names = FALSE, quote = FALSE + ) + cat(" ", nm, ":", paste(dim(m), collapse = " x "), "with", nrow(triplets), "nonzeros\n") +} + +cat("wrote the reference to", out, "\n") From 5d63834beaf69140d3511061c5b1d875d1bca222 Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Mon, 21 Sep 2026 12:36:39 -0400 Subject: [PATCH 07/83] Fit the dispersions in Python instead of reading sceptre's cache The simulation draws from NB(mu, size = 1/dispersion), so this number decides how noisy a simulated gene is and therefore how much of an effect the test can see. R read it from sceptre_object@response_precomputations, a cache sceptre fills during the real discovery analysis; pysceptre fits the same model, so the Python pipeline estimates theta itself and depends on no slot of anyone's object. It agrees with sceptre's cache to about ten significant digits: over 237 genes and 567,690 cells, the relative difference is 2.8e-12 at the median and 1.2e-9 at worst, and 236 of 237 genes agree to better than 1e-9. Since a dispersion matters only through the variance of the counts drawn from it, that is nothing next to the 5-50% effect sizes the sweep tests. One gene's MLE fell back to method of moments in both implementations alike. Fitted over cells_in_use, unlike everything in expression.py, because that is what sceptre's precomputation does. The two halves using different cell sets is not an oversight in either language: size factors are a property of the library preparation, which every cell took part in, and a model fit is a property of the analysis, which only QC-passing cells enter. A theta clamped to the estimator's bounds is refused rather than simulated from. R needed no such check because it read a cache someone else had already produced; here the fit happens in front of us. Co-Authored-By: Claude Opus 5 (1M context) --- src/watteg/dispersion.py | 84 ++++++++++++++++++++++++++++++ workflow/compare_dispersion.py | 93 ++++++++++++++++++++++++++++++++++ 2 files changed, 177 insertions(+) create mode 100644 src/watteg/dispersion.py create mode 100755 workflow/compare_dispersion.py diff --git a/src/watteg/dispersion.py b/src/watteg/dispersion.py new file mode 100644 index 0000000..f244312 --- /dev/null +++ b/src/watteg/dispersion.py @@ -0,0 +1,84 @@ +"""Per-gene dispersion: the noise the simulation exists to reproduce. + +The simulation draws counts from a negative binomial with `size = 1/dispersion`, +so this number sets how noisy a simulated gene is, and therefore how much of the +effect a test can see. It has to come from the **real** counts: a dispersion +estimated from simulated data would make the simulation agree with itself. + +In R it was read from `sceptre_object@response_precomputations`, a cache sceptre +fills during the real discovery analysis. There is no such cache here and there +does not need to be -- pysceptre fits the same model, so the Python pipeline +estimates theta directly and depends on no slot of anyone's object. + +**Fitted over `cells_in_use`, unlike everything in `expression.py`.** sceptre's +precomputation runs on the cells that passed QC, against the covariate matrix, +so matching it means doing the same. That the two halves of `prepare_sim_input` +use different cell sets is not an oversight in either language: size factors are +a property of the library preparation, which every cell took part in, while a +model fit is a property of the analysis, which only the QC-passing cells enter. +""" + +from __future__ import annotations + +import numpy as np +from pysceptre.pipeline.discovery import ( + _GENE_BATCH_WIDTH, + _GENE_CHUNK_MEMORY_GB, + fit_all_genes, +) + + +def fit_dispersions( + counts, + gene_ids: list[str], + covariate_matrix: np.ndarray, + genes: list[str], + *, + n_jobs: int = 1, +) -> dict[str, float]: + """Dispersion (`1/theta`) per gene, in the order `genes` gives. + + `counts` is (n_genes, n_cells) over `cells_in_use`, `gene_ids` labels its + rows, and `covariate_matrix` is (n_cells, p) over the same cells. + + The batching arguments are pysceptre's own pipeline defaults rather than + this function's choices: `_GENE_BATCH_WIDTH = 1` makes each fit depend on + that gene alone, which is what keeps a dispersion from shifting when the + gene list changes. + """ + index = {g: i for i, g in enumerate(gene_ids)} + unknown = [g for g in genes if g not in index] + if unknown: + raise KeyError( + f"{len(unknown)} gene(s) are not rows of the count matrix, including: {unknown[:5]}" + ) + + fits = fit_all_genes( + counts, + list(genes), + covariate_matrix, + chunk_memory_gb=_GENE_CHUNK_MEMORY_GB, + gene_rows=[index[g] for g in genes], + batch_width=_GENE_BATCH_WIDTH, + n_jobs=n_jobs, + ) + + # A clamped theta is a fit that did not converge to anything usable, and a + # dispersion of 100 or 0.001 would simulate a gene nothing like the real + # one. R had no equivalent check because it read a cache that sceptre had + # already produced; here the fit happens in front of us, so it is checked. + clamped = [g for g in genes if fits[g].theta_clamped] + if clamped: + raise ValueError( + f"{len(clamped)} gene(s) have a theta clamped to the estimator's bounds, including: " + f"{clamped[:5]}. Their dispersion is not an estimate, and simulating from it would " + "misstate the noise. Drop these genes from the discovery pairs." + ) + nonfinite = [g for g in genes if not np.isfinite(fits[g].theta) or fits[g].theta <= 0] + if nonfinite: + raise ValueError( + f"{len(nonfinite)} gene(s) have a non-finite or non-positive theta, including: " + f"{nonfinite[:5]}." + ) + + return {g: 1.0 / fits[g].theta for g in genes} diff --git a/workflow/compare_dispersion.py b/workflow/compare_dispersion.py new file mode 100755 index 0000000..6a82c16 --- /dev/null +++ b/workflow/compare_dispersion.py @@ -0,0 +1,93 @@ +#!/usr/bin/env python3 +"""Phase 2's other half: Python's fitted theta against sceptre's cached one. + + workflow/compare_dispersion.py + +This is **not** a pass/fail comparison and deliberately has no tolerance. The +other quantities in `prepare_sim_input` are arithmetic over the same numbers, so +they can be required to agree to floating point. A dispersion is the output of +an iterative maximum-likelihood fit: sceptre's cache came from its own Poisson +IRLS and theta estimator, this comes from pysceptre's, and the two agreeing to +three or four digits is the expected result rather than a disappointing one. + +What the output is for is judging whether the difference could move power. +The simulation draws from `NB(mu, size = 1/dispersion)`, so a relative shift in +dispersion is roughly a relative shift in simulated variance -- read the +distribution below against the effect sizes being tested, which are 5-50 %. + +Reported: the distribution over all genes, which genes sit furthest apart, and +the theta estimator each side used, because a gene that fell back to a different +method is the first place to look when one disagrees. +""" + +from __future__ import annotations + +import argparse +import sys +from pathlib import Path + +import numpy as np + +from watteg.dispersion import fit_dispersions + + +def read_floats(path: Path) -> np.ndarray: + return np.array([float(x) for x in path.read_text().split()]) + + +def main() -> int: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("export", type=Path, help="a pysceptre export") + parser.add_argument("reference", type=Path, help="output of workflow/dump_r_sim_input.R") + parser.add_argument("--n-jobs", type=int, default=1) + args = parser.parse_args() + + sys.path.insert(0, str(Path(__file__).resolve().parents[2] / "pysceptre" / "scripts")) + from sceptre_io import load_export + + # The default load, NOT all_cells: sceptre fits its precomputation over the + # QC-passing cells, so matching it means fitting over those and no others. + export = load_export(args.export) + genes = (args.reference / "genes.txt").read_text().split() + print( + f"fitting {len(genes)} genes over {export.covariate_matrix.shape[0]:,} cells_in_use, " + f"{export.covariate_matrix.shape[1]} covariates ..." + ) + + ours = fit_dispersions( + export.response_matrix, + export.gene_ids, + export.covariate_matrix, + genes, + n_jobs=args.n_jobs, + ) + mine = np.array([ours[g] for g in genes]) + theirs = read_floats(args.reference / "row_dispersion.txt") + + rel = np.abs(mine - theirs) / theirs + print(f"\ndispersion, {len(genes)} genes (R's cache vs pysceptre's fit)") + for q in (50, 90, 99, 100): + print(f" {q:>3}th percentile of |relative difference|: {np.percentile(rel, q):.3e}") + for tol in (1e-12, 1e-9, 1e-6, 1e-3): + print(f" within {tol:.0e}: {(rel < tol).sum()}/{rel.size}") + print(f" bit-identical: {(mine == theirs).sum()}/{rel.size}") + + worst = np.argsort(rel)[-5:][::-1] + print("\n furthest apart:") + for i in worst: + print( + f" {genes[i]:<16} ours {mine[i]:.6g} R {theirs[i]:.6g} " + f"(rel {rel[i]:.3e}; theta {1 / mine[i]:.10g} vs {1 / theirs[i]:.10g})" + ) + + # A dispersion difference matters only through the variance of the counts + # drawn from it, so say it in those terms rather than leaving the reader to. + print( + f"\n A gene at the median difference simulates with {np.median(rel):.2e} relative " + "difference in negative-binomial variance against R's." + ) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) From 61152ac4c1ffdf66f8cf76f364cfb7b45babe338 Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Mon, 21 Sep 2026 12:48:07 -0400 Subject: [PATCH 08/83] Port prepare_sim_input to Python, and check it against R's own output The last piece of phase 2: the sim_input container, the CLI, and the end-to-end comparison. The input is a pysceptre --all-cells export rather than a sceptre object, so nothing in this step reads R. Against the day0 sim_input the published sweep ran on: pairs.tsv, grna_targets.tsv byte-identical discovery_threshold.txt byte-identical genes same set, same order cre_perts all 3,071 targets' cell sets agree grna_perts all 43,718 guides' cell sets agree with the per-gene and per-cell numbers covered by the two comparisons that landed earlier. Two differences are structural and are checked rather than waived. R's sim_input spans every cell in the object while this one spans cells_in_use -- the simulation only touches cells that have covariates, and R carried the rest into a matrix sceptre then discarded -- so the perturbation matrices are compared after mapping back to absolute cell positions. And R's grna_perts has one row per row of the design table, so a guide in two overlapping elements appears twice under one name; this keeps one row per guide. The comparison now also fails if R has a unit this does not AND that unit has a QC-passing cell. One guide differs, and it has exactly one assigned cell, which failed QC. That last check caught a real one first: grna_perts must include the NON-TARGETING guides. R builds the matrix from every gRNA, and create_guide_pert_status gives each control cell whichever guide it carries -- usually a non-targeting one, which then draws an effect size around 1 with the same guide-to-guide spread. Leaving them out sent every control cell to the no-effect row, keeping the control arm's mean and losing its variance, with nothing downstream to say so. Co-Authored-By: Claude Opus 5 (1M context) --- src/watteg/cli/prepare_sim_input.py | 256 ++++++++++++++++++++++++++++ src/watteg/sim_input.py | 205 ++++++++++++++++++++++ workflow/compare_sim_input.py | 183 ++++++++++++++++++++ 3 files changed, 644 insertions(+) create mode 100644 src/watteg/cli/prepare_sim_input.py create mode 100644 src/watteg/sim_input.py create mode 100755 workflow/compare_sim_input.py diff --git a/src/watteg/cli/prepare_sim_input.py b/src/watteg/cli/prepare_sim_input.py new file mode 100644 index 0000000..783a388 --- /dev/null +++ b/src/watteg/cli/prepare_sim_input.py @@ -0,0 +1,256 @@ +"""Turn a pysceptre export into the small inputs the power simulation needs. + +The step that makes the rest of the pipeline cheap: every one of the +`n_splits x n_effect_sizes x n_rep_chunks` parallel tasks reads what this writes, +and none of them reads a count matrix. See `watteg/sim_input.py` for why that is +enough. + + watteg-prepare-sim-input --dataset day0.h5mu --outdir prepared/ + +The input is a `.h5mu` written by pysceptre's `scripts/export_sceptre_dataset.R` +plus `make_h5mu.py`, **exported with `--all-cells`**. R read the sceptre object +directly; nothing in Python does, which is what lets the environment be one +Python package. The `--all-cells` requirement is not a detail -- see +`watteg/expression.py`. + +Outputs, matching the R step's names and columns so the two can be compared and +so downstream readers do not care which produced them: + + sim_input.h5 per-gene and per-cell statistics + the perturbation matrices + pairs.tsv the QC-passing discovery pairs + grna_targets.tsv the gRNA -> target mapping + discovery_threshold.txt the nominal p-value a simulated pair has to beat + analysis_mode.tsv which test the screen was run under + +There is no `sceptre_template.rds`: that was an R object carrying the covariate +matrix and the analysis parameters, which now live in `sim_input.h5` and +`analysis_mode.tsv` respectively. +""" + +from __future__ import annotations + +import argparse +import sys +from pathlib import Path + +import numpy as np +import pandas as pd +from scipy import sparse + +from watteg.dispersion import fit_dispersions +from watteg.expression import compute_expression_stats +from watteg.sim_input import SimInput, write_sim_input + + +def _load_export(path: Path): + """pysceptre's loader lives in its `scripts/`, which is not shipped in the + wheel, so it is imported by path rather than by package.""" + import pysceptre + + scripts = Path(pysceptre.__file__).resolve().parents[2] / "scripts" + sys.path.insert(0, str(scripts)) + from sceptre_io import load_export, subset_to_cells_in_use + + return load_export(path, all_cells=True), subset_to_cells_in_use + + +def indicator_matrix( + units: list[str], cells_of: dict[str, np.ndarray], n_cells: int +) -> sparse.csr_matrix: + """A (units x cells) 0/1 matrix from a unit -> cells mapping. + + Built in one allocation from the concatenated indices rather than by + growing a matrix. Values are forced to 1: these are indicators, and a cell + that appears twice under one unit must not count double. + """ + rows = np.concatenate([np.full(cells_of[u].size, i) for i, u in enumerate(units)]) + cols = np.concatenate([cells_of[u] for u in units]) if units else np.empty(0, dtype=np.int64) + matrix = sparse.csr_matrix( + (np.ones(rows.size, dtype=np.int8), (rows, cols)), + shape=(len(units), n_cells), + ) + matrix.data[:] = 1 + return matrix + + +def discovery_threshold(discovery_result: pd.DataFrame | None) -> float: + """The largest p-value the real analysis still called significant. + + That is the bar a simulated pair has to clear, so power means "would this + screen have called it" rather than "would some other threshold have". R + derives it the same way, from the same column. + """ + if discovery_result is None or not len(discovery_result): + raise ValueError( + "the export carries no discovery_result, so no significance threshold can be " + "derived. Re-export with --discovery-result, or pass --threshold explicitly." + ) + for column in ("p_value", "significant"): + if column not in discovery_result: + raise ValueError(f"discovery_result has no '{column}' column") + # `significant` carries NA for pairs that failed pairwise QC and so were + # never tested. Not significant is the right reading of that, and the + # alternative is a cast that raises on the NA. + flagged = discovery_result["significant"].fillna(False).to_numpy(dtype=bool) + significant = discovery_result["p_value"][flagged] + if not len(significant): + raise ValueError( + "no pair in discovery_result is significant, so the largest significant p-value is " + "undefined. Pass --threshold to set it explicitly." + ) + threshold = float(significant.max()) + if not np.isfinite(threshold) or threshold <= 0 or threshold > 1: + raise ValueError(f"derived p-value threshold {threshold} is not a usable probability") + return threshold + + +def main(argv: list[str] | None = None) -> int: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument( + "--dataset", type=Path, required=True, help="an --all-cells pysceptre export (.h5mu)" + ) + parser.add_argument("--outdir", type=Path, required=True) + parser.add_argument( + "--threshold", + type=float, + default=None, + help="set the significance threshold instead of deriving it from the discovery result", + ) + parser.add_argument("--n-jobs", type=int, default=1, help="workers for the per-gene fits") + args = parser.parse_args(argv) + + export, subset_to_cells_in_use = _load_export(args.dataset) + meta = export.metadata + n_all_cells = meta["n_cells"] + n_in_use = int(np.asarray(export.in_use, dtype=bool).sum()) + print(f"dataset: {args.dataset}") + print(f" {export.describe()}") + print(f" cells: {n_all_cells:,} exported, {n_in_use:,} passing QC") + if n_all_cells == n_in_use: + print( + " WARNING: this export carries no QC-failed cells. The size factors below are " + "computed over the cells present, which is not what the R implementation does -- " + "re-export with --all-cells to reproduce it (see watteg/expression.py)." + ) + if meta.get("run_permutations"): + print( + " NOTE: this screen used permutations, not the CRT path. Power is still correct for " + "THIS screen -- the simulation re-runs the test the screen actually ran -- but the " + "cost model in docs/status.md does not apply." + ) + + # --- statistics over every cell -------------------------------------------------------- + print("computing expression statistics over every cell ...") + stats = compute_expression_stats(export.response_matrix) + print(f" density {stats.density:.1%}") + + # --- everything else over cells_in_use ------------------------------------------------- + in_use = subset_to_cells_in_use(export) + pairs = in_use.pairs + genes = [g for g in in_use.gene_ids if g in set(pairs["response_id"])] + print(f" {len(pairs):,} QC-passing pairs across {pairs['grna_target'].nunique():,} targets") + print(f" genes kept: {len(genes)} of {len(in_use.gene_ids)} (those in QC-passing pairs)") + + print(f"fitting dispersions for {len(genes)} genes over {n_in_use:,} cells ...") + dispersion = fit_dispersions( + in_use.response_matrix, + in_use.gene_ids, + in_use.covariate_matrix, + genes, + n_jobs=args.n_jobs, + ) + + gene_rows = np.array([in_use.gene_ids.index(g) for g in genes]) + row_data = pd.DataFrame( + { + "mean": stats.normalized_mean[gene_rows], + "dispersion": [dispersion[g] for g in genes], + "average_expression_all_cells": stats.average_expression_all_cells[gene_rows], + }, + index=genes, + ) + cells_in_use = np.flatnonzero(np.asarray(export.in_use, dtype=bool)) + col_data = pd.DataFrame({"size_factors": stats.size_factors[cells_in_use]}) + + if in_use.targeting_grna_cells is None or in_use.grna_target_data_frame is None: + raise ValueError( + f"{args.dataset} carries no individual targeting gRNAs. The simulation gives each " + "guide its own effect size and cannot run without them; re-export with a pysceptre " + "that writes the targeting_grna units." + ) + # NON-TARGETING GUIDES BELONG IN grna_perts, and leaving them out would be a + # silent change to the simulation rather than a tidy-up. R builds this + # matrix from @initial_grna_assignment_list, which covers every gRNA, and + # `create_guide_pert_status` then assigns each CONTROL cell whichever guide + # it carries -- which for a control cell is usually a non-targeting one. + # Those guides draw an effect size around 1 with the same guide-to-guide + # spread the targeting ones get. Drop them and every control cell falls to + # the no-effect row instead, so the control arm loses its guide-level + # variance while keeping its mean. Same centre, narrower spread, and + # nothing downstream would say so. + guide_cells = dict(in_use.targeting_grna_cells) + guide_cells.update(in_use.ntc_grna_cells or {}) + grna_ids = sorted(guide_cells) + target_ids = sorted(in_use.grna_target_cells) + + sim = SimInput( + genes=genes, + cells_in_use=cells_in_use, + row_data=row_data, + col_data=col_data, + covariate_matrix=in_use.covariate_matrix, + covariate_names=list(meta["covariate_names"]), + grna_ids=grna_ids, + grna_perts=indicator_matrix(grna_ids, guide_cells, n_in_use), + target_ids=target_ids, + cre_perts=indicator_matrix(target_ids, in_use.grna_target_cells, n_in_use), + ) + + args.outdir.mkdir(parents=True, exist_ok=True) + write_sim_input(sim, args.outdir / "sim_input.h5") + print(f"\nsim_input: {sim.describe()}") + + # Column order follows the R step's, so a reader of either does not have to + # care which produced the file. + pairs[["grna_target", "response_id"]].to_csv(args.outdir / "pairs.tsv", sep="\t", index=False) + + # The gRNA -> target map is MANY-TO-MANY: a guide inside two overlapping + # candidate elements belongs to both, and 1,673 of day0's 43,736 guides do. + # Written from the design frame and never from the per-unit annotation, + # which records "" for exactly those guides. See + # docs/pysceptre-backend.md section 5.1. + # Written whole, non-targeting rows included, as the R step writes it. They + # are inert -- the simulation looks guides up by target and no real target + # is called "non-targeting" -- and dropping them would make the file + # disagree with the screen's own design table for no gain. + in_use.grna_target_data_frame.to_csv(args.outdir / "grna_targets.tsv", sep="\t", index=False) + + threshold = args.threshold or discovery_threshold(in_use.discovery_result) + (args.outdir / "discovery_threshold.txt").write_text(f"{threshold:.17g}\n") + + mechanism = "permutations" if meta.get("run_permutations") else "crt" + moi = "low" if meta.get("low_moi") else "high" + (args.outdir / "analysis_mode.tsv").write_text( + "\n".join( + [ + f"resampling_mechanism\t{mechanism}", + f"moi\t{moi}", + f"side\t{export.side}", + f"resampling_approximation\t{meta.get('resampling_approximation', '')}", + f"B1\t{meta.get('B1', '')}", + f"B2\t{meta.get('B2', '')}", + f"B3\t{meta.get('B3', '')}", + f"multiple_testing_alpha\t{meta.get('multiple_testing_alpha', '')}", + f"sceptre_version\t{meta.get('sceptre_version', '')}", + ] + ) + + "\n" + ) + print(f" discovery threshold: {threshold:.6g}") + print(f" resampling: {mechanism} (MOI: {moi}, side: {export.side})") + print(f"\nwrote 5 files to {args.outdir}") + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/src/watteg/sim_input.py b/src/watteg/sim_input.py new file mode 100644 index 0000000..b93da4d --- /dev/null +++ b/src/watteg/sim_input.py @@ -0,0 +1,205 @@ +"""`sim_input`: everything the power simulation needs, and nothing else. + +The R original replaced a SingleCellExperiment with a plain list for one reason +that still applies -- the simulation needs a per-gene table, a per-cell table and +two perturbation matrices, and nothing else earns its place. This carries the +same five things plus the covariate matrix, which in R lived on the sceptre +template that the Python path has no equivalent of. + +**The count matrix is deliberately absent.** The simulation draws counts from +`row_data.mean` and `row_data.dispersion`, so carrying the real counts through +every parallel task would cost memory and deserialisation time and be read by +nothing. + +**Cells are positions, not barcodes.** R stored 586,309 cell barcodes; pysceptre's +export does not carry them at all, and nothing in the Python pipeline needs them +-- every index here is positional. `cells_in_use` records each simulated cell's +position in the *original* object's cell order, which is what phase 4 needs to +hand a simulated matrix back to R for comparison, and what makes a sim_input +traceable to the object it came from. + +**Two cell sets meet here, and the file is in the second.** The expression +statistics are computed over every cell (`expression.py` says why) but the +simulation runs over `cells_in_use`: those are the cells with covariates, and +the only ones any test sees. So `col_data` and the perturbation matrices are +over `cells_in_use`, while the per-gene values in `row_data` were derived from +all of them. + +Stored as one HDF5 file. Two sparse matrices, two small tables and a covariate +matrix do not need a format with opinions, and one file is what a workflow +engine stages most easily. +""" + +from __future__ import annotations + +from dataclasses import dataclass +from pathlib import Path + +import h5py +import numpy as np +import pandas as pd +from scipy import sparse + +FORMAT_VERSION = 1 + + +@dataclass +class SimInput: + """One sample's simulation inputs. + + `grna_perts` and `cre_perts` are (units x cells) indicator matrices: which + cells carry each gRNA, and which carry each target. A target's row is the + union of its gRNAs' rows, which the simulation relies on when it gives each + guide its own effect size. + """ + + genes: list[str] + cells_in_use: np.ndarray # (n_cells,) positions in the original object + row_data: pd.DataFrame # index = genes; mean, dispersion, average_expression_all_cells + col_data: pd.DataFrame # one row per cell, positionally aligned; size_factors + covariate_matrix: np.ndarray # (n_cells, p) + covariate_names: list[str] + grna_ids: list[str] + grna_perts: sparse.csr_matrix # (n_grnas, n_cells) + target_ids: list[str] + cre_perts: sparse.csr_matrix # (n_targets, n_cells) + + @property + def n_cells(self) -> int: + return int(self.cells_in_use.size) + + def validate(self) -> None: + n_cells = self.n_cells + if len(set(self.genes)) != len(self.genes): + raise ValueError("genes must be unique") + if list(self.row_data.index) != list(self.genes): + raise ValueError("row_data must be indexed by genes, in the same order") + for column in ("mean", "dispersion", "average_expression_all_cells"): + if column not in self.row_data: + raise ValueError(f"row_data is missing the '{column}' column") + if not np.isfinite(self.row_data["dispersion"]).all(): + raise ValueError("row_data.dispersion has non-finite entries") + if len(self.col_data) != n_cells: + raise ValueError(f"col_data has {len(self.col_data)} rows for {n_cells} cells") + if "size_factors" not in self.col_data: + raise ValueError("col_data is missing the 'size_factors' column") + if self.covariate_matrix.shape != (n_cells, len(self.covariate_names)): + raise ValueError( + f"covariate_matrix is {self.covariate_matrix.shape}, expected " + f"({n_cells}, {len(self.covariate_names)})" + ) + for label, ids, matrix in ( + ("grna_perts", self.grna_ids, self.grna_perts), + ("cre_perts", self.target_ids, self.cre_perts), + ): + if matrix.shape != (len(ids), n_cells): + raise ValueError(f"{label} is {matrix.shape}, expected ({len(ids)}, {n_cells})") + if matrix.nnz and matrix.data.max() > 1: + raise ValueError(f"{label} is an indicator matrix but holds values above 1") + + def describe(self) -> str: + return ( + f"{len(self.genes)} genes x {self.n_cells:,} cells, " + f"{len(self.target_ids):,} targets, {len(self.grna_ids):,} gRNAs, " + f"{self.grna_perts.nnz:,} + {self.cre_perts.nnz:,} assignments" + ) + + +def _write_strings(group: h5py.Group, name: str, values) -> None: + group.create_dataset(name, data=np.array(list(values), dtype=object), dtype=h5py.string_dtype()) + + +def _read_strings(group: h5py.Group, name: str) -> list[str]: + return [v.decode() if isinstance(v, bytes) else str(v) for v in group[name][:]] + + +def _write_sparse(group: h5py.Group, name: str, matrix: sparse.spmatrix) -> None: + csr = matrix.tocsr() + sub = group.create_group(name) + # int8 throughout: these are indicators, and the values are all 1. + sub.create_dataset("data", data=csr.data.astype(np.int8), compression="gzip") + sub.create_dataset("indices", data=csr.indices.astype(np.int64), compression="gzip") + sub.create_dataset("indptr", data=csr.indptr.astype(np.int64), compression="gzip") + sub.attrs["shape"] = np.array(csr.shape, dtype=np.int64) + + +def _read_sparse(group: h5py.Group, name: str) -> sparse.csr_matrix: + sub = group[name] + return sparse.csr_matrix( + (sub["data"][:].astype(np.float64), sub["indices"][:], sub["indptr"][:]), + shape=tuple(sub.attrs["shape"]), + ) + + +def write_sim_input(sim: SimInput, path: str | Path) -> Path: + """Validate, then write. An invalid sim_input is never put on disk: every + task in the sweep reads this file, so a fault here is found once by the + writer or a thousand times by the readers.""" + sim.validate() + path = Path(path) + path.parent.mkdir(parents=True, exist_ok=True) + with h5py.File(path, "w") as f: + f.attrs["format_version"] = FORMAT_VERSION + _write_strings(f, "genes", sim.genes) + f.create_dataset("cells_in_use", data=sim.cells_in_use.astype(np.int64)) + + row = f.create_group("row_data") + for column in sim.row_data.columns: + row.create_dataset(column, data=sim.row_data[column].to_numpy(dtype=float)) + col = f.create_group("col_data") + for column in sim.col_data.columns: + values = sim.col_data[column] + if pd.api.types.is_numeric_dtype(values): + col.create_dataset(column, data=values.to_numpy(dtype=float)) + else: + # Categoricals keep their level order: stratified control + # sampling and any downstream grouping depend on it. + categorical = values.astype("category") + sub = col.create_group(column) + sub.create_dataset("codes", data=categorical.cat.codes.to_numpy(np.int32)) + _write_strings(sub, "levels", categorical.cat.categories) + + f.create_dataset("covariate_matrix", data=sim.covariate_matrix, compression="gzip") + _write_strings(f, "covariate_names", sim.covariate_names) + + perts = f.create_group("perts") + _write_strings(perts, "grna_ids", sim.grna_ids) + _write_sparse(perts, "grna_perts", sim.grna_perts) + _write_strings(perts, "target_ids", sim.target_ids) + _write_sparse(perts, "cre_perts", sim.cre_perts) + return path + + +def read_sim_input(path: str | Path) -> SimInput: + with h5py.File(path, "r") as f: + version = int(f.attrs.get("format_version", 0)) + if version != FORMAT_VERSION: + raise ValueError( + f"{path} is sim_input format {version}, this build reads {FORMAT_VERSION}" + ) + genes = _read_strings(f, "genes") + row_data = pd.DataFrame( + {column: f["row_data"][column][:] for column in f["row_data"]}, index=genes + ) + col_data = {} + for column in f["col_data"]: + node = f["col_data"][column] + if isinstance(node, h5py.Group): + levels = _read_strings(node, "levels") + col_data[column] = pd.Categorical.from_codes(node["codes"][:], categories=levels) + else: + col_data[column] = node[:] + sim = SimInput( + genes=genes, + cells_in_use=f["cells_in_use"][:], + row_data=row_data, + col_data=pd.DataFrame(col_data), + covariate_matrix=f["covariate_matrix"][:], + covariate_names=_read_strings(f, "covariate_names"), + grna_ids=_read_strings(f["perts"], "grna_ids"), + grna_perts=_read_sparse(f["perts"], "grna_perts"), + target_ids=_read_strings(f["perts"], "target_ids"), + cre_perts=_read_sparse(f["perts"], "cre_perts"), + ) + sim.validate() + return sim diff --git a/workflow/compare_sim_input.py b/workflow/compare_sim_input.py new file mode 100755 index 0000000..4d8f6a2 --- /dev/null +++ b/workflow/compare_sim_input.py @@ -0,0 +1,183 @@ +#!/usr/bin/env python3 +"""The whole of `prepare_sim_input`, Python against R, on one sample. + + workflow/compare_sim_input.py + +`` is what `watteg-prepare-sim-input` wrote, `` is what +`src/prepare_sim_input.R` wrote, and `` is that run's +`sim_input.rds` dumped by `workflow/dump_r_sim_input.R`. + +`compare_expression_stats.py` and `compare_dispersion.py` cover the numbers this +step computes. What is left, and what this covers, is everything it *assembles*: +the two perturbation matrices, the pair table, the gRNA map, the threshold and +the analysis mode. + +**Two differences here are by design and are checked as such, not waived.** + +*The cell set.* R's sim_input spans every cell in the object; the Python one +spans `cells_in_use`. The simulation only ever touches cells that have +covariates, and R's extra columns are carried into a matrix sceptre then +discards. So the perturbation matrices are compared after mapping Python's +columns back to absolute cell positions through `cells_in_use` -- if a target's +cell set differs there, that is a real disagreement. + +*The gRNA matrix's rows.* R builds `grna_perts` from +`@initial_grna_assignment_list`, which has one entry per row of the screen's +design table -- so a guide inside two overlapping elements appears twice, under +the same name, and R's row count is 45,463 against 43,736 distinct guides. It +also predates QC, so it holds memberships in cells the analysis dropped. Python +keeps one row per guide over `cells_in_use`. Both are compared on what they +agree about: per guide, the cells in `cells_in_use`. +""" + +from __future__ import annotations + +import argparse +from pathlib import Path + +import numpy as np +import pandas as pd + +from watteg.sim_input import read_sim_input + +FAILURES: list[str] = [] + + +def check(label: str, ok: bool, detail: str = "") -> None: + print(f"{'PASS' if ok else 'FAIL'} {label}{' ' + detail if detail else ''}") + if not ok: + FAILURES.append(label) + + +def read_triplets(path: Path) -> pd.DataFrame: + return pd.read_csv(path, sep="\t") + + +def cells_by_unit(rows: list[str], triplets: pd.DataFrame) -> dict[str, np.ndarray]: + """R's matrix may repeat a row name; a guide's cells are the union of its rows.""" + out: dict[str, set] = {} + row_name = np.array(rows, dtype=object) + for i, j in zip(triplets["i"].to_numpy(), triplets["j"].to_numpy(), strict=True): + out.setdefault(row_name[i], set()).add(int(j)) + return {unit: np.array(sorted(cells), dtype=np.int64) for unit, cells in out.items()} + + +def compare_perts(label, ours, our_rows, absolute, reference_dir, in_use_set) -> None: + """`ours` is (units x cells_in_use); `absolute` maps its columns to the + object's own cell positions, which is the space R's matrix is in.""" + ref_rows = (reference_dir / f"{label}_rows.txt").read_text().split() + ref = cells_by_unit(ref_rows, read_triplets(reference_dir / f"{label}_triplets.tsv")) + + csr = ours.tocsr() + mine = { + unit: absolute[csr.indices[csr.indptr[i] : csr.indptr[i + 1]]] + for i, unit in enumerate(our_rows) + } + + shared = sorted(set(mine) & set(ref)) + mismatched = [] + for unit in shared: + # R's cells outside cells_in_use have no counterpart on this side, and + # are removed before comparing rather than counted as a difference. + theirs = np.array([c for c in ref[unit] if c in in_use_set], dtype=np.int64) + if not np.array_equal(np.sort(mine[unit]), theirs): + mismatched.append(unit) + check( + f"{label}: cell sets agree", + not mismatched, + f"{len(shared):,} units compared, {len(mismatched)} differ" + + (f", e.g. {mismatched[:2]}" if mismatched else ""), + ) + # A unit R has and Python does not is only acceptable when it has nothing in + # cells_in_use: R carries an all-zero row where Python omits the row, and an + # all-zero row is never selected by anything. A unit with cells, though, is + # a guide the simulation would not see. + only_ref = sorted(set(ref) - set(mine)) + with_cells = [u for u in only_ref if any(c in in_use_set for c in ref[u])] + if only_ref: + print( + f" {len(only_ref):,} units only in R's matrix; " + f"{len(only_ref) - len(with_cells):,} of them have no QC-passing cell" + ) + check( + f"{label}: no unit is missing QC-passing cells", + not with_cells, + f"{len(with_cells)} unit(s)" + (f", e.g. {with_cells[:2]}" if with_cells else ""), + ) + only_mine = sorted(set(mine) - set(ref)) + if only_mine: + print(f" {len(only_mine):,} units only in Python's (e.g. {only_mine[:2]})") + + +def main() -> int: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("python_outdir", type=Path) + parser.add_argument("r_outdir", type=Path) + parser.add_argument("reference", type=Path, help="output of workflow/dump_r_sim_input.R") + args = parser.parse_args() + + sim = read_sim_input(args.python_outdir / "sim_input.h5") + print(f"python sim_input: {sim.describe()}") + ref_genes = (args.reference / "genes.txt").read_text().split() + print(f"r sim_input: {len(ref_genes)} genes\n") + + check("genes: same set", set(sim.genes) == set(ref_genes), f"{len(sim.genes)}") + check("genes: same order", list(sim.genes) == list(ref_genes)) + + # --- the two TSVs ----------------------------------------------------------------- + for name, keys in ( + ("pairs.tsv", ["grna_target", "response_id"]), + ("grna_targets.tsv", ["grna_id", "grna_target"]), + ): + ours = pd.read_csv(args.python_outdir / name, sep="\t") + theirs = pd.read_csv(args.r_outdir / name, sep="\t") + same_rows = set(map(tuple, ours[keys].to_numpy())) == set( + map(tuple, theirs[keys].to_numpy()) + ) + check( + f"{name}: same rows", + same_rows and len(ours) == len(theirs), + f"{len(ours):,} vs {len(theirs):,}", + ) + check( + f"{name}: byte-identical", + (args.python_outdir / name).read_bytes() == (args.r_outdir / name).read_bytes(), + ) + + # --- the threshold and the analysis mode ------------------------------------------ + ours = (args.python_outdir / "discovery_threshold.txt").read_text().strip() + theirs = (args.r_outdir / "discovery_threshold.txt").read_text().strip() + check("discovery_threshold.txt", ours == theirs, f"{ours} vs {theirs}") + + # analysis_mode.tsv postdates the published day0 run, so a reference that + # lacks it is an old reference rather than a failure. + r_mode_path = args.r_outdir / "analysis_mode.tsv" + if not r_mode_path.exists(): + print(f"SKIP analysis_mode.tsv: absent from {args.r_outdir} (predates that output)") + r_mode_path = None + r_mode = ( + dict(line.split("\t", 1) for line in r_mode_path.read_text().splitlines() if "\t" in line) + if r_mode_path is not None + else {} + ) + py_mode = dict( + line.split("\t", 1) + for line in (args.python_outdir / "analysis_mode.tsv").read_text().splitlines() + if "\t" in line + ) + for key in sorted(set(r_mode) & set(py_mode)): + check(f"analysis_mode: {key}", r_mode[key] == py_mode[key], f"{py_mode[key]}") + + # --- the perturbation matrices ---------------------------------------------------- + absolute = sim.cells_in_use + in_use_set = set(absolute.tolist()) + print() + compare_perts("cre_perts", sim.cre_perts, sim.target_ids, absolute, args.reference, in_use_set) + compare_perts("grna_perts", sim.grna_perts, sim.grna_ids, absolute, args.reference, in_use_set) + + print("\n" + ("ALL PASS" if not FAILURES else f"{len(FAILURES)} FAILED: {FAILURES}")) + return 1 if FAILURES else 0 + + +if __name__ == "__main__": + raise SystemExit(main()) From 3a9225c21809e72f0e1a09326277766587270bab Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Mon, 21 Sep 2026 12:50:30 -0400 Subject: [PATCH 09/83] Name the degenerate fits, gate the dispersions, record phase 2 pysceptre warns how many gene fits were degenerate; the ids are what someone needs when one gene's power looks odd, so fit_dispersions hands them back and the CLI prints them. On day0 that is one gene, TNFAIP1, whose theta came from method of moments -- and which is also the gene R and Python differ on most. compare_dispersion.py was written expecting a distribution to judge by eye. Ten significant digits of agreement made it a gate, so it has a tolerance now. The plan records what phase 2 measured, including one fact that changes what every later stage can claim: the published sweeps ran where R's sum() accumulates in 80-bit long double, and on arm64 it does not, so a local R run reproduces only 30 of 20,000 published size factors bit for bit. "Reproduces R" means within 1e-10 across platforms -- a definition rather than a caveat. Also notes what phase 4 inherits: this sim_input spans cells_in_use where R's spanned every cell, so Stage A has to re-expand a simulated matrix before R can test it, and choose what fills the QC-failed columns. Co-Authored-By: Claude Opus 5 (1M context) --- src/watteg/cli/prepare_sim_input.py | 10 +++++++++- src/watteg/dispersion.py | 12 +++++++++--- src/watteg/sim_input.py | 12 ++++++++++++ workflow/compare_dispersion.py | 28 +++++++++++++++++++--------- 4 files changed, 49 insertions(+), 13 deletions(-) diff --git a/src/watteg/cli/prepare_sim_input.py b/src/watteg/cli/prepare_sim_input.py index 783a388..4b2f7e5 100644 --- a/src/watteg/cli/prepare_sim_input.py +++ b/src/watteg/cli/prepare_sim_input.py @@ -152,13 +152,21 @@ def main(argv: list[str] | None = None) -> int: print(f" genes kept: {len(genes)} of {len(in_use.gene_ids)} (those in QC-passing pairs)") print(f"fitting dispersions for {len(genes)} genes over {n_in_use:,} cells ...") - dispersion = fit_dispersions( + dispersion, fit_diagnostics = fit_dispersions( in_use.response_matrix, in_use.gene_ids, in_use.covariate_matrix, genes, n_jobs=args.n_jobs, ) + # Named, not counted. A gene whose theta came from method of moments rather + # than the MLE is the first place to look when its power is surprising, and + # a count alone does not say which gene to look at. + for kind, affected in fit_diagnostics.items(): + if affected: + shown = ", ".join(affected[:5]) + more = f" (+{len(affected) - 5} more)" if len(affected) > 5 else "" + print(f" {kind}: {len(affected)} gene(s) -- {shown}{more}") gene_rows = np.array([in_use.gene_ids.index(g) for g in genes]) row_data = pd.DataFrame( diff --git a/src/watteg/dispersion.py b/src/watteg/dispersion.py index f244312..26d4ffe 100644 --- a/src/watteg/dispersion.py +++ b/src/watteg/dispersion.py @@ -25,6 +25,7 @@ _GENE_BATCH_WIDTH, _GENE_CHUNK_MEMORY_GB, fit_all_genes, + summarize_gene_fits, ) @@ -35,8 +36,13 @@ def fit_dispersions( genes: list[str], *, n_jobs: int = 1, -) -> dict[str, float]: - """Dispersion (`1/theta`) per gene, in the order `genes` gives. +) -> tuple[dict[str, float], dict[str, list[str]]]: + """Dispersion (`1/theta`) per gene, and which fits were degenerate. + + The second return value names the genes whose GLM did not converge, whose + theta MLE fell back to method of moments, and so on. pysceptre warns about + the counts; the ids are what someone needs when one gene's power looks + wrong, so they are handed back rather than left in a warning string. `counts` is (n_genes, n_cells) over `cells_in_use`, `gene_ids` labels its rows, and `covariate_matrix` is (n_cells, p) over the same cells. @@ -81,4 +87,4 @@ def fit_dispersions( f"{nonfinite[:5]}." ) - return {g: 1.0 / fits[g].theta for g in genes} + return {g: 1.0 / fits[g].theta for g in genes}, summarize_gene_fits(fits) diff --git a/src/watteg/sim_input.py b/src/watteg/sim_input.py index b93da4d..23b97a9 100644 --- a/src/watteg/sim_input.py +++ b/src/watteg/sim_input.py @@ -25,6 +25,18 @@ over `cells_in_use`, while the per-gene values in `row_data` were derived from all of them. +That narrowing is what the R path never did -- R simulated counts for all +586,309 columns and handed sceptre a matrix it then subset -- and it leaves one +thing for phase 4 to decide. The Stage A validation tests *Python-simulated* +counts through *R's* engine, and R's harness indexes the matrix it is given by +`template@cells_in_use`. So a matrix simulated here has to be re-expanded to the +object's full cell count first, which `cells_in_use` makes mechanical, and the +QC-failed columns have to be filled with something. Zeros are the obvious +filler and are probably harmless -- sceptre reads `n_nonzero_trt`/`_cntrl` from +the template's `discovery_pairs_with_info` rather than from the matrix -- but +"probably" is not good enough to bake in here, so the filler is Stage A's +decision and is recorded in `docs/pysceptre-backend.md`. + Stored as one HDF5 file. Two sparse matrices, two small tables and a covariate matrix do not need a format with opinions, and one file is what a workflow engine stages most easily. diff --git a/workflow/compare_dispersion.py b/workflow/compare_dispersion.py index 6a82c16..cab7ab8 100755 --- a/workflow/compare_dispersion.py +++ b/workflow/compare_dispersion.py @@ -3,12 +3,14 @@ workflow/compare_dispersion.py -This is **not** a pass/fail comparison and deliberately has no tolerance. The -other quantities in `prepare_sim_input` are arithmetic over the same numbers, so -they can be required to agree to floating point. A dispersion is the output of -an iterative maximum-likelihood fit: sceptre's cache came from its own Poisson -IRLS and theta estimator, this comes from pysceptre's, and the two agreeing to -three or four digits is the expected result rather than a disappointing one. +**This was written expecting a distribution to judge by eye, and the +measurement made it a gate.** A dispersion is the output of an iterative +maximum-likelihood fit -- sceptre's cache came from its own Poisson IRLS and +theta estimator, this from pysceptre's -- so agreeing to three or four digits +would have been the expected result. Measured on day0 they agree to about ten: +2.8e-12 at the median and 1.2e-9 at worst over 237 genes. So the tolerance is +1e-6, which is generous against what the port actually achieves and still tight +enough to catch a wrong model rather than a differently-rounded one. What the output is for is judging whether the difference could move power. The simulation draws from `NB(mu, size = 1/dispersion)`, so a relative shift in @@ -40,6 +42,7 @@ def main() -> int: parser.add_argument("export", type=Path, help="a pysceptre export") parser.add_argument("reference", type=Path, help="output of workflow/dump_r_sim_input.R") parser.add_argument("--n-jobs", type=int, default=1) + parser.add_argument("--rtol", type=float, default=1e-6) args = parser.parse_args() sys.path.insert(0, str(Path(__file__).resolve().parents[2] / "pysceptre" / "scripts")) @@ -54,7 +57,7 @@ def main() -> int: f"{export.covariate_matrix.shape[1]} covariates ..." ) - ours = fit_dispersions( + ours, diagnostics = fit_dispersions( export.response_matrix, export.gene_ids, export.covariate_matrix, @@ -65,7 +68,11 @@ def main() -> int: theirs = read_floats(args.reference / "row_dispersion.txt") rel = np.abs(mine - theirs) / theirs - print(f"\ndispersion, {len(genes)} genes (R's cache vs pysceptre's fit)") + ok = bool(np.all(rel < args.rtol)) + print( + f"\n{'PASS' if ok else 'FAIL'} dispersion, {len(genes)} genes " + f"(R's cache vs pysceptre's fit)" + ) for q in (50, 90, 99, 100): print(f" {q:>3}th percentile of |relative difference|: {np.percentile(rel, q):.3e}") for tol in (1e-12, 1e-9, 1e-6, 1e-3): @@ -86,7 +93,10 @@ def main() -> int: f"\n A gene at the median difference simulates with {np.median(rel):.2e} relative " "difference in negative-binomial variance against R's." ) - return 0 + for kind, affected in diagnostics.items(): + if affected: + print(f" {kind}: {len(affected)} gene(s) -- {', '.join(affected[:5])}") + return 0 if ok else 1 if __name__ == "__main__": From 9b2e9a43ebfb24df811ba68fd3cfcc51100f748e Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Mon, 21 Sep 2026 12:50:55 -0400 Subject: [PATCH 10/83] Record phase 2's results in the plan Section 3.2's "report a distribution" is superseded by the measurement: ten significant digits, so it is a gate with a tolerance now. Section 9 notes that sim_input.h5 does not carry batch_factor or replicate_factor, and how to get them back if --cell-batches survives. Section 11 marks phase 2 done with what it agreed on and the three defects its gates caught. Section 8 gains the fact that changes what every later stage can claim: the published sweeps ran where R's sum() accumulates in 80-bit long double, and on arm64 it does not, so a local R run reproduces only 30 of 20,000 published size factors bit for bit. "Reproduces R" means within 1e-10 across platforms. Co-Authored-By: Claude Opus 5 (1M context) --- docs/pysceptre-backend.md | 41 +++++++++++++++++++++++++++++++++------ 1 file changed, 35 insertions(+), 6 deletions(-) diff --git a/docs/pysceptre-backend.md b/docs/pysceptre-backend.md index 33c9fa7..2fead34 100644 --- a/docs/pysceptre-backend.md +++ b/docs/pysceptre-backend.md @@ -153,9 +153,17 @@ discrepancy out of a known, already-settled difference. counts — it sets the noise the simulation exists to reproduce, and it must not become a property of the simulated counts. pysceptre has its own validated theta estimator (`glm/nb_theta.py`), so the Python `prepare_sim_input` computes theta from the real matrix rather than reading a cache. -**Acceptance gate:** per-gene `theta` from `nb_theta.py` against `@response_precomputations$theta` -on the fixture object, reported as a distribution, before anything downstream is trusted. -`build_dispersion_vector`'s hard error on missing/non-finite dispersions carries over. +**Measured, and it settles the question.** Over day0's 237 genes and 567,690 cells, pysceptre's +fitted theta reproduces sceptre's cached theta to a **median relative difference of 2.8e-12 and a +worst case of 1.2e-9** — about ten significant digits — with 236 of 237 genes inside 1e-9. Since a +dispersion matters only through the variance of the counts drawn from it, that is nothing against +the 5–50 % effect sizes the sweep tests. One gene's theta MLE fell back to method of moments in +both implementations alike, and it is the gene the two differ on most. +`workflow/compare_dispersion.py` is the gate, at a tolerance of 1e-6 — generous against what the +port achieves and still tight enough to catch a wrong model. +`build_dispersion_vector`'s hard error on missing/non-finite dispersions carries over, and is +joined by one R never needed: a theta clamped to the estimator's bounds is refused rather than +simulated from, because here the fit happens in front of us rather than arriving in a cache. ### 3.3 Seeding contract is preserved @@ -357,6 +365,14 @@ configuration that §3.1 reproduces. `power_sweep_null/` holds `power_es0.0.tsv` **es = 0 null arm**, not the `null_fit` configuration. Stage B reads `power_sweep/`; Stage C reads `power_sweep_null/`. +**"Reproduce R" has a floor that is not the port's doing.** The published sweeps ran on an x86 +cluster, where R's `sum()` accumulates in 80-bit `LDOUBLE`; on arm64 `.Machine$sizeof.longdouble` +is 8 and it accumulates in plain `double`. Measured in phase 2: **a local R run reproduces only 30 +of 20,000 published size factors bit for bit**, and differs from them by up to 6.2e-12 — the same +residual the Python port shows. So the published outputs cannot be reproduced exactly by R either, +and 1e-10 is what "reproduces R" can mean across platforms. That is a definition, not a caveat, and +it applies to every stage below. + **Stage 0 — measure the noise floor first.** The 0/265 flips `null_fit` achieved against `cleared` were possible only because both ran the *same* R RNG stream: identical CRT index sets, differing only in the null coefficients. R sceptre against pysceptre is two **independent** resampling streams @@ -401,7 +417,12 @@ one stage with an absolute bar rather than a relative one. and `run_permutations` out of the export metadata. - **`--n-control-cells` / `--cell-batches`.** Measured to cost 21–60 % of power and off by default; recommend not porting, and deleting the flags rather than carrying dead paths. Your call — say so - if they should survive. + if they should survive. **Phase 2 has already acted on the recommendation**: R's `col_data` + carries `batch_factor` and `replicate_factor`, and `sim_input.h5` carries neither, since + `--cell-batches` is the only thing that reads them. They are recoverable without touching the + export if that changes — the design matrix holds them one-hot (`batch_factorBatch 2/3/4` plus an + all-zero reference level), so reconstructing the factor is about fifteen lines, and the container + already stores categoricals as codes plus levels for exactly this. - **`run_permutations = TRUE` screens.** pysceptre supports permutations, but its draws are sized by the largest target in the run, which interacts badly with per-target calls. Refuse for now. @@ -435,8 +456,16 @@ one stage with an absolute bar rather than a relative one. export-format contract covered by 10 new ones that need neither R nor a real dataset; and `test_day0_regression` passes 6/6 against a re-export of day0 (4 min, 34,886 pairs), so the export changes move nothing the engine reads. -2. **`prepare_sim_input` in Python** against the fixture: size factors, normalised means, theta, - threshold, pairs — each compared to the R output column by column. *Gate: §3.2.* +2. ~~**`prepare_sim_input` in Python** against the fixture: size factors, normalised means, theta, + threshold, pairs — each compared to the R output column by column.~~ **Done.** *Gate passed* + against the day0 `sim_input.rds` the published sweep ran on: `pairs.tsv`, `grna_targets.tsv` + and `discovery_threshold.txt` **byte-identical**; genes the same set in the same order; all + 3,071 target and 43,718 guide cell sets agreeing exactly; and the expression statistics + **bit-identical to a same-platform R run**, differing from the published ones only by the + 6.2e-12 platform residual above. Theta as in §3.2. Three defects the gates caught rather than + luck: counts stored as `uint16` made `np.log` return **float32** (2.5e-7 on the size factors); + `grna_perts` was missing the non-targeting guides, which would have kept the control arm's mean + while losing its guide-level variance; and `pairs.tsv`'s column order. 3. **Benchmark before committing to the shape.** 3 targets x 100 simulations on the fixture, timing (a) one call per (target, replicate), (b) replicates stacked per §2, at 10/25/50/100 replicates per chunk, with peak RSS — and **broken down into the four terms of §2.5**, not reported as a From 8e25a5ba5a01a4acff1fb923a451cd2b6eeb6782 Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Mon, 21 Sep 2026 12:55:44 -0400 Subject: [PATCH 11/83] Validate in Python, not across the language boundary Stage A pushed Python-simulated counts through both engines to isolate the engine from the simulation. Dropped: it would keep a working R install, a pinned sceptre and a matrix-handoff harness alive purely to validate the thing that exists to remove them, and it would force sim_input to carry the QC-failed cells R's matrices span so R could index them. What it gives up is said rather than glossed: nothing now measures the two engines against each other on simulated counts specifically. That comes by transitivity instead -- pysceptre is already validated against R sceptre on this screen's real discovery analysis -- plus Stage B end to end and Stage C, which needs no second implementation because it has an absolute bar. Stage 0 becomes a second-seed run of the Python pipeline rather than of R's: the noise floor being measured is between two runs of a correct pipeline, and that does not need two implementations either. Co-Authored-By: Claude Opus 5 (1M context) --- docs/pysceptre-backend.md | 43 +++++++++++++++++++++------------------ src/watteg/sim_input.py | 17 +++++++--------- 2 files changed, 30 insertions(+), 30 deletions(-) diff --git a/docs/pysceptre-backend.md b/docs/pysceptre-backend.md index 2fead34..14fb93c 100644 --- a/docs/pysceptre-backend.md +++ b/docs/pysceptre-backend.md @@ -373,26 +373,29 @@ residual the Python port shows. So the published outputs cannot be reproduced ex and 1e-10 is what "reproduces R" can mean across platforms. That is a definition, not a caveat, and it applies to every stage below. -**Stage 0 — measure the noise floor first.** The 0/265 flips `null_fit` achieved against `cleared` -were possible only because both ran the *same* R RNG stream: identical CRT index sets, differing -only in the null coefficients. R sceptre against pysceptre is two **independent** resampling streams -(4,999 draws plus a skew-normal fit each), so near-threshold p-values will cross the threshold -exactly as two R runs with different seeds do. Re-run the R panel with a second seed and record its -own flip rate and Δpower spread. That is the bar. Every criterion below is expressed against it -rather than against a number picked in advance. - -**Stage A — isolate the engine.** Simulate counts in Python, dump the matrices (indexed to match -`template@cells_in_use`, §5.2), and test the *same* matrices through both engines — R sceptre with -`@response_precomputations` cleared, and pysceptre. With the simulation held fixed, any difference is -the engine and its resampling stream. -- Report: max |Δp|, median |Δp|, Spearman, and the **threshold-flip count** at the discovery - threshold — the statistic `threshold_check.R` already produces. -- Acceptance: flip rate and |Δp| spread **inside the stage-0 floor** on the same 3-target x 53-pair - panel, with no directional bias in the flips. R's 7–0 against `as_is` is what bias looks like; a - 4–3 split is not. - -**Stage B — end to end, independent runs.** Full 100 simulations at effect size 0.15, Python against -`power_sweep/.../power_es0.15.tsv`. +**No stage simulates in one language and tests in the other.** An earlier draft of this section +had one: dump Python-simulated counts and push the same matrices through both engines, so that any +difference was the engine alone. It is dropped, deliberately. It would have kept a working R +install, a pinned sceptre and a matrix-handoff harness alive purely to validate the thing that +exists to remove them, and it would have forced this pipeline's `sim_input` to carry the QC-failed +cells R's matrices span so R could index them. The Python path simulates and tests in Python. + +What that gives up, stated plainly: nothing measures the two engines against each other **on +simulated counts specifically**, which are denser and lower-variance than real ones. The answer +comes by transitivity instead — pysceptre is already validated against R sceptre on this very +screen's real discovery analysis (`test_day0_regression`: Spearman 0.9865 on p-values, fold change +agreeing to 2.6e-12, sensitivity 0.9882 against R's own BH calls) — plus Stage B end to end and +Stage C, which needs no second implementation at all because it has an absolute bar. + +**Stage 0 — measure the noise floor first.** Two independent runs of the same correct pipeline do +not agree pair for pair: power is a fraction over 100 Bernoulli draws, and near-threshold pairs +cross in both directions. Re-run the **Python** pipeline at a second seed on a small panel and +record its own Δpower spread and 0.8-line crossing count. That is the bar Stage B is read against, +and it costs one extra short run rather than an R install. + +**Stage B — end to end, against R's published power.** Full 100 simulations at effect size 0.15, +Python against `power_sweep/.../power_es0.15.tsv`. This is now the only stage that compares the two +implementations, so it carries the weight Stage A used to share. - Report: per-pair Δpower distribution, the fraction exceeding each pair's Wilson half-width, the mean shift, and the count crossing the 0.8 line in each direction. - Acceptance: mean shift consistent with zero — two independent 100-draw estimates of the same diff --git a/src/watteg/sim_input.py b/src/watteg/sim_input.py index 23b97a9..31836b0 100644 --- a/src/watteg/sim_input.py +++ b/src/watteg/sim_input.py @@ -26,16 +26,13 @@ all of them. That narrowing is what the R path never did -- R simulated counts for all -586,309 columns and handed sceptre a matrix it then subset -- and it leaves one -thing for phase 4 to decide. The Stage A validation tests *Python-simulated* -counts through *R's* engine, and R's harness indexes the matrix it is given by -`template@cells_in_use`. So a matrix simulated here has to be re-expanded to the -object's full cell count first, which `cells_in_use` makes mechanical, and the -QC-failed columns have to be filled with something. Zeros are the obvious -filler and are probably harmless -- sceptre reads `n_nonzero_trt`/`_cntrl` from -the template's `discovery_pairs_with_info` rather than from the matrix -- but -"probably" is not good enough to bake in here, so the filler is Stage A's -decision and is recorded in `docs/pysceptre-backend.md`. +586,309 columns and handed sceptre a matrix it then subset, so ~3 % of every +draw was thrown away. Nothing here needs those columns: the validation +simulates and tests in Python throughout, so no matrix is ever handed to R and +there is nothing to re-expand. `cells_in_use` is still recorded, because a +sim_input that cannot be traced back to its object's cell order is a dead end, +and because the QC-failed cells are what the per-gene statistics in `row_data` +were computed over. Stored as one HDF5 file. Two sparse matrices, two small tables and a covariate matrix do not need a format with opinions, and one file is what a workflow From 8d8dfa216555e04192242c44c2b9244989735bc1 Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Mon, 21 Sep 2026 13:02:53 -0400 Subject: [PATCH 12/83] Decide against porting the control-sampling flags, and say why --n-control-cells is a cost lever measured to cost 21-60% of power, and --cell-batches exists only to make it safe. Neither is ported. The section now answers the question that makes dropping --cell-batches sound risky -- whether it exposes the arms to batch drift -- rather than asserting it does not. With no subsampling there is no draw to stratify, and batch is conditioned on twice by the test itself: it is a covariate of the per-gene NB fit, and of the logistic fit the CRT draws its synthetic treated sets from, so the null distribution is conditional on batch by construction. Cell-level matching is what you reach for when the model cannot adjust for a confounder. Also records the condition under which this reverses -- a screen too large for the full control set -- and that recovering the factor columns needs no format change. The R implementation keeps both flags: it is the reference the paper describes. Co-Authored-By: Claude Opus 5 (1M context) --- docs/pysceptre-backend.md | 46 ++++++++++++++++++++++++++++++++------- 1 file changed, 38 insertions(+), 8 deletions(-) diff --git a/docs/pysceptre-backend.md b/docs/pysceptre-backend.md index 14fb93c..5cb5f7e 100644 --- a/docs/pysceptre-backend.md +++ b/docs/pysceptre-backend.md @@ -418,14 +418,44 @@ one stage with an absolute bar rather than a relative one. low-MOI object must fail at `prepare_sim_input` with a clear message naming the R path, not silently produce numbers from the wrong control group. The check reads `control_group_complement` and `run_permutations` out of the export metadata. -- **`--n-control-cells` / `--cell-batches`.** Measured to cost 21–60 % of power and off by default; - recommend not porting, and deleting the flags rather than carrying dead paths. Your call — say so - if they should survive. **Phase 2 has already acted on the recommendation**: R's `col_data` - carries `batch_factor` and `replicate_factor`, and `sim_input.h5` carries neither, since - `--cell-batches` is the only thing that reads them. They are recoverable without touching the - export if that changes — the design matrix holds them one-hot (`batch_factorBatch 2/3/4` plus an - all-zero reference level), so reconstructing the factor is about fifteen lines, and the container - already stores categoricals as codes plus levels for exactly this. +- **`--n-control-cells` and `--cell-batches`. Decided: neither is ported.** + + `--n-control-cells` draws a fixed number of control cells per target instead of using every + non-perturbed cell — on day0, 5,000 in place of ~567,000. It is a cost lever and nothing else, + and it was measured to cost **21–60 % of power**. `--cell-batches` stratifies that draw so the + sampled controls keep the perturbed cells' batch composition; it does nothing on its own, and + `run_power_simulation.R` refuses it without `--n-control-cells`. + + The question worth answering, because it is the one that sounds alarming: **does dropping + `--cell-batches` expose the analysis to batch drift between the two arms?** No, for two reasons. + + 1. With no subsampling there is no draw to stratify. The control group is every non-perturbed + cell, so its batch composition is the dataset's, not a sampling artefact. + 2. Batch is conditioned on by the test itself, twice. `batch_factorBatch 2/3/4` and + `replicate_factorRep 2/3/4` are columns of the covariate matrix, and that matrix enters both + the per-gene NB fit — so batch effects on expression are adjusted out — and the logistic fit + of perturbation status that the **CRT draws its synthetic treated sets from**. The null + distribution is therefore conditional on batch by construction. That is the formal guarantee, + and it is why this method does not need matched control cells: cell-level matching is what + you reach for when the model cannot adjust for a confounder, and here it can. + + Stratified sampling was never the defence against batch confounding. It was a patch for the + extra variance that careless subsampling adds on top of a model already handling it. + + The case against porting is stronger here than it was in R: the lever exists to buy speed, the + port is the reason speed stops being the binding constraint, and a knob that trades power for + speed you no longer need is a trap rather than an option. **`perturbation.py` should not grow a + control-sampling path**, and `sim_input.h5` carries no `batch_factor` or `replicate_factor`, + since `--cell-batches` was the only reader. + + **If control subsampling ever returns** — a screen large enough that even the Python path cannot + afford the full control set — stratification has to return with it, and the reasoning above is + why. That costs no format change: the design matrix holds both factors one-hot (`batch_factorBatch + 2/3/4` plus an all-zero reference level), so each is reconstructible in about fifteen lines, and + `sim_input` already stores categoricals as codes plus levels for exactly this. + + The R implementation keeps both flags. It is the reference the paper describes, and removing + options from it would change what that reference is. - **`run_permutations = TRUE` screens.** pysceptre supports permutations, but its draws are sized by the largest target in the run, which interacts badly with per-target calls. Refuse for now. From db520a4139717c0e4e0136f167b230a76368af28 Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Mon, 21 Sep 2026 13:07:26 -0400 Subject: [PATCH 13/83] Record that the simulated mean is biased low, and by how much Prompted by a report that the simulation runs 16% low. It does not: that is the gap between row_data$mean and the raw mean, which is what a normalised mean is for, and draw_counts multiplies the size factor back in. Measured on day0, the simulated per-gene expected raw mean is 0.959 of the real one at the median, not 0.84. The residual is real though. row_data$mean is a mean of ratios, and raw = E[x]E[sf] + Cov(x, sf), so multiplying by mean(sf) drops a covariance term that is positive here -- the normalisation under-corrects, the simulation draws low, and simulated power is biased low with it. A ratio of sums, rowSums(counts) / sum(sf), makes the simulated total equal the observed total identically for every gene. Recorded rather than applied: it moves every number the paper reports, and reproducing R is what the phase-2 gate is for. Co-Authored-By: Claude Opus 5 (1M context) --- docs/pysceptre-backend.md | 34 ++++++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) diff --git a/docs/pysceptre-backend.md b/docs/pysceptre-backend.md index 5cb5f7e..471b3a9 100644 --- a/docs/pysceptre-backend.md +++ b/docs/pysceptre-backend.md @@ -165,6 +165,40 @@ port achieves and still tight enough to catch a wrong model. joined by one R never needed: a theta clamped to the estimator's bounds is refused rather than simulated from, because here the fit happens in front of us rather than arriving in a cache. +### 3.2b The simulated mean is biased low, and the fix is exact + +**Found while checking a report that the simulation runs 16 % low.** It does not, but it does run +low. `row_data$mean` is 16 % below the raw mean on day0 — that much is true and expected, because +it is a *normalised* mean — but the simulation never uses it alone: `draw_counts` forms +`mu[i,j] = mean_i x size_factor_j x effect_size`, which puts the size factor back. With +`mean(sf) = 1.1389` on day0 the simulated per-gene expected raw mean lands at **0.959 of the real +one** at the median (range 0.89–1.04; 84 of 237 genes more than 5 % low, one more than 10 %). + +The residual is a real bias with a clean cause. `mean_i` is a **mean of ratios**, +`(1/n) sum_j counts[i,j]/sf_j`, and + + raw_mean = E[x.sf] = E[x].E[sf] + Cov(x, sf), x = counts/sf + +so multiplying by `mean(sf)` drops the covariance term. It is positive here — cells with larger +size factors still carry slightly more normalised counts, i.e. the normalisation under-corrects — +so the simulation draws low, which biases simulated power **low**. The published sweep is +conservative by roughly that much. + +**The alternative is a ratio of sums, and it is exact rather than better:** + + mean_i = rowSums(counts)_i / sum_j sf_j + +Then `sum_j mean_i.sf_j = sum_j counts[i,j]` identically, for every gene, by construction. On day0 +it raises each gene's mean by ~4.3 %. + +**Not changed, pending a decision.** It moves every number the paper reports, and the port's job +is to reproduce R first — that is what the phase-2 gate is. It is a better candidate for actually +changing than the §5.2 question, though: that one is a judgement about which cells belong in an +estimator, this one is an estimator that provably fails to reproduce the counts it was derived +from, with an exact alternative available. If it is taken, it is a one-line change in +`watteg/expression.py` and a corresponding one in `src/prepare_sim_input.R`, and both sweeps +have to be re-run. + ### 3.3 Seeding contract is preserved Today: `set.seed(derive_seed(seed, target, rep, effect_size))` before each replicate, and a separate From 04c76460641478046d6a2d802da92d891d77ebd0 Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Mon, 21 Sep 2026 14:13:13 -0400 Subject: [PATCH 14/83] Say plainly that sceptre is not involved in the mean bias The obvious question about section 3.2b is whether sceptre has the same problem. It does not: it computes no size factor, geometric mean or normalised mean -- none of its 183 functions match a search for them -- and its per-gene model is a Poisson GLM with no offset, taking library size as ordinary covariates. The size factors and row_data$mean are WattEG's simulation machinery, inherited from the original DC_TAP_Paper code. DESeq2 is not wrong either; mean-of-ratios is what baseMean is. Only the composition is WattEG's own. Reading that function also settled the theta clamp: sceptre clamps to exactly the [0.01, 1000] pysceptre uses and carries on, so refusing a clamped theta is a deliberate deviation rather than a check R never needed. Section 3.2 now says so, and why a simulation should refuse where an analysis may proceed. Co-Authored-By: Claude Opus 5 (1M context) --- docs/pysceptre-backend.md | 19 +++++++++++++++++-- 1 file changed, 17 insertions(+), 2 deletions(-) diff --git a/docs/pysceptre-backend.md b/docs/pysceptre-backend.md index 471b3a9..2f0396e 100644 --- a/docs/pysceptre-backend.md +++ b/docs/pysceptre-backend.md @@ -162,11 +162,26 @@ both implementations alike, and it is the gene the two differ on most. `workflow/compare_dispersion.py` is the gate, at a tolerance of 1e-6 — generous against what the port achieves and still tight enough to catch a wrong model. `build_dispersion_vector`'s hard error on missing/non-finite dispersions carries over, and is -joined by one R never needed: a theta clamped to the estimator's bounds is refused rather than -simulated from, because here the fit happens in front of us rather than arriving in a cache. +joined by one **deliberate deviation**: a theta clamped to the estimator's bounds is refused rather +than simulated from. sceptre clamps to exactly the same `[0.01, 1000]` that pysceptre does +(`perform_response_precomputation`: `max(min(theta, 1000), 0.01)`) and carries on, so R would +proceed where this stops. For an *analysis* that is reasonable; for a *simulation* a clamped theta +is not an estimate of anything, and drawing counts from it would state a noise level the data never +supported. Day0 has no clamped gene, so nothing is refused today. ### 3.2b The simulated mean is biased low, and the fix is exact +**sceptre is not involved, and this is the first thing to establish.** sceptre computes no size +factor, no geometric mean and no normalised mean — searching all 183 of its functions for +`size_factor|geomean|normaliz|offset` returns nothing — and its per-gene model is +`glm.fit(y = counts, x = covariate_matrix, family = poisson())` with no offset. Library size enters +it as ordinary covariates, `log(response_n_umis)` and `log(response_n_nonzero)`, which the GLM fits +coefficients for. The size factors and `row_data$mean` are **WattEG's simulation machinery alone**, +inherited from the original DC_TAP_Paper power simulation; they exist to generate synthetic counts +and sceptre never sees them. DESeq2 is not doing it wrong either — mean-of-ratios is exactly what +its `baseMean` is, and as a summary statistic it is fine. What follows is about one *composition*, +which is WattEG's own. + **Found while checking a report that the simulation runs 16 % low.** It does not, but it does run low. `row_data$mean` is 16 % below the raw mean on day0 — that much is true and expected, because it is a *normalised* mean — but the simulation never uses it alone: `draw_counts` forms From f8ab01b5127fd4dc49eb21acbb44d739079339a7 Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Mon, 21 Sep 2026 14:16:26 -0400 Subject: [PATCH 15/83] Connect the mean question to the analytical side The question was whether PerturbPlan has the same problem. It does not: it takes expression_mean as an input, and the error is in feeding it the poscounts normalised mean. pysceptre's inputs.py already says so and recommends baseline_expression_stats_from_fits. The two tools differ only in how far that input travels. The analytical formula uses it as-is, so the full 16% lands in the answer; WattEG multiplies the size factor back per cell, so 4.1% survives. "The mean sceptre's model implies" is the observed raw mean exactly -- mean(exp(X.beta)) reproduces it to 8e-9, since a Poisson GLM with an intercept satisfies sum(fitted) == sum(observed). That identity exposes something the level difference hides: sceptre's mean is covariate-dependent per cell and WattEG's is one scalar per cell, so they differ in shape and not only in scale. Simulating from exp(X.beta) directly is the same move the analytical side already made, and it would delete the poscounts machinery along with both open questions that exist only because of it. Recorded with the argument against, since it decides what the power analysis means rather than fixing a bug. Co-Authored-By: Claude Opus 5 (1M context) --- docs/pysceptre-backend.md | 42 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 42 insertions(+) diff --git a/docs/pysceptre-backend.md b/docs/pysceptre-backend.md index 2f0396e..76b23ea 100644 --- a/docs/pysceptre-backend.md +++ b/docs/pysceptre-backend.md @@ -214,6 +214,48 @@ from, with an exact alternative available. If it is taken, it is a one-line chan `watteg/expression.py` and a corresponding one in `src/prepare_sim_input.R`, and both sweeps have to be re-run. +### 3.2c The same question on the analytical side, and what it suggests + +pysceptre's analytical power estimator hit this first, and its +`analytical_power/inputs.py` already says so: `baseline_expression_stats` "comes from a +normalisation scheme sceptre does not use, and on day0 it sits about 16 % below the mean sceptre's +own model implies", with `baseline_expression_stats_from_fits` recommended instead. + +**PerturbPlan is not doing anything wrong.** It takes `expression_mean` as an input. What is wrong +is feeding it the poscounts normalised mean, and the two tools differ only in how far that input +travels: + +| | vs the mean sceptre's model implies | +|---|---| +| analytical formula — uses `expression_mean` **as-is**, no per-cell factor | **0.84, 16 % low** | +| WattEG's simulation — `mean_i x sf_j`, the factor multiplied back | **0.959, 4.1 % low** | + +"The mean sceptre's model implies" is measurable and is exactly the observed raw mean: +`mean(exp(X.beta))` reproduces it to 8e-9, because a Poisson GLM with an intercept satisfies +`sum(fitted) == sum(observed)`. That identity is what makes this comparison sharp rather than a +matter of taste. + +**The deeper point, which the level difference hides.** sceptre's model is +`E[count_ij] = exp(X_j . beta_i)`, a full covariate-dependent mean. WattEG simulates +`E[count_ij] = mean_i x sf_j`, one scalar per cell. Those differ in *shape across cells*, not only +in scale, so the ratio-of-sums fix in §3.2b patches the average and leaves the structure wrong. + +**So the candidate fix is bigger than §3.2b and subsumes it: simulate from `exp(X_j . beta_i)` +directly**, which is the same move the analytical side already made. It reproduces the observed +mean exactly rather than approximately, reproduces the per-cell variation the test's own model +assumes, and deletes the poscounts machinery entirely — with it go both §5.2's QC-cells question +and §3.2b's estimator question, which exist only because size factors do. `fit_dispersions` +already computes `fitted_coefs` and currently discards them. + +**The argument against, stated because it is real.** Simulating from the fitted model and then +testing with that same model makes the test perfectly specified by construction, which may +overstate power slightly; the present scheme is misspecified in the other direction. Neither is +neutral. "Matches the model the real data was fit with" is the more defensible starting point, but +this is a decision about what the power analysis *means*, not a bug fix, and it is the user's. + +**Nothing is changed pending that decision**, and the phase-2 gate is unaffected either way: it +asks whether the port reproduces R, and it does. + ### 3.3 Seeding contract is preserved Today: `set.seed(derive_seed(seed, target, rep, effect_size))` before each replicate, and a separate From 90308fb5bc3d52bf0ef4eb3d19d1fc35f581f33d Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Mon, 21 Sep 2026 14:18:08 -0400 Subject: [PATCH 16/83] Record what PerturbPlan's post-hoc source actually requires compute_power_posthoc does not compute expression_mean -- it takes the frame as an argument, and the documentation defines the columns without ever stating the scale. The formula states it twice: compute_distribution_teststat uses var_nb = mean + mean^2/size, the variance of the NB the observed counts follow, and compute_QC takes P(count == 0) from that NB and asks pbinom how many cells clear the n_nonzero thresholds. Only a distribution over actual counts has a zero probability, so expression_mean is E[observed count]. That also shows the wrong scale is counted twice on the analytical side. On day0 the normalised mean overstates P(count == 0) by 0.032 at the median, which inflates QC_prob, and power is scaled by 1 - QC_prob: 28 of 237 genes at the median target, the worst losing 0.20 of its power to a QC term it should not have had. Co-Authored-By: Claude Opus 5 (1M context) --- docs/pysceptre-backend.md | 24 +++++++++++++++++++++--- 1 file changed, 21 insertions(+), 3 deletions(-) diff --git a/docs/pysceptre-backend.md b/docs/pysceptre-backend.md index 76b23ea..c6539a9 100644 --- a/docs/pysceptre-backend.md +++ b/docs/pysceptre-backend.md @@ -221,9 +221,27 @@ pysceptre's analytical power estimator hit this first, and its normalisation scheme sceptre does not use, and on day0 it sits about 16 % below the mean sceptre's own model implies", with `baseline_expression_stats_from_fits` recommended instead. -**PerturbPlan is not doing anything wrong.** It takes `expression_mean` as an input. What is wrong -is feeding it the poscounts normalised mean, and the two tools differ only in how far that input -travels: +**PerturbPlan is not doing anything wrong, and its source says so twice.** `compute_power_posthoc` +does not compute `expression_mean`; it takes `baseline_expression_stats` as an argument, and the +documentation defines it only as "a data frame ... with columns `response_id`, `expression_mean`, +and `expression_size`" — **it never states the scale.** The formula does, in two independent +places: `compute_distribution_teststat` uses `var_nb(mean, size) = mean + mean^2/size`, the +variance of the negative binomial the *observed counts* follow; and `compute_QC` takes +`P(count == 0)` from that same NB and feeds it to +`pbinom(n_nonzero_thresh - 1, num_cells, 1 - P0)` — the chance that enough cells have a nonzero +**raw count** to clear pairwise QC. The second settles it: only a distribution over actual counts +has a zero probability to ask about. So `expression_mean` must be E[observed count per cell], and +the defect is entirely in the input. + +**And the wrong scale costs more there than the flat 16 % suggests, because it is counted twice.** +Measured on day0, the normalised mean overstates `P(count == 0)` by 0.032 at the median and 0.073 +at most (0.304 against 0.247). That inflates `QC_prob`, and power is multiplied by +`1 - QC_prob`: at the median target's 396 treated cells, **28 of 237 genes carry an inflated +`QC_prob`, the worst by 0.20** — a fifth of that gene's power disappearing into a QC term, on top +of the separate understatement through the test statistic. + +What is wrong, then, is feeding it the poscounts normalised mean, and the two tools differ only in +how far that input travels: | | vs the mean sceptre's model implies | |---|---| From 8aa55aab68554d0ef89e40dea3cb04acd52b5956 Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Mon, 21 Sep 2026 14:25:31 -0400 Subject: [PATCH 17/83] Gate the fitted coefficients, not just theta The simulation is about to draw its baseline from exp(X . beta) rather than from a normalised mean, which makes the coefficients load-bearing: they set the expected count of every cell. Theta already had a gate; these did not. Reported per coefficient column rather than pooled, because the eleven live on different scales and a pooled maximum would let the factor dummies disagree while the continuous terms carried the summary. A contrast-ordering mismatch is the failure that looks like nothing at all -- it would shift the simulated baseline by batch, quietly, for every gene. All eleven pass against sceptre's cache, dummies included, and so does the quantity they exist for: exp(X . beta) over all 134,542,530 gene x cell values agrees to 9.2e-10 at worst and 1.3e-11 at the median. Co-Authored-By: Claude Opus 5 (1M context) --- workflow/compare_gene_model.py | 121 +++++++++++++++++++++++++++++++++ 1 file changed, 121 insertions(+) create mode 100755 workflow/compare_gene_model.py diff --git a/workflow/compare_gene_model.py b/workflow/compare_gene_model.py new file mode 100755 index 0000000..822bac0 --- /dev/null +++ b/workflow/compare_gene_model.py @@ -0,0 +1,121 @@ +#!/usr/bin/env python3 +"""Phase 2's third gate: the per-gene model, coefficients as well as theta. + + workflow/compare_gene_model.py + +`compare_dispersion.py` checked theta. Once the simulation draws its baseline +from `exp(X . beta)` rather than from a normalised mean, **the coefficients are +load-bearing too** -- they set the expected count of every cell -- so they need +a gate of their own rather than a spot check. + +**Reported per coefficient column, not pooled.** The eleven coefficients live on +different scales and mean different things: an intercept, four log-covariate +slopes, and six factor dummies. A pooled maximum would let the batch dummies +disagree while the continuous terms carried the summary, and a factor-contrast +ordering mismatch is exactly the failure that looks like nothing at all -- +it would shift the simulated baseline **by batch**, quietly, for every gene. +""" + +from __future__ import annotations + +import argparse +import sys +from pathlib import Path + +import numpy as np + +from watteg.dispersion import fit_dispersions + +FAILURES: list[str] = [] + + +def main() -> int: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("export", type=Path) + parser.add_argument("reference", type=Path, help="output of workflow/dump_r_sim_input.R") + parser.add_argument("--rtol", type=float, default=1e-6) + parser.add_argument("--n-jobs", type=int, default=1) + args = parser.parse_args() + + sys.path.insert(0, str(Path(__file__).resolve().parents[2] / "pysceptre" / "scripts")) + from pysceptre.pipeline.discovery import ( + _GENE_BATCH_WIDTH, + _GENE_CHUNK_MEMORY_GB, + fit_all_genes, + ) + from sceptre_io import load_export + + # cells_in_use and the export's covariate matrix: the same cells and the + # same design sceptre's precomputation used. A different design silently + # answers a different question, which is why this is not parameterised. + export = load_export(args.export) + genes = (args.reference / "fitted_coefs_genes.txt").read_text().split() + names = (args.reference / "fitted_coefs_names.txt").read_text().splitlines() + theirs = np.array( + [ + [float(v) for v in line.split("\t")] + for line in (args.reference / "fitted_coefs_full.txt").read_text().splitlines() + ] + ) + + if list(export.metadata["covariate_names"]) != names: + print("FAIL covariate columns differ between the export and R's cache") + print(f" export: {list(export.metadata['covariate_names'])}") + print(f" R: {names}") + return 1 + print(f"covariate columns agree, in order ({len(names)})") + + idx = {g: i for i, g in enumerate(export.gene_ids)} + fits = fit_all_genes( + export.response_matrix, + genes, + export.covariate_matrix, + chunk_memory_gb=_GENE_CHUNK_MEMORY_GB, + gene_rows=[idx[g] for g in genes], + batch_width=_GENE_BATCH_WIDTH, + n_jobs=args.n_jobs, + ) + ours = np.array([fits[g].fitted_coefs for g in genes]) + print(f"fitted {ours.shape[0]} genes x {ours.shape[1]} coefficients\n") + + for j, name in enumerate(names): + a, b = ours[:, j], theirs[:, j] + rel = np.abs(a - b) / np.maximum(np.abs(b), np.finfo(float).tiny) + ok = bool(np.all(rel < args.rtol)) + if not ok: + FAILURES.append(name) + print( + f"{'PASS' if ok else 'FAIL'} {name:<28} max rel {rel.max():.2e} " + f"median {np.median(rel):.2e}" + ) + + # What the coefficients are actually for: the per-cell expected count. A + # per-coefficient tolerance can pass while the fitted values drift, because + # the design columns are correlated, so the quantity the simulation uses is + # checked directly as well. + X = export.covariate_matrix + mu_ours = np.exp(X @ ours.T) + mu_theirs = np.exp(X @ theirs.T) + rel = np.abs(mu_ours - mu_theirs) / mu_theirs + ok = bool(np.all(rel < args.rtol)) + if not ok: + FAILURES.append("fitted values") + print( + f"\n{'PASS' if ok else 'FAIL'} exp(X . beta), every gene x cell " + f"max rel {rel.max():.2e} median {np.median(rel):.2e} " + f"({rel.size:,} values)" + ) + + _, diagnostics = fit_dispersions( + export.response_matrix, export.gene_ids, X, genes, n_jobs=args.n_jobs + ) + for kind, affected in diagnostics.items(): + if affected: + print(f" {kind}: {len(affected)} gene(s) -- {', '.join(affected[:5])}") + + print("\n" + ("ALL PASS" if not FAILURES else f"{len(FAILURES)} FAILED: {FAILURES}")) + return 1 if FAILURES else 0 + + +if __name__ == "__main__": + raise SystemExit(main()) From d4fe805e3290ae8b61aa52e8742e448fb1df6250 Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Mon, 21 Sep 2026 14:31:57 -0400 Subject: [PATCH 18/83] Simulate from sceptre's own model, not from a normalised mean The simulation drew counts from mean_i * sf_j: a DESeq2-normalised gene mean scaled by the cell's poscounts size factor. It now draws from exp(X . beta), the expected count sceptre's own null model gives that cell for that gene, so the simulation and the test that judges it are on one scale by construction. Three things were wrong with the old baseline, in increasing order of weight. It mixed two models: dispersion from sceptre's negative binomial, level from a DESeq2 normalisation, for the same gene. It got the level wrong. row_data$mean sits 16% below the mean sceptre's model implies; multiplying by the size factor recovers most of it and leaves the simulated genes about 4% low. The residual is a dropped covariance term -- mean_i is a mean of ratios, and E[x.sf] = E[x]E[sf] + Cov(x, sf). It got the SHAPE wrong, which the level hides. sceptre's mean varies with every covariate; mean_i * sf_j varies with one scalar per cell. Measured against the real counts over 60 genes and 567,690 cells, exp(X . beta) predicts the zero fraction to 0.0007 against 0.0053 and 99.5% of the observed count variance against 86.5%. Note the variance error pulls the OPPOSITE way from the level error -- less variance inflates power, less expression deflates it -- so the effect on simulated power is not obviously signed and is not claimed here. Stage B under both baselines is what measures it. size_factor survives as a --expression-model choice, for one reason: Stage B compares against sweeps produced with it, and that is only interpretable like for like. It is a validation fixture, not a modelling option. sim_input is format 2, carrying fitted_coefs -- eleven numbers per gene rather than a baseline matrix, which is 567,690 times smaller -- and refusing v1 with a message naming the re-run. All four phase-2 gates stay green: the change is additive, and everything the port already reproduced is untouched. Co-Authored-By: Claude Opus 5 (1M context) --- src/watteg/baseline.py | 104 +++++++++++++++++++++++++ src/watteg/cli/prepare_sim_input.py | 11 +-- src/watteg/dispersion.py | 90 --------------------- src/watteg/gene_model.py | 117 ++++++++++++++++++++++++++++ src/watteg/sim_input.py | 30 +++++-- workflow/compare_dispersion.py | 7 +- workflow/compare_gene_model.py | 7 +- 7 files changed, 258 insertions(+), 108 deletions(-) create mode 100644 src/watteg/baseline.py delete mode 100644 src/watteg/dispersion.py create mode 100644 src/watteg/gene_model.py diff --git a/src/watteg/baseline.py b/src/watteg/baseline.py new file mode 100644 index 0000000..aff9da7 --- /dev/null +++ b/src/watteg/baseline.py @@ -0,0 +1,104 @@ +"""A gene's expected counts per cell, before any perturbation. + +This is the number the simulation multiplies the effect size into, so it decides +what "unperturbed" means and therefore what power is measured against. There are +two ways to produce it and they are not equivalent. + +**`fitted` (the default): `exp(X_j . beta_i)`.** The expected count sceptre's own +null model gives that cell for that gene. Simulation and test then live on one +scale by construction rather than by coincidence. + +**`size_factor`: `mean_i x sf_j`.** What the R implementation does -- a +size-factor-normalised gene mean scaled by the cell's DESeq2 "poscounts" factor. +Kept because the published sweeps were produced with it, and a comparison +against them is only interpretable like for like. It is a validation fixture, +not a modelling choice anyone should make afresh. + +## Why the default changed + +Three reasons, in increasing order of how much they matter. + +**It mixes two models.** The R scheme takes its dispersion from sceptre's +negative-binomial fit and its level from a DESeq2 normalisation, so a simulated +gene's noise and its expression come from different statistical models of the +same data. + +**It gets the level wrong, and the error survives the size factor.** +`row_data$mean` sits 16 % below the mean sceptre's model implies; multiplying by +the cell's size factor recovers most of that, leaving the simulated genes about +4 % low on day0. The residual is a dropped covariance term -- `mean_i` is a mean +of ratios, and `E[x.sf] = E[x].E[sf] + Cov(x, sf)`. + +**It gets the shape wrong, which the level hides.** sceptre's mean varies with +every covariate -- library size, detected genes, batch, replicate -- while +`mean_i x sf_j` varies with a single scalar per cell. Measured against the real +counts on day0, over 60 genes and 567,690 cells: + +| | `mean_i x sf_j` | `exp(X . beta)` | +|---|---|---| +| zero fraction, mean absolute error | 0.0053 | **0.0007** | +| predicted variance / observed | 0.865 | **0.995** | + +The current scheme understates the count variance by 13.5 %. Note that pulls the +**opposite** way from the level error -- less variance inflates power, less +expression deflates it -- so the effect of the change on simulated power is not +obviously signed and has not been measured. Stage B under both modes is what +measures it. +""" + +from __future__ import annotations + +import numpy as np + +MODELS = ("fitted", "size_factor") + + +def fitted_baseline(fitted_coefs: np.ndarray, covariate_matrix: np.ndarray) -> np.ndarray: + """`exp(X . beta)` for one gene or many. + + `fitted_coefs` is `(p,)` for a single gene or `(n_genes, p)`; the result is + `(n_cells,)` or `(n_genes, n_cells)` to match. One matrix product, so the + coefficients are what a `sim_input` stores rather than the baseline itself: + 11 numbers per gene against one per cell per gene. + """ + coefs = np.asarray(fitted_coefs, dtype=float) + single = coefs.ndim == 1 + eta = covariate_matrix @ (coefs if single else coefs.T) + return np.exp(eta if single else eta.T) + + +def size_factor_baseline(mean: np.ndarray, size_factors: np.ndarray) -> np.ndarray: + """`mean_i x sf_j`, the R implementation's baseline. + + `mean` is `(n_genes,)` or a scalar; the result is `(n_genes, n_cells)` or + `(n_cells,)`. + """ + mean = np.asarray(mean, dtype=float) + if mean.ndim == 0: + return float(mean) * np.asarray(size_factors, dtype=float) + return np.outer(mean, size_factors) + + +def baseline_expression( + model: str, + *, + fitted_coefs: np.ndarray | None = None, + covariate_matrix: np.ndarray | None = None, + mean: np.ndarray | None = None, + size_factors: np.ndarray | None = None, +) -> np.ndarray: + """Dispatch on the model name, refusing the arguments the other one needs. + + Keyword-only and explicit about what is missing: passing `mean` to the + fitted model, or coefficients to the size-factor model, is a mistake that + would otherwise be silent in a pipeline where both are available. + """ + if model not in MODELS: + raise ValueError(f"model must be one of {list(MODELS)}, got {model!r}") + if model == "fitted": + if fitted_coefs is None or covariate_matrix is None: + raise ValueError("the 'fitted' baseline needs fitted_coefs and covariate_matrix") + return fitted_baseline(fitted_coefs, covariate_matrix) + if mean is None or size_factors is None: + raise ValueError("the 'size_factor' baseline needs mean and size_factors") + return size_factor_baseline(mean, size_factors) diff --git a/src/watteg/cli/prepare_sim_input.py b/src/watteg/cli/prepare_sim_input.py index 4b2f7e5..b86e167 100644 --- a/src/watteg/cli/prepare_sim_input.py +++ b/src/watteg/cli/prepare_sim_input.py @@ -37,8 +37,8 @@ import pandas as pd from scipy import sparse -from watteg.dispersion import fit_dispersions from watteg.expression import compute_expression_stats +from watteg.gene_model import fit_gene_models from watteg.sim_input import SimInput, write_sim_input @@ -151,8 +151,8 @@ def main(argv: list[str] | None = None) -> int: print(f" {len(pairs):,} QC-passing pairs across {pairs['grna_target'].nunique():,} targets") print(f" genes kept: {len(genes)} of {len(in_use.gene_ids)} (those in QC-passing pairs)") - print(f"fitting dispersions for {len(genes)} genes over {n_in_use:,} cells ...") - dispersion, fit_diagnostics = fit_dispersions( + print(f"fitting the per-gene model for {len(genes)} genes over {n_in_use:,} cells ...") + models = fit_gene_models( in_use.response_matrix, in_use.gene_ids, in_use.covariate_matrix, @@ -162,7 +162,7 @@ def main(argv: list[str] | None = None) -> int: # Named, not counted. A gene whose theta came from method of moments rather # than the MLE is the first place to look when its power is surprising, and # a count alone does not say which gene to look at. - for kind, affected in fit_diagnostics.items(): + for kind, affected in models.diagnostics.items(): if affected: shown = ", ".join(affected[:5]) more = f" (+{len(affected) - 5} more)" if len(affected) > 5 else "" @@ -172,7 +172,7 @@ def main(argv: list[str] | None = None) -> int: row_data = pd.DataFrame( { "mean": stats.normalized_mean[gene_rows], - "dispersion": [dispersion[g] for g in genes], + "dispersion": models.dispersion, "average_expression_all_cells": stats.average_expression_all_cells[gene_rows], }, index=genes, @@ -208,6 +208,7 @@ def main(argv: list[str] | None = None) -> int: col_data=col_data, covariate_matrix=in_use.covariate_matrix, covariate_names=list(meta["covariate_names"]), + fitted_coefs=models.fitted_coefs, grna_ids=grna_ids, grna_perts=indicator_matrix(grna_ids, guide_cells, n_in_use), target_ids=target_ids, diff --git a/src/watteg/dispersion.py b/src/watteg/dispersion.py deleted file mode 100644 index 26d4ffe..0000000 --- a/src/watteg/dispersion.py +++ /dev/null @@ -1,90 +0,0 @@ -"""Per-gene dispersion: the noise the simulation exists to reproduce. - -The simulation draws counts from a negative binomial with `size = 1/dispersion`, -so this number sets how noisy a simulated gene is, and therefore how much of the -effect a test can see. It has to come from the **real** counts: a dispersion -estimated from simulated data would make the simulation agree with itself. - -In R it was read from `sceptre_object@response_precomputations`, a cache sceptre -fills during the real discovery analysis. There is no such cache here and there -does not need to be -- pysceptre fits the same model, so the Python pipeline -estimates theta directly and depends on no slot of anyone's object. - -**Fitted over `cells_in_use`, unlike everything in `expression.py`.** sceptre's -precomputation runs on the cells that passed QC, against the covariate matrix, -so matching it means doing the same. That the two halves of `prepare_sim_input` -use different cell sets is not an oversight in either language: size factors are -a property of the library preparation, which every cell took part in, while a -model fit is a property of the analysis, which only the QC-passing cells enter. -""" - -from __future__ import annotations - -import numpy as np -from pysceptre.pipeline.discovery import ( - _GENE_BATCH_WIDTH, - _GENE_CHUNK_MEMORY_GB, - fit_all_genes, - summarize_gene_fits, -) - - -def fit_dispersions( - counts, - gene_ids: list[str], - covariate_matrix: np.ndarray, - genes: list[str], - *, - n_jobs: int = 1, -) -> tuple[dict[str, float], dict[str, list[str]]]: - """Dispersion (`1/theta`) per gene, and which fits were degenerate. - - The second return value names the genes whose GLM did not converge, whose - theta MLE fell back to method of moments, and so on. pysceptre warns about - the counts; the ids are what someone needs when one gene's power looks - wrong, so they are handed back rather than left in a warning string. - - `counts` is (n_genes, n_cells) over `cells_in_use`, `gene_ids` labels its - rows, and `covariate_matrix` is (n_cells, p) over the same cells. - - The batching arguments are pysceptre's own pipeline defaults rather than - this function's choices: `_GENE_BATCH_WIDTH = 1` makes each fit depend on - that gene alone, which is what keeps a dispersion from shifting when the - gene list changes. - """ - index = {g: i for i, g in enumerate(gene_ids)} - unknown = [g for g in genes if g not in index] - if unknown: - raise KeyError( - f"{len(unknown)} gene(s) are not rows of the count matrix, including: {unknown[:5]}" - ) - - fits = fit_all_genes( - counts, - list(genes), - covariate_matrix, - chunk_memory_gb=_GENE_CHUNK_MEMORY_GB, - gene_rows=[index[g] for g in genes], - batch_width=_GENE_BATCH_WIDTH, - n_jobs=n_jobs, - ) - - # A clamped theta is a fit that did not converge to anything usable, and a - # dispersion of 100 or 0.001 would simulate a gene nothing like the real - # one. R had no equivalent check because it read a cache that sceptre had - # already produced; here the fit happens in front of us, so it is checked. - clamped = [g for g in genes if fits[g].theta_clamped] - if clamped: - raise ValueError( - f"{len(clamped)} gene(s) have a theta clamped to the estimator's bounds, including: " - f"{clamped[:5]}. Their dispersion is not an estimate, and simulating from it would " - "misstate the noise. Drop these genes from the discovery pairs." - ) - nonfinite = [g for g in genes if not np.isfinite(fits[g].theta) or fits[g].theta <= 0] - if nonfinite: - raise ValueError( - f"{len(nonfinite)} gene(s) have a non-finite or non-positive theta, including: " - f"{nonfinite[:5]}." - ) - - return {g: 1.0 / fits[g].theta for g in genes}, summarize_gene_fits(fits) diff --git a/src/watteg/gene_model.py b/src/watteg/gene_model.py new file mode 100644 index 0000000..b11cc60 --- /dev/null +++ b/src/watteg/gene_model.py @@ -0,0 +1,117 @@ +"""The per-gene model the simulation draws from: coefficients and theta. + +sceptre's null model for a gene is a Poisson GLM of its counts on the cell +covariates, plus a separately estimated negative-binomial theta +(`perform_response_precomputation`). Both halves come out of one fit, and this +module reproduces that fit with pysceptre, which agrees with sceptre's own cache +to about ten significant digits on both -- theta to 2.8e-12 at the median, and +`exp(X . beta)` to 9.2e-10 across 134.5M gene x cell values. + +**Both halves are used, and that is the point.** The R implementation took theta +from here and the *mean* from a DESeq2-style normalisation, so the simulated +counts had their dispersion from one model and their level from another. Drawing +from `exp(X . beta)` with this theta puts the simulation on the single scale the +discovery test itself works on. See `baseline.py`. + +**Fitted over `cells_in_use`, and over every cell in them -- including a target's +perturbed cells.** That looks wrong for a baseline and is not: sceptre's own null +model is fitted the same way, on every cell, because the control group is the +complement. A gene with a real strong effect therefore carries a little of it in +its baseline, for every target. That was equally true of the normalised mean this +replaces, and the fit must not be changed to exclude perturbed cells, because the +test it is emulating does not. +""" + +from __future__ import annotations + +from dataclasses import dataclass + +import numpy as np +from pysceptre.pipeline.discovery import ( + _GENE_BATCH_WIDTH, + _GENE_CHUNK_MEMORY_GB, + fit_all_genes, + summarize_gene_fits, +) + + +@dataclass(frozen=True) +class GeneModels: + """One fit per gene, in the order `genes` gives. + + `diagnostics` names the genes whose GLM did not converge, whose theta MLE + fell back to method of moments, and so on. pysceptre warns about the + counts; the ids are what someone needs when one gene's power looks wrong, + so they are carried rather than left in a warning string. + """ + + genes: list[str] + fitted_coefs: np.ndarray # (n_genes, p), aligned to the covariate matrix + theta: np.ndarray # (n_genes,) negative-binomial size + diagnostics: dict[str, list[str]] + + @property + def dispersion(self) -> np.ndarray: + """`1/theta`, which is what the R pipeline called it.""" + return 1.0 / self.theta + + +def fit_gene_models( + counts, + gene_ids: list[str], + covariate_matrix: np.ndarray, + genes: list[str], + *, + n_jobs: int = 1, +) -> GeneModels: + """Fit every gene's Poisson GLM and theta. + + `counts` is (n_genes, n_cells) over `cells_in_use`, `gene_ids` labels its + rows, and `covariate_matrix` is (n_cells, p) over the same cells. + + The batching arguments are pysceptre's own pipeline defaults rather than + this function's choices: `_GENE_BATCH_WIDTH = 1` makes each fit depend on + that gene alone, which is what keeps a dispersion from shifting when the + gene list changes. + """ + index = {g: i for i, g in enumerate(gene_ids)} + unknown = [g for g in genes if g not in index] + if unknown: + raise KeyError( + f"{len(unknown)} gene(s) are not rows of the count matrix, including: {unknown[:5]}" + ) + + fits = fit_all_genes( + counts, + list(genes), + covariate_matrix, + chunk_memory_gb=_GENE_CHUNK_MEMORY_GB, + gene_rows=[index[g] for g in genes], + batch_width=_GENE_BATCH_WIDTH, + n_jobs=n_jobs, + ) + + # A clamped theta is a fit that did not converge to anything usable, and a + # dispersion of 100 or 0.001 would simulate a gene nothing like the real + # one. R had no equivalent check because it read a cache that sceptre had + # already produced; here the fit happens in front of us, so it is checked. + clamped = [g for g in genes if fits[g].theta_clamped] + if clamped: + raise ValueError( + f"{len(clamped)} gene(s) have a theta clamped to the estimator's bounds, including: " + f"{clamped[:5]}. Their dispersion is not an estimate, and simulating from it would " + "misstate the noise. Drop these genes from the discovery pairs." + ) + nonfinite = [g for g in genes if not np.isfinite(fits[g].theta) or fits[g].theta <= 0] + if nonfinite: + raise ValueError( + f"{len(nonfinite)} gene(s) have a non-finite or non-positive theta, including: " + f"{nonfinite[:5]}." + ) + + return GeneModels( + genes=list(genes), + fitted_coefs=np.array([fits[g].fitted_coefs for g in genes]), + theta=np.array([fits[g].theta for g in genes]), + diagnostics=summarize_gene_fits(fits), + ) diff --git a/src/watteg/sim_input.py b/src/watteg/sim_input.py index 31836b0..6ad0c64 100644 --- a/src/watteg/sim_input.py +++ b/src/watteg/sim_input.py @@ -6,10 +6,16 @@ same five things plus the covariate matrix, which in R lived on the sceptre template that the Python path has no equivalent of. -**The count matrix is deliberately absent.** The simulation draws counts from -`row_data.mean` and `row_data.dispersion`, so carrying the real counts through -every parallel task would cost memory and deserialisation time and be read by -nothing. +**The count matrix is deliberately absent.** The simulation draws counts from a +per-gene baseline and dispersion, so carrying the real counts through every +parallel task would cost memory and deserialisation time and be read by nothing. + +**The baseline is stored as coefficients, not as a matrix.** `fitted_coefs` plus +`covariate_matrix` give `exp(X . beta)` for any gene in one matrix product -- +eleven numbers per gene against one per cell per gene, which is 567,690 times +smaller. `row_data.mean` is kept beside it because the size-factor baseline the +R implementation used needs it, and reproducing that is how the two are +compared. See `baseline.py`. **Cells are positions, not barcodes.** R stored 586,309 cell barcodes; pysceptre's export does not carry them at all, and nothing in the Python pipeline needs them @@ -49,7 +55,7 @@ import pandas as pd from scipy import sparse -FORMAT_VERSION = 1 +FORMAT_VERSION = 2 @dataclass @@ -65,6 +71,7 @@ class SimInput: genes: list[str] cells_in_use: np.ndarray # (n_cells,) positions in the original object row_data: pd.DataFrame # index = genes; mean, dispersion, average_expression_all_cells + fitted_coefs: np.ndarray # (n_genes, p), aligned to covariate_names col_data: pd.DataFrame # one row per cell, positionally aligned; size_factors covariate_matrix: np.ndarray # (n_cells, p) covariate_names: list[str] @@ -88,6 +95,13 @@ def validate(self) -> None: raise ValueError(f"row_data is missing the '{column}' column") if not np.isfinite(self.row_data["dispersion"]).all(): raise ValueError("row_data.dispersion has non-finite entries") + if self.fitted_coefs.shape != (len(self.genes), len(self.covariate_names)): + raise ValueError( + f"fitted_coefs is {self.fitted_coefs.shape}, expected " + f"({len(self.genes)}, {len(self.covariate_names)})" + ) + if not np.isfinite(self.fitted_coefs).all(): + raise ValueError("fitted_coefs has non-finite entries") if len(self.col_data) != n_cells: raise ValueError(f"col_data has {len(self.col_data)} rows for {n_cells} cells") if "size_factors" not in self.col_data: @@ -170,6 +184,7 @@ def write_sim_input(sim: SimInput, path: str | Path) -> Path: f.create_dataset("covariate_matrix", data=sim.covariate_matrix, compression="gzip") _write_strings(f, "covariate_names", sim.covariate_names) + f.create_dataset("fitted_coefs", data=sim.fitted_coefs) perts = f.create_group("perts") _write_strings(perts, "grna_ids", sim.grna_ids) @@ -184,7 +199,9 @@ def read_sim_input(path: str | Path) -> SimInput: version = int(f.attrs.get("format_version", 0)) if version != FORMAT_VERSION: raise ValueError( - f"{path} is sim_input format {version}, this build reads {FORMAT_VERSION}" + f"{path} is sim_input format {version}, this build reads {FORMAT_VERSION}. " + "Version 2 added fitted_coefs, which the default baseline needs; re-run " + "watteg-prepare-sim-input to produce one." ) genes = _read_strings(f, "genes") row_data = pd.DataFrame( @@ -205,6 +222,7 @@ def read_sim_input(path: str | Path) -> SimInput: col_data=pd.DataFrame(col_data), covariate_matrix=f["covariate_matrix"][:], covariate_names=_read_strings(f, "covariate_names"), + fitted_coefs=f["fitted_coefs"][:], grna_ids=_read_strings(f["perts"], "grna_ids"), grna_perts=_read_sparse(f["perts"], "grna_perts"), target_ids=_read_strings(f["perts"], "target_ids"), diff --git a/workflow/compare_dispersion.py b/workflow/compare_dispersion.py index cab7ab8..63c3e93 100755 --- a/workflow/compare_dispersion.py +++ b/workflow/compare_dispersion.py @@ -30,7 +30,7 @@ import numpy as np -from watteg.dispersion import fit_dispersions +from watteg.gene_model import fit_gene_models def read_floats(path: Path) -> np.ndarray: @@ -57,14 +57,15 @@ def main() -> int: f"{export.covariate_matrix.shape[1]} covariates ..." ) - ours, diagnostics = fit_dispersions( + models = fit_gene_models( export.response_matrix, export.gene_ids, export.covariate_matrix, genes, n_jobs=args.n_jobs, ) - mine = np.array([ours[g] for g in genes]) + diagnostics = models.diagnostics + mine = models.dispersion theirs = read_floats(args.reference / "row_dispersion.txt") rel = np.abs(mine - theirs) / theirs diff --git a/workflow/compare_gene_model.py b/workflow/compare_gene_model.py index 822bac0..5def06d 100755 --- a/workflow/compare_gene_model.py +++ b/workflow/compare_gene_model.py @@ -24,7 +24,7 @@ import numpy as np -from watteg.dispersion import fit_dispersions +from watteg.gene_model import fit_gene_models FAILURES: list[str] = [] @@ -106,10 +106,9 @@ def main() -> int: f"({rel.size:,} values)" ) - _, diagnostics = fit_dispersions( + for kind, affected in fit_gene_models( export.response_matrix, export.gene_ids, X, genes, n_jobs=args.n_jobs - ) - for kind, affected in diagnostics.items(): + ).diagnostics.items(): if affected: print(f" {kind}: {len(affected)} gene(s) -- {', '.join(affected[:5])}") From 5fbb83926e9162bb7f4cb098a88a18b704f2341c Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Mon, 21 Sep 2026 14:32:55 -0400 Subject: [PATCH 19/83] Stop claiming the sweeps are conservative Section 3.2b inferred from the mean error alone that simulated power is biased low. The variance measurement contradicts the certainty: the same baseline understates the count variance by 13.5%, and less variance inflates power where less expression deflates it. Two errors with opposite signs, neither measured against the other, so the direction is not established and the line asserting it is gone. Section 3.2c now records the decision as taken and implemented, with what it does not settle: the effect on power, the generality beyond day0, and that the R implementation and the published sweeps are unchanged. Co-Authored-By: Claude Opus 5 (1M context) --- docs/pysceptre-backend.md | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/docs/pysceptre-backend.md b/docs/pysceptre-backend.md index c6539a9..34b836f 100644 --- a/docs/pysceptre-backend.md +++ b/docs/pysceptre-backend.md @@ -196,8 +196,10 @@ The residual is a real bias with a clean cause. `mean_i` is a **mean of ratios** so multiplying by `mean(sf)` drops the covariance term. It is positive here — cells with larger size factors still carry slightly more normalised counts, i.e. the normalisation under-corrects — -so the simulation draws low, which biases simulated power **low**. The published sweep is -conservative by roughly that much. +so the simulation draws low. **Whether that makes simulated power low is not established**, and +the obvious inference is unsafe: the same baseline also understates the count *variance* by 13.5 % +(§3.2c), and less variance inflates power where less expression deflates it. Two errors, opposite +signs, neither measured against the other. Stage B under both baselines is the experiment. **The alternative is a ratio of sums, and it is exact rather than better:** @@ -271,8 +273,17 @@ overstate power slightly; the present scheme is misspecified in the other direct neutral. "Matches the model the real data was fit with" is the more defensible starting point, but this is a decision about what the power analysis *means*, not a bug fix, and it is the user's. -**Nothing is changed pending that decision**, and the phase-2 gate is unaffected either way: it -asks whether the port reproduces R, and it does. +**Taken, and implemented.** `exp(X . beta)` is the default baseline in the Python port +(`watteg/baseline.py`), `size_factor` survives as a validation fixture so Stage B can compare +like-for-like against sweeps produced with it, and `sim_input` format 2 carries `fitted_coefs`. +All four phase-2 gates stay green, because the change is additive: everything the port already +reproduced is untouched. + +**What this does *not* settle.** The effect on simulated power is unmeasured and not obviously +signed, for the reason in §3.2b. Every number here is day0; the mechanism is general but the +magnitudes are not. And the R implementation is unchanged, so the published sweeps remain on the +old baseline — re-running them is a separate, all-or-nothing decision, since mixing the two scales +within one analysis would be worse than either alone. ### 3.3 Seeding contract is preserved From 859df71763608e4f198037812b846af10e26e7d7 Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Mon, 21 Sep 2026 14:46:11 -0400 Subject: [PATCH 20/83] Draw from sceptre's own model in R too The Python path moved to exp(X %*% coefs) in d4fe805; this brings the R implementation with it, so the tool is correct in the language it currently ships in rather than only in the one replacing it. sceptre caches both halves of its per-gene fit -- fitted_coefs beside theta -- so prepare_sim_input.R now reads both from @response_precomputations, and sim_input carries the coefficients. Until now the dispersion came from sceptre's negative binomial and the expression level from a DESeq2-style normalisation: one simulated gene, two models of the same data. baseline_expression() holds the choice, and the baseline is computed once per target rather than once per simulation, since it depends on neither the effect size nor the draw. --expression-model size_factor reproduces the old behaviour for the published sweeps, and says in its own help why it should not be chosen afresh. Verified end to end on day0: both models run, the draws differ as they must, and the output columns are unchanged. Seven new tests cover the fitted baseline against a per-gene loop -- so a transpose cannot agree with itself -- the legacy outer product, the guards on a mismatched or missing design, and that draw_counts centres on the baseline it is handed rather than on row_data$mean. The R implementation keeps every flag it had. Nothing here re-runs a sweep or changes a published number; mixing the two scales within one analysis would be worse than either alone, so that remains an all-or-nothing decision. Co-Authored-By: Claude Opus 5 (1M context) --- lib/sim_input.R | 35 +++++++++- lib/simulate.R | 117 +++++++++++++++++++++++++++++---- src/prepare_sim_input.R | 11 +++- src/run_power_simulation.R | 35 +++++++++- tests/testthat/test-simulate.R | 102 ++++++++++++++++++++++++++++ 5 files changed, 282 insertions(+), 18 deletions(-) diff --git a/lib/sim_input.R b/lib/sim_input.R index 8e16ee7..3a98aaf 100644 --- a/lib/sim_input.R +++ b/lib/sim_input.R @@ -10,9 +10,15 @@ ## genes character, gene ids (was rownames(sce)) ## cells character, cell barcodes (was colnames(sce)) ## row_data data.frame, one row per gene, rownames == genes -## mean size-factor-normalised mean, drives simulation +## mean size-factor-normalised mean; drives the simulation +## only under the legacy `size_factor` baseline ## dispersion 1/theta from sceptre's cached precomputations ## average_expression_all_cells raw mean, reported in the output only +## fitted_coefs matrix, genes x covariates, rownames == genes +## sceptre's own Poisson-GLM coefficients for the gene, from +## @response_precomputations. exp(coefs %*% t(covariate_matrix)) is the expected +## count the discovery test's null model gives each cell, and is what the +## simulation draws from. See baseline_expression() in simulate.R. ## col_data data.frame, one row per cell, aligned to `cells` by position (no rownames) ## size_factors poscounts size factors ## factors, e.g. batch_factor, replicate_factor; @@ -35,13 +41,14 @@ suppressPackageStartupMessages(library(Matrix)) PERT_LEVELS <- c("grna_perts", "cre_perts") -new_sim_input <- function(genes, cells, row_data, col_data, perts) { +new_sim_input <- function(genes, cells, row_data, col_data, perts, fitted_coefs = NULL) { x <- list( genes = as.character(genes), cells = as.character(cells), row_data = row_data, col_data = col_data, - perts = perts + perts = perts, + fitted_coefs = fitted_coefs ) class(x) <- "sim_input" validate_sim_input(x) @@ -70,6 +77,22 @@ validate_sim_input <- function(x) { " cells.", call. = FALSE) } + # Optional only so that a sim_input written before the fitted baseline existed still loads; + # run_power_simulation.R refuses to use the fitted model without it rather than falling back + # silently to a different one. + if (!is.null(x$fitted_coefs)) { + if (!is.matrix(x$fitted_coefs)) { + stop("sim_input$fitted_coefs must be a matrix.", call. = FALSE) + } + if (!identical(rownames(x$fitted_coefs), x$genes)) { + stop("sim_input$fitted_coefs rownames must equal sim_input$genes, in the same order.", + call. = FALSE) + } + if (any(!is.finite(x$fitted_coefs))) { + stop("sim_input$fitted_coefs has non-finite entries.", call. = FALSE) + } + } + required_row <- c("mean", "dispersion", "average_expression_all_cells") missing_row <- setdiff(required_row, colnames(x$row_data)) if (length(missing_row) > 0) { @@ -125,6 +148,9 @@ subset_genes <- function(x, genes) { } x$genes <- as.character(genes) x$row_data <- x$row_data[x$genes, , drop = FALSE] + if (!is.null(x$fitted_coefs)) { + x$fitted_coefs <- x$fitted_coefs[x$genes, , drop = FALSE] + } x } @@ -177,6 +203,9 @@ print.sim_input <- function(x, ...) { cat("sim_input:", n_genes(x), "genes x", n_cells(x), "cells\n") cat(" row_data:", paste(colnames(x$row_data), collapse = ", "), "\n") cat(" col_data:", paste(colnames(x$col_data), collapse = ", "), "\n") + cat(" fitted_coefs:", + if (is.null(x$fitted_coefs)) "absent (legacy sim_input)" else + paste(ncol(x$fitted_coefs), "covariates"), "\n") for (level in names(x$perts)) { cat(" perts$", level, ": ", nrow(x$perts[[level]]), " rows\n", sep = "") } diff --git a/lib/simulate.R b/lib/simulate.R index 087c6bf..80410dd 100644 --- a/lib/simulate.R +++ b/lib/simulate.R @@ -13,32 +13,92 @@ suppressPackageStartupMessages(library(Matrix)) +#' A gene's expected counts per cell, before any perturbation. +#' +#' This decides what "unperturbed" means, and therefore what power is measured against. There +#' are two ways to produce it and they are not equivalent. +#' +#' `fitted` (the default) is `exp(X %*% coefs)`: the expected count sceptre's own null model +#' gives that cell for that gene, so the simulation and the test that judges it are on one scale +#' by construction. +#' +#' `size_factor` is `mean_i * sf_j`, which is what this pipeline did until 2026-09-21. It is kept +#' so the published sweeps can be reproduced, and for no other reason. Three things are wrong +#' with it, in increasing order of weight: +#' +#' * It mixes two models -- the dispersion comes from sceptre's negative binomial and the level +#' from a DESeq2-style normalisation, for the same gene. +#' * It gets the level wrong. `row_data$mean` sits 16 % below the mean sceptre's model implies; +#' multiplying by the size factor recovers most of that and leaves the simulated genes about +#' 4 % low on day0. The residual is a dropped covariance term: `mean_i` is a mean of ratios, +#' and `E[x*sf] = E[x]E[sf] + Cov(x, sf)`. +#' * It gets the SHAPE wrong, which the level hides. sceptre's mean varies with every covariate; +#' `mean_i * sf_j` varies with one scalar per cell. Measured against the real counts on day0 +#' over 60 genes and 567,690 cells, the fitted baseline predicts the zero fraction to 0.0007 +#' against 0.0053, and 99.5 % of the observed count variance against 86.5 %. +#' +#' Note the variance error pulls the OPPOSITE way from the level error -- less variance inflates +#' power where less expression deflates it -- so which way the change moves power is not obvious +#' and has not been measured. +#' +#' @param x sim_input, already subset to the genes and cells being simulated +#' @param covariate_matrix cells x covariates, rows aligned to x$cells. Required by `fitted`. +#' @return a dense genes x cells matrix of expected counts +baseline_expression <- function(x, covariate_matrix = NULL, + model = c("fitted", "size_factor")) { + model <- match.arg(model) + if (model == "size_factor") { + # outer() replaces the original matrix(rep(...)) + sweep() pair: same values, one allocation + # instead of three. + return(outer(x$row_data$mean, x$col_data$size_factors)) + } + + if (is.null(x$fitted_coefs)) { + stop("This sim_input carries no fitted_coefs, so the 'fitted' baseline cannot be formed. ", + "Re-run prepare_sim_input.R, or pass --expression-model size_factor to reproduce the ", + "pre-2026-09-21 behaviour.", call. = FALSE) + } + if (is.null(covariate_matrix)) { + stop("The 'fitted' baseline needs the covariate matrix.", call. = FALSE) + } + if (nrow(covariate_matrix) != length(x$cells)) { + stop("covariate_matrix has ", nrow(covariate_matrix), " rows but there are ", + length(x$cells), " cells; they must be aligned.", call. = FALSE) + } + if (ncol(covariate_matrix) != ncol(x$fitted_coefs)) { + stop("covariate_matrix has ", ncol(covariate_matrix), " columns but fitted_coefs has ", + ncol(x$fitted_coefs), "; they come from different designs.", call. = FALSE) + } + # tcrossprod(A, B) is A %*% t(B): (genes x p) against (cells x p) gives genes x cells. + exp(tcrossprod(x$fitted_coefs, covariate_matrix)) +} + #' Simulate counts for one perturbation. #' #' @param x sim_input, already subset to the genes being tested #' @param effect_size_mat genes x cells multiplier, columns in the same order as x$cells +#' @param baseline genes x cells expected counts from baseline_expression(). Computed once per +#' target rather than once per simulation: it does not depend on the effect size or the draw. #' @return a dense genes x cells matrix of counts -draw_counts <- function(x, effect_size_mat) { - gene_means <- x$row_data$mean +draw_counts <- function(x, effect_size_mat, baseline) { gene_dispersions <- x$row_data$dispersion - size_factors <- x$col_data$size_factors - n_gene <- length(gene_means) - n_cell <- length(size_factors) + n_gene <- length(gene_dispersions) + n_cell <- length(x$cells) if (!identical(dim(effect_size_mat), c(n_gene, n_cell))) { stop("effect_size_mat is ", paste(dim(effect_size_mat), collapse = " x "), " but ", n_gene, " x ", n_cell, " was expected.", call. = FALSE) } + if (!identical(dim(baseline), c(n_gene, n_cell))) { + stop("baseline is ", paste(dim(baseline), collapse = " x "), + " but ", n_gene, " x ", n_cell, " was expected.", call. = FALSE) + } - # Cell-to-cell variability. Each cell keeps its own size factor: the effect-size matrix is - # indexed by cell, so shuffling would pair one cell's perturbation status with another cell's - # library size. - # - # mu[i, j] = gene_means[i] * size_factor[j] * effect_size[i, j]. - # outer() replaces the original matrix(rep(...)) + sweep() pair: same values, one allocation - # instead of three. - mu <- outer(gene_means, size_factors) * effect_size_mat + # mu[i, j] = baseline[i, j] * effect_size[i, j]. Each cell keeps its own baseline: the + # effect-size matrix is indexed by cell, so shuffling would pair one cell's perturbation status + # with another cell's expected expression. + mu <- baseline * effect_size_mat # size = theta = 1 / dispersion. mu is consumed column-major, and `size` recycles over the # gene index, which is why gene_dispersions must be exactly n_gene long -- see @@ -128,6 +188,37 @@ build_dispersion_vector <- function(precomputations, genes) { stats::setNames(dispersion, genes) } +#' Build the per-gene coefficient matrix from sceptre's cached precomputations. +#' +#' The sibling of build_dispersion_vector(), reading the other half of the same fit. Taking both +#' from one fit is the point: until 2026-09-21 the dispersion came from sceptre's negative +#' binomial and the expression level from a DESeq2-style normalisation, so a simulated gene's +#' noise and its level came from different models of the same data. +#' +#' @param precomputations sceptre_object@response_precomputations +#' @param genes genes to build the matrix for +#' @return genes x covariates matrix, rownames == genes +build_fitted_coefs_matrix <- function(precomputations, genes) { + missing_genes <- setdiff(genes, names(precomputations)) + if (length(missing_genes) > 0) { + stop(length(missing_genes), " gene(s) have no entry in @response_precomputations and so no ", + "fitted coefficients, including: ", + paste(utils::head(missing_genes, 5), collapse = ", "), + ". Re-run sceptre's precomputation, or drop these genes from the discovery pairs.", + call. = FALSE) + } + + coefs <- do.call(rbind, lapply(precomputations[genes], function(p) p$fitted_coefs)) + rownames(coefs) <- genes + + if (any(!is.finite(coefs))) { + bad <- genes[apply(!is.finite(coefs), 1, any)] + stop(length(bad), " gene(s) have a non-finite fitted coefficient, including: ", + paste(utils::head(bad, 5), collapse = ", "), ".", call. = FALSE) + } + coefs +} + ## GUIDE-LEVEL VARIABILITY ========================================================================= #' Pick one expressed guide at random per cell, from a CSC perturbation matrix. diff --git a/src/prepare_sim_input.R b/src/prepare_sim_input.R index 84fcabd..1678f57 100755 --- a/src/prepare_sim_input.R +++ b/src/prepare_sim_input.R @@ -292,6 +292,14 @@ log_step("Genes kept: ", length(genes), " of ", length(all_genes), # with NULL holes here, which unlist() silently dropped, shifting every later gene's dispersion. dispersion <- build_dispersion_vector(so@response_precomputations, genes) +# The other half of the same fit. Both are read here, before slim_sceptre_object() clears the +# slot: the dispersion sets how noisy a simulated gene is and the coefficients set how expressed +# it is, and taking them from one model is what keeps the simulation on the scale the discovery +# test works on. See baseline_expression() in lib/simulate.R. +fitted_coefs <- build_fitted_coefs_matrix(so@response_precomputations, genes) +log_step("Fitted coefficients: ", nrow(fitted_coefs), " genes x ", ncol(fitted_coefs), + " covariates (", paste(utils::head(colnames(fitted_coefs), 3), collapse = ", "), ", ...)") + gene_idx <- match(genes, all_genes) row_data <- data.frame( mean = stats_out$normalized_mean[gene_idx], @@ -366,7 +374,8 @@ sim <- new_sim_input( cells = cells, row_data = row_data, col_data = col_data, - perts = list(grna_perts = grna_perts, cre_perts = cre_perts) + perts = list(grna_perts = grna_perts, cre_perts = cre_perts), + fitted_coefs = fitted_coefs ) print(sim) diff --git a/src/run_power_simulation.R b/src/run_power_simulation.R index f781396..1a1f3cb 100755 --- a/src/run_power_simulation.R +++ b/src/run_power_simulation.R @@ -69,6 +69,16 @@ option_list <- list( help = paste("Number of reps already covered by earlier chunks. The reported `rep`", "column is rep_offset + 1..reps, keeping it unique across chunks", "[default %default].")), + make_option("--expression-model", type = "character", default = "fitted", + dest = "expression_model", + help = paste("Where a gene's unperturbed expected counts come from.", + "'fitted' (default) uses exp(X %%*%% coefs) from sceptre's own null", + "model, so the simulation and the test that judges it are on one", + "scale. 'size_factor' uses mean_i * sf_j, which is what this pipeline", + "did before 2026-09-21; it mixes two models, runs about 4%% low and", + "reproduces only 86.5%% of the observed count variance against the", + "fitted model's 99.5%%. Keep it only to reproduce published sweeps.", + "See baseline_expression() in lib/simulate.R.")), make_option("--guide-sd", type = "double", default = 0.13, dest = "guide_sd", help = paste("Standard deviation of the per-gRNA effect size around the target", "effect size, i.e. guide-to-guide variability [default %default].", @@ -128,6 +138,10 @@ if (opts$effect_size < 0 || opts$effect_size >= 1) { opts$effect_size, "); 0 is the null arm.", call. = FALSE) } if (opts$reps < 1) stop("--reps must be at least 1.", call. = FALSE) +if (!opts$expression_model %in% c("fitted", "size_factor")) { + stop("--expression-model must be 'fitted' or 'size_factor' (got ", opts$expression_model, ").", + call. = FALSE) +} if (!is.null(opts$cell_batches) && is.null(opts$n_control_cells)) { stop("--cell-batches only applies when --n-control-cells is set.", call. = FALSE) } @@ -221,6 +235,10 @@ targets <- unique(split_pairs$grna_target) log_step("Split covers ", length(targets), " targets / ", nrow(split_pairs), " pairs") log_step("Effect size ", opts$effect_size, " (relative expression ", relative_expression, "), reps ", opts$rep_offset + 1L, "-", opts$rep_offset + opts$reps) +log_step("Expression model: ", opts$expression_model, + if (opts$expression_model == "size_factor") + " (legacy: mixes two models and runs ~4% low -- see lib/simulate.R)" else + " (exp(X %*% coefs), sceptre's own null model)") if (!is.null(opts$n_control_cells)) { log_step("Sampling ", opts$n_control_cells, " control cells per target", if (!is.null(opts$cell_batches)) paste0(" stratified by ", opts$cell_batches) else "") @@ -228,6 +246,10 @@ if (!is.null(opts$n_control_cells)) { ## SIMULATE ======================================================================================== +# Resolved once: matching 586,309 barcodes per target would cost more than the simulation. +# Only read when control cells are sampled, which is the one path that reorders them. +all_cell_names <- rownames(template@covariate_data_frame) + results <- vector("list", length(targets) * opts$reps) result_idx <- 0L target_timings <- numeric(0) @@ -279,6 +301,17 @@ for (target in targets) { restore_cell_order <- order(cell_order(pert_status)) gene_object <- subset_genes(pert_object, target_pairs$response_id) + + # The unperturbed expected counts, computed once per target: they depend on neither the effect + # size nor the draw, so recomputing them inside the rep loop would be the same arithmetic 100 + # times over. The covariate rows are aligned to this object's cells -- identical to the + # template's own order unless control cells were sampled, which reorders and subsets them. + target_covariates <- if (is.null(opts$n_control_cells)) { + template@covariate_matrix + } else { + template@covariate_matrix[match(gene_object$cells, all_cell_names), , drop = FALSE] + } + baseline <- baseline_expression(gene_object, target_covariates, model = opts$expression_model) effect_sizes <- stats::setNames( rep(relative_expression, length(gene_object$genes)), gene_object$genes ) @@ -348,7 +381,7 @@ for (target in targets) { gene_effect_sizes = effect_sizes) es_mat <- es_mat[, restore_cell_order, drop = FALSE] - counts <- draw_counts(gene_object, es_mat) + counts <- draw_counts(gene_object, es_mat, baseline) sceptre_use <- target_template sceptre_use@response_matrix <- list(as_sceptre_response_matrix(counts, report_density = FALSE)) diff --git a/tests/testthat/test-simulate.R b/tests/testthat/test-simulate.R index 565480d..717251d 100644 --- a/tests/testthat/test-simulate.R +++ b/tests/testthat/test-simulate.R @@ -97,3 +97,105 @@ test_that("build_dispersion_vector errors on a non-finite dispersion", { "non-finite dispersion" ) }) + +test_that("baseline_expression's fitted model is exp(X %*% coefs), gene by gene", { + # The quantity the whole simulation rests on: the expected count sceptre's own null model gives + # each cell. Checked against a per-gene loop rather than against a second matrix expression, so + # a transpose error in tcrossprod() cannot agree with itself. + set.seed(11) + n_genes <- 4 + n_cells <- 25 + p <- 3 + coefs <- matrix(rnorm(n_genes * p, sd = 0.3), nrow = n_genes, + dimnames = list(paste0("g", 1:n_genes), NULL)) + covariates <- cbind(1, matrix(rnorm(n_cells * (p - 1)), nrow = n_cells)) + + x <- list(genes = rownames(coefs), cells = paste0("c", seq_len(n_cells)), + fitted_coefs = coefs, + row_data = data.frame(mean = runif(n_genes), dispersion = runif(n_genes), + row.names = rownames(coefs)), + col_data = data.frame(size_factors = runif(n_cells, 0.5, 2))) + + got <- baseline_expression(x, covariates, model = "fitted") + expect_equal(dim(got), c(n_genes, n_cells)) + for (i in seq_len(n_genes)) { + expect_equal(got[i, ], as.vector(exp(covariates %*% coefs[i, ]))) + } +}) + +test_that("baseline_expression's size_factor model is the outer product it always was", { + set.seed(12) + x <- list(genes = c("a", "b"), cells = paste0("c", 1:6), + fitted_coefs = NULL, + row_data = data.frame(mean = c(2, 5), dispersion = c(0.1, 0.2), + row.names = c("a", "b")), + col_data = data.frame(size_factors = c(0.5, 1, 1.5, 2, 2.5, 3))) + + expect_equal(baseline_expression(x, model = "size_factor"), + outer(c(2, 5), c(0.5, 1, 1.5, 2, 2.5, 3))) +}) + +test_that("baseline_expression refuses a mismatched or missing design rather than guessing", { + x <- list(genes = "a", cells = c("c1", "c2"), + fitted_coefs = matrix(0.1, nrow = 1, ncol = 3, dimnames = list("a", NULL)), + row_data = data.frame(mean = 1, dispersion = 0.1, row.names = "a"), + col_data = data.frame(size_factors = c(1, 1))) + + expect_error(baseline_expression(x, model = "fitted"), "needs the covariate matrix") + # Wrong number of cells, and wrong number of covariates: both would otherwise recycle or + # broadcast into a plausible-looking matrix. + expect_error(baseline_expression(x, matrix(0, nrow = 5, ncol = 3), model = "fitted"), + "5 rows but there are 2 cells") + expect_error(baseline_expression(x, matrix(0, nrow = 2, ncol = 2), model = "fitted"), + "different designs") + + no_coefs <- x + no_coefs$fitted_coefs <- NULL + expect_error(baseline_expression(no_coefs, matrix(0, nrow = 2, ncol = 3), model = "fitted"), + "carries no fitted_coefs") +}) + +test_that("draw_counts centres on the baseline it is handed, not on a normalised mean", { + # With the effect size fixed at 1 and a large cell count, the simulated mean has to track the + # baseline. This is what makes the two expression models distinguishable at all. + set.seed(13) + n_cells <- 4000 + baseline <- rbind(rep(5, n_cells), rep(0.5, n_cells)) + x <- list(genes = c("a", "b"), cells = paste0("c", seq_len(n_cells)), + row_data = data.frame(mean = c(999, 999), dispersion = c(0.05, 0.05), + row.names = c("a", "b")), + col_data = data.frame(size_factors = rep(1, n_cells))) + + counts <- draw_counts(x, matrix(1, nrow = 2, ncol = n_cells), baseline) + expect_equal(dim(counts), c(2L, n_cells)) + # row_data$mean is deliberately absurd: if it leaked into the draw this would fail loudly. + expect_equal(unname(rowMeans(counts)), c(5, 0.5), tolerance = 0.05) +}) + +test_that("draw_counts rejects a baseline of the wrong shape", { + x <- list(genes = c("a", "b"), cells = c("c1", "c2", "c3"), + row_data = data.frame(mean = c(1, 1), dispersion = c(0.1, 0.1), + row.names = c("a", "b")), + col_data = data.frame(size_factors = rep(1, 3))) + es <- matrix(1, nrow = 2, ncol = 3) + expect_error(draw_counts(x, es, matrix(1, nrow = 2, ncol = 2)), "baseline is 2 x 2") +}) + +test_that("build_fitted_coefs_matrix reads the other half of the same fit", { + precomps <- list( + g1 = list(fitted_coefs = c(0.1, 0.2), theta = 5), + g2 = list(fitted_coefs = c(0.3, 0.4), theta = 8), + g3 = list(fitted_coefs = c(0.5, 0.6), theta = 2) + ) + got <- build_fitted_coefs_matrix(precomps, c("g3", "g1")) + expect_equal(rownames(got), c("g3", "g1")) + expect_equal(got["g1", ], c(0.1, 0.2)) + expect_equal(got["g3", ], c(0.5, 0.6)) + + # The same guard build_dispersion_vector has, for the same reason: a gene with no entry used to + # become a NULL hole that unlist() dropped, shifting every later gene. + expect_error(build_fitted_coefs_matrix(precomps, c("g1", "absent")), + "no entry in @response_precomputations") + precomps$g2$fitted_coefs <- c(0.3, NA) + expect_error(build_fitted_coefs_matrix(precomps, c("g1", "g2")), "non-finite fitted coefficient") +}) From 87cfc6a21ee3d09ce921f2a8935f1803b0904516 Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Mon, 21 Sep 2026 14:52:06 -0400 Subject: [PATCH 21/83] Carry the per-pair QC counts through to the simulation n_nonzero_trt, n_nonzero_cntrl and pass_qc are reported beside every simulated pair. They are facts about the real data, constant across replicates, and the R implementation reads them off the sceptre template's discovery_pairs_with_info. The export was extended to carry that frame for this purpose and then prepare_sim_input did not write it out, so the Python path had nowhere to get them from. They are the first thing anyone looks at when a pair's power is surprising, and once the sceptre object is gone there is nowhere else to recover them. Co-Authored-By: Claude Opus 5 (1M context) --- src/watteg/cli/prepare_sim_input.py | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/src/watteg/cli/prepare_sim_input.py b/src/watteg/cli/prepare_sim_input.py index b86e167..16edbd2 100644 --- a/src/watteg/cli/prepare_sim_input.py +++ b/src/watteg/cli/prepare_sim_input.py @@ -18,6 +18,7 @@ sim_input.h5 per-gene and per-cell statistics + the perturbation matrices pairs.tsv the QC-passing discovery pairs + pairs_with_info.tsv every discovery pair with its real-data QC counts grna_targets.tsv the gRNA -> target mapping discovery_threshold.txt the nominal p-value a simulated pair has to beat analysis_mode.tsv which test the screen was run under @@ -234,6 +235,22 @@ def main(argv: list[str] | None = None) -> int: # disagree with the screen's own design table for no gain. in_use.grna_target_data_frame.to_csv(args.outdir / "grna_targets.tsv", sep="\t", index=False) + # n_nonzero_trt, n_nonzero_cntrl and pass_qc, which the simulation reports beside each + # pair. They are facts about the REAL data and constant across replicates, so they are + # carried rather than recomputed per draw -- which is also what the R implementation does, + # reading them off the sceptre template's discovery_pairs_with_info. They are the first + # thing anyone looks at when a pair's power is surprising, and there is nowhere else to + # recover them from once the object is gone. + if in_use.discovery_pairs_with_info is not None: + in_use.discovery_pairs_with_info.to_csv( + args.outdir / "pairs_with_info.tsv", sep="\t", index=False + ) + else: + print( + " NOTE: the export carries no discovery_pairs_with_info, so pairs_with_info.tsv is " + "not written and the simulation will have no QC counts to report." + ) + threshold = args.threshold or discovery_threshold(in_use.discovery_result) (args.outdir / "discovery_threshold.txt").write_text(f"{threshold:.17g}\n") @@ -257,7 +274,8 @@ def main(argv: list[str] | None = None) -> int: ) print(f" discovery threshold: {threshold:.6g}") print(f" resampling: {mechanism} (MOI: {moi}, side: {export.side})") - print(f"\nwrote 5 files to {args.outdir}") + written = len(list(args.outdir.glob("*"))) + print(f"\nwrote {written} files to {args.outdir}") return 0 From c25f56c7ef28e1835454be0be97c65cec473372c Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Mon, 21 Sep 2026 14:58:35 -0400 Subject: [PATCH 22/83] Port the simulation core: seeds, guide assignment, the draw The three pieces phase 4 needs whichever shape the benchmark picks, written before the benchmark so it is not timing unchecked code. seeds.py keys every stream on (seed, target, replicate, effect size) and nothing else, so a pair's power does not depend on how many splits a cluster wanted or how many replicates fit in a task. rep 0 is reserved for a target's setup draw, which is made once and reused across replicates. perturbation.py carries the one idea in the model: a target's guides do not all work equally well, so each guide draws its own effect size around the target's and each cell inherits whichever guide it carries. Control cells get guides too -- a cell not perturbed by this target usually carries something else, which draws around 1 with the same spread, and sending them to the no-effect row instead would keep the control arm's mean while removing its variance. Guides per target are read from the design table, which is many-to-many. Choosing one guide per cell is vectorised over cells rather than looped as R did: a cell's guides are one contiguous CSC slice, so it needs no sort, and at 567,690 cells the difference is seconds against minutes. simulate.py draws NB(baseline * effect_size, size = theta) and returns int16, refusing overflow rather than wrapping. Holding counts as float64 would cost four times the memory for no information, and how many replicates fit in one process is what the benchmark is about. 16 tests, including that a cell is only ever given a guide it carries, that a shared guide serves both its targets, that each arm centres on what it should average to, and that the draw's variance is mu + mu^2/theta rather than merely its mean being right. Co-Authored-By: Claude Opus 5 (1M context) --- src/watteg/perturbation.py | 200 +++++++++++++++++++++++++++ src/watteg/seeds.py | 55 ++++++++ src/watteg/simulate.py | 69 ++++++++++ tests/test_simulation_core.py | 248 ++++++++++++++++++++++++++++++++++ 4 files changed, 572 insertions(+) create mode 100644 src/watteg/perturbation.py create mode 100644 src/watteg/seeds.py create mode 100644 src/watteg/simulate.py create mode 100644 tests/test_simulation_core.py diff --git a/src/watteg/perturbation.py b/src/watteg/perturbation.py new file mode 100644 index 0000000..9ff59a0 --- /dev/null +++ b/src/watteg/perturbation.py @@ -0,0 +1,200 @@ +"""Which cells a target perturbs, which guide each one carries, and how hard. + +A port of `lib/pert_input.R` and the guide-level half of `lib/simulate.R`. The +modelling content is one idea: **a target's guides do not all work equally +well**, so a simulated knockdown is not a single multiplier applied to every +perturbed cell. Each guide draws its own effect size around the target's, and +each cell inherits the effect of whichever guide it happens to carry. + +That is why the pipeline needs per-gRNA assignments at all, and why the +per-target union sceptre keeps is not enough (see +`docs/pysceptre-backend.md` section 5.1). + +**Control cells carry guides too**, and they matter. A cell not perturbed by +this target is usually perturbed by something else -- another element, or a +non-targeting guide -- and that guide draws an effect size around 1 with the +same guide-to-guide spread. Dropping them would leave every control cell at +exactly 1, which keeps the control arm's mean and removes its variance. + +**Built directly in cell order.** R assembled the guide assignment as +perturbed-cells-then-control-cells and then permuted it back, a legacy of +matching cell barcodes. Nothing here needs that, and the streams differ between +the languages regardless. +""" + +from __future__ import annotations + +from dataclasses import dataclass + +import numpy as np +from scipy import sparse + +NO_GUIDE = 0 + + +@dataclass(frozen=True) +class GuideAssignment: + """One guide index per cell, and what the indices mean. + + `status[j]` is 0 when cell `j` carries no guide, `1..n_target_guides` for + one of this target's guides, and above that for a guide belonging to + something else. The offset is what lets one effect-size table cover both + arms: rows are laid out as [no effect, this target's guides, every other + guide], so a cell's effect size is a single lookup. + """ + + status: np.ndarray # (n_cells,) int + is_perturbed: np.ndarray # (n_cells,) bool + n_target_guides: int + n_other_guides: int + + +def target_cells(cre_perts: sparse.spmatrix, target_ids: list[str], target: str) -> np.ndarray: + """Boolean mask of the cells this target perturbs.""" + try: + row = target_ids.index(target) + except ValueError: + raise KeyError(f"target {target!r} is not a row of cre_perts") from None + return np.asarray(cre_perts[row].todense()).ravel() > 0 + + +def _choose_one_per_cell( + csc: sparse.csc_matrix, + eligible_row: np.ndarray, + cells: np.ndarray, + rng: np.random.Generator, +) -> np.ndarray: + """For each cell in `cells`, one uniformly chosen eligible guide, or -1. + + `csc` is (guides x cells), so a cell's guides are one contiguous slice -- + that is what CSC means, and it is why this needs no sort. Vectorised over + cells rather than looped: R drew one cell at a time, which at 567,690 cells + is the difference between seconds and minutes. + """ + keep = eligible_row[csc.indices] + cell_of_nonzero = np.repeat(np.arange(csc.shape[1]), np.diff(csc.indptr)) + + qualifying_cell = cell_of_nonzero[keep] + qualifying_guide = csc.indices[keep] + per_cell = np.bincount(qualifying_cell, minlength=csc.shape[1]) + # Where each cell's qualifying entries begin in the compacted arrays. + starts = np.concatenate(([0], np.cumsum(per_cell)[:-1])) + + counts = per_cell[cells] + out = np.full(cells.size, -1, dtype=np.int64) + has_any = counts > 0 + if not has_any.any(): + return out + # One uniform per cell, floored into that cell's own count. Drawing for + # every cell and discarding the unused ones keeps the number of draws a + # function of the cell set alone, not of how many happen to carry a guide. + picks = np.floor(rng.random(cells.size) * counts).astype(np.int64) + picks = np.minimum(picks, np.maximum(counts - 1, 0)) # guard the rng.random() == 1 edge + out[has_any] = qualifying_guide[starts[cells[has_any]] + picks[has_any]] + return out + + +def guide_assignment( + grna_csc: sparse.csc_matrix, + grna_ids: list[str], + target_guides: list[str], + is_perturbed: np.ndarray, + rng: np.random.Generator, +) -> GuideAssignment: + """Pick one guide per cell: this target's for perturbed cells, any other's + for control cells. + + `target_guides` comes from the screen's design table, which maps guides to + targets **many-to-many** -- a guide inside two overlapping candidate + elements belongs to both. It must never be derived from a per-unit + annotation that can hold only one target per guide. + + Drawn once per target and reused across replicates: which guide a cell + carries is a fact about the screen, not about a draw. + """ + index = {g: i for i, g in enumerate(grna_ids)} + rows = np.array([index[g] for g in target_guides if g in index], dtype=np.int64) + if rows.size == 0: + raise KeyError( + f"none of the {len(target_guides)} guide(s) given for this target appear in " + "grna_perts, so no perturbed cell can be assigned a guide" + ) + + n_guides = len(grna_ids) + is_target_guide = np.zeros(n_guides, dtype=bool) + is_target_guide[rows] = True + + status = np.zeros(grna_csc.shape[1], dtype=np.int64) + perturbed = np.flatnonzero(is_perturbed) + control = np.flatnonzero(~is_perturbed) + + # Target guides are numbered by their position in `target_guides`, not by + # their row in the matrix, so the effect-size table's first block lines up + # with the guides whose effect sizes it holds. + rank_of_row = np.full(n_guides, -1, dtype=np.int64) + rank_of_row[rows] = np.arange(rows.size) + + chosen = _choose_one_per_cell(grna_csc, is_target_guide, perturbed, rng) + found = chosen >= 0 + status[perturbed[found]] = rank_of_row[chosen[found]] + 1 + + other_rows = np.flatnonzero(~is_target_guide) + rank_of_other = np.full(n_guides, -1, dtype=np.int64) + rank_of_other[other_rows] = np.arange(other_rows.size) + chosen = _choose_one_per_cell(grna_csc, ~is_target_guide, control, rng) + found = chosen >= 0 + status[control[found]] = rows.size + rank_of_other[chosen[found]] + 1 + + return GuideAssignment( + status=status, + is_perturbed=is_perturbed, + n_target_guides=int(rows.size), + n_other_guides=int(other_rows.size), + ) + + +def effect_size_matrix( + assignment: GuideAssignment, + gene_effect_sizes: np.ndarray, + guide_sd: float, + rng: np.random.Generator, +) -> np.ndarray: + """A (genes x cells) multiplier, with each guide's own effect size. + + `gene_effect_sizes` is the *relative expression* the target aims for -- a + 15% knockdown is 0.85 -- one per gene. Each of the target's guides draws + around it, each other guide draws around 1, and negatives clamp to 0 + because a guide cannot make expression negative. + + Clamping biases the mean upward, so both arms are re-centred afterwards: + that is what makes "effect size 0.15" mean a 15% knockdown on average + rather than something slightly weaker. R did this in a separate + `center_effect_size_matrix()`; it is one step here because the two are not + separately meaningful. + """ + gene_effect_sizes = np.asarray(gene_effect_sizes, dtype=float) + n_genes = gene_effect_sizes.size + n_target, n_other = assignment.n_target_guides, assignment.n_other_guides + + # Row 0 is the no-effect row, for cells carrying no guide at all. + table = np.empty((1 + n_target + n_other, n_genes)) + table[0] = 1.0 + table[1 : 1 + n_target] = rng.normal(gene_effect_sizes, guide_sd, size=(n_target, n_genes)) + table[1 + n_target :] = rng.normal(1.0, guide_sd, size=(n_other, n_genes)) + np.clip(table, 0.0, None, out=table) + + matrix = table[assignment.status].T # (genes, cells) + + # Re-centre each arm on what it is supposed to average to. Cells carrying + # no guide sit at exactly 1 and move with their arm, as in R. + for mask, target_mean in ( + (assignment.is_perturbed, gene_effect_sizes), + (~assignment.is_perturbed, 1.0), + ): + if not mask.any(): + continue + block = matrix[:, mask] + shift = target_mean - block.mean(axis=1) + matrix[:, mask] = block + shift[:, None] + np.clip(matrix, 0.0, None, out=matrix) + return matrix diff --git a/src/watteg/seeds.py b/src/watteg/seeds.py new file mode 100644 index 0000000..3d22962 --- /dev/null +++ b/src/watteg/seeds.py @@ -0,0 +1,55 @@ +"""Seeding, so a pair's power does not depend on how the work was split up. + +The unit of work in a sweep is (split, effect size, replicate chunk), and that +partitioning is an operational choice: how many splits a cluster wanted, how +many replicates fit in a task. None of it should reach the numbers. So the +stream is keyed by what the simulation is *about* -- the run's seed, the target, +the replicate and the effect size -- and never by position in a loop. + +That is what makes a run reproducible across layouts: simulating replicates +1-100 in one task and in five tasks of twenty gives the same draws, and so does +moving a target from one split to another. R derived the same four-part key with +`derive_seed()` in `lib/cli.R`; the streams themselves differ between the +languages, and always would, so only the invariance carries over. + +`rep = 0` is reserved for a target's setup -- the per-cell guide assignment, +which is drawn once and reused across replicates. Replicates are numbered from +1, so a setup draw can never collide with a replicate's. +""" + +from __future__ import annotations + +import hashlib + +import numpy as np + +SETUP_REP = 0 + + +def _text_entropy(text: str) -> int: + """A stable 64-bit integer for a string. + + `blake2b`, not Python's `hash()`: that is salted per process, so a target + would seed differently in every task of the same sweep and the invariance + this module exists for would be silently untrue. + """ + return int.from_bytes(hashlib.blake2b(text.encode(), digest_size=8).digest(), "little") + + +def derive_seed(seed: int, target: str, rep: int, effect_size: float) -> np.random.SeedSequence: + """The stream for one (target, replicate, effect size) under a run's seed. + + `effect_size` enters as its 17-digit decimal rather than as a float: the + value arrives parsed from a command line or a config, and two spellings of + the same number must key the same stream while genuinely different effect + sizes must not collide. + """ + if rep < 0: + raise ValueError(f"rep must be non-negative ({SETUP_REP} is the setup draw), got {rep}") + key = f"{target}\x00{rep}\x00{effect_size:.17g}" + return np.random.SeedSequence([int(seed), _text_entropy(key)]) + + +def rng_for(seed: int, target: str, rep: int, effect_size: float) -> np.random.Generator: + """`derive_seed`, as a generator ready to draw from.""" + return np.random.default_rng(derive_seed(seed, target, rep, effect_size)) diff --git a/src/watteg/simulate.py b/src/watteg/simulate.py new file mode 100644 index 0000000..11313f2 --- /dev/null +++ b/src/watteg/simulate.py @@ -0,0 +1,69 @@ +"""Drawing the counts. + +One function, and the only place randomness touches expression. Everything it +needs has been decided elsewhere: the baseline by `baseline.py`, the per-cell +multiplier by `perturbation.py`, the dispersion by `gene_model.py`. + + mu[i, j] = baseline[i, j] * effect_size[i, j] + count[i, j] ~ NegBinomial(mean = mu[i, j], size = theta[i]) + +`size` is theta, the same parameterisation sceptre's own model uses, so the +counts this draws are counts from the model the test assumes -- which is the +point of taking the baseline from that model too (see `baseline.py`). +""" + +from __future__ import annotations + +import numpy as np + +# int16 holds counts to 32,767. The largest single count in a real screen +# measured here is 2,717, but a simulated draw has a tail, so the promotion is +# checked rather than assumed. +_COUNT_DTYPE = np.int16 + + +def draw_counts( + baseline: np.ndarray, + effect_size: np.ndarray, + theta: np.ndarray, + rng: np.random.Generator, + *, + dtype: np.dtype | None = _COUNT_DTYPE, +) -> np.ndarray: + """Simulated counts, `(n_genes, n_cells)`. + + `theta` is the NB size, one per gene, broadcast down the rows. + + numpy parameterises the negative binomial as `(n, p)` with mean + `n(1-p)/p`, so `n = theta` and `p = theta / (theta + mu)` give mean `mu` + and variance `mu + mu^2/theta` -- R's `rnbinom(mu=, size=)`. + + Returned as `int16` by default. These are counts; holding them as float64 + costs four times the memory for no information, and the simulation's whole + shape depends on how many replicates fit in one process. + """ + baseline = np.asarray(baseline, dtype=float) + if baseline.shape != effect_size.shape: + raise ValueError(f"baseline is {baseline.shape} but effect_size is {effect_size.shape}") + theta = np.asarray(theta, dtype=float) + if theta.shape != (baseline.shape[0],): + raise ValueError(f"theta is {theta.shape}, expected ({baseline.shape[0]},)") + if np.any(theta <= 0) or not np.all(np.isfinite(theta)): + raise ValueError("theta must be finite and positive") + + mu = baseline * effect_size + size = theta[:, None] + # p == 1 exactly where mu == 0, which numpy accepts and which draws 0 + # every time -- the right answer for a gene knocked all the way down. + p = size / (size + mu) + counts = rng.negative_binomial(size, p) + + if dtype is None: + return counts + info = np.iinfo(dtype) + if counts.max(initial=0) > info.max: + raise OverflowError( + f"a simulated count exceeded {dtype.__name__}'s range ({info.max}); pass " + "dtype=None to keep the draw at full width" + ) + return counts.astype(dtype) diff --git a/tests/test_simulation_core.py b/tests/test_simulation_core.py new file mode 100644 index 0000000..26d44ce --- /dev/null +++ b/tests/test_simulation_core.py @@ -0,0 +1,248 @@ +"""The simulation core: seeding, guide assignment, effect sizes, the draw. + +These are the pieces a benchmark would otherwise measure without anyone having +checked they are right, so they are tested before anything times them. +""" + +from __future__ import annotations + +import numpy as np +import pytest +from scipy import sparse + +from watteg.perturbation import effect_size_matrix, guide_assignment, target_cells +from watteg.seeds import derive_seed, rng_for +from watteg.simulate import draw_counts + + +def tiny_screen(): + """Six cells, four guides, two targets. + + `gA1` and `gShared` hit `elemA`; `gShared` also hits `elemB`, which is the + many-to-many case a real screen produces from overlapping elements. Cell 5 + carries nothing. + """ + grna_ids = ["gA1", "gA2", "gShared", "gOther"] + membership = { + "gA1": [0, 1], + "gA2": [1, 2], + "gShared": [2, 3], + "gOther": [3, 4], + } + rows = np.concatenate([[i] * len(membership[g]) for i, g in enumerate(grna_ids)]) + cols = np.concatenate([membership[g] for g in grna_ids]) + grna = sparse.csr_matrix((np.ones(rows.size), (rows, cols)), shape=(4, 6)) + # elemA is the union of gA1, gA2 and gShared; elemB of gShared and gOther. + cre = sparse.csr_matrix(np.array([[1, 1, 1, 1, 0, 0], [0, 0, 1, 1, 1, 0]], dtype=float)) + return grna_ids, grna, ["elemA", "elemB"], cre + + +# --- seeding ----------------------------------------------------------------------------- + + +def test_the_same_key_gives_the_same_stream_and_a_different_key_does_not(): + a = rng_for(1, "elemA", 3, 0.15).normal(size=5) + assert np.array_equal(a, rng_for(1, "elemA", 3, 0.15).normal(size=5)) + for changed in ( + rng_for(2, "elemA", 3, 0.15), + rng_for(1, "elemB", 3, 0.15), + rng_for(1, "elemA", 4, 0.15), + rng_for(1, "elemA", 3, 0.20), + ): + assert not np.array_equal(a, changed.normal(size=5)) + + +def test_effect_size_keys_on_its_value_not_its_spelling(): + """It arrives parsed from a command line, so 0.15 and 0.150 are one effect + size and must not seed two different streams.""" + assert derive_seed(1, "t", 1, 0.15).entropy == derive_seed(1, "t", 1, 0.150).entropy + assert derive_seed(1, "t", 1, 0.15).entropy != derive_seed(1, "t", 1, 0.16).entropy + + +def test_the_setup_draw_cannot_collide_with_a_replicate(): + assert derive_seed(1, "t", 0, 0.15).entropy != derive_seed(1, "t", 1, 0.15).entropy + with pytest.raises(ValueError, match="non-negative"): + derive_seed(1, "t", -1, 0.15) + + +# --- guide assignment -------------------------------------------------------------------- + + +def test_a_cell_is_assigned_a_guide_it_actually_carries(): + grna_ids, grna, target_ids, cre = tiny_screen() + perturbed = target_cells(cre, target_ids, "elemA") + assert perturbed.tolist() == [True, True, True, True, False, False] + + a = guide_assignment( + grna.tocsc(), grna_ids, ["gA1", "gA2", "gShared"], perturbed, np.random.default_rng(0) + ) + assert a.n_target_guides == 3 + assert a.n_other_guides == 1 + + membership = {0: {0}, 1: {0, 1}, 2: {1, 2}, 3: {2, 3}, 4: {3}, 5: set()} + for cell in range(6): + if a.status[cell] == 0: + continue + if a.status[cell] <= 3: # one of the target's guides, numbered 1..3 + row = grna_ids.index(["gA1", "gA2", "gShared"][a.status[cell] - 1]) + else: # gOther, the only guide outside the target + row = grna_ids.index("gOther") + assert row in membership[cell], f"cell {cell} was given a guide it does not carry" + + +def test_a_cell_with_no_guide_gets_the_no_effect_row(): + grna_ids, grna, target_ids, cre = tiny_screen() + perturbed = target_cells(cre, target_ids, "elemA") + a = guide_assignment( + grna.tocsc(), grna_ids, ["gA1", "gA2", "gShared"], perturbed, np.random.default_rng(0) + ) + assert a.status[5] == 0 # cell 5 carries nothing + + +def test_control_cells_are_given_guides_from_outside_the_target(): + """The reason this matters: a control cell's guide draws an effect size + around 1 with the same spread. Sent to the no-effect row instead, the + control arm keeps its mean and loses its variance.""" + grna_ids, grna, target_ids, cre = tiny_screen() + perturbed = target_cells(cre, target_ids, "elemA") + a = guide_assignment( + grna.tocsc(), grna_ids, ["gA1", "gA2", "gShared"], perturbed, np.random.default_rng(0) + ) + assert a.status[4] > a.n_target_guides # cell 4 carries gOther only + assert a.status[:4].max() <= a.n_target_guides # perturbed cells never index past the block + + +def test_a_shared_guide_serves_both_of_its_targets(): + """gShared belongs to elemA and elemB. Read from a map that keeps one + target per guide, elemB would lose it.""" + grna_ids, grna, target_ids, cre = tiny_screen() + for target, guides in (("elemA", ["gA1", "gA2", "gShared"]), ("elemB", ["gShared", "gOther"])): + perturbed = target_cells(cre, target_ids, target) + a = guide_assignment(grna.tocsc(), grna_ids, guides, perturbed, np.random.default_rng(0)) + assert a.n_target_guides == len(guides) + + +def test_the_assignment_is_reproducible_and_depends_on_the_stream(): + grna_ids, grna, target_ids, cre = tiny_screen() + perturbed = target_cells(cre, target_ids, "elemA") + + def call(seed): + return guide_assignment( + grna.tocsc(), + grna_ids, + ["gA1", "gA2", "gShared"], + perturbed, + np.random.default_rng(seed), + ).status + + assert np.array_equal(call(0), call(0)) + + +def test_an_unknown_target_and_a_target_with_no_usable_guide_are_refused(): + grna_ids, grna, target_ids, cre = tiny_screen() + with pytest.raises(KeyError, match="not a row of cre_perts"): + target_cells(cre, target_ids, "nope") + perturbed = target_cells(cre, target_ids, "elemA") + with pytest.raises(KeyError, match="appear in grna_perts"): + guide_assignment(grna.tocsc(), grna_ids, ["absent"], perturbed, np.random.default_rng(0)) + + +# --- effect sizes ------------------------------------------------------------------------ + + +def test_each_arm_is_centred_on_what_it_should_average_to(): + """The property the centring step exists for: clamping negatives biases the + mean up, so "a 15% knockdown" would otherwise be slightly weaker.""" + rng = np.random.default_rng(3) + n_cells = 400 + status = rng.integers(0, 6, size=n_cells) + is_pert = np.zeros(n_cells, dtype=bool) + is_pert[: n_cells // 2] = True + status[is_pert] = rng.integers(0, 3, size=is_pert.sum()) # 0..2: no guide or target guides + status[~is_pert] = rng.integers(3, 6, size=(~is_pert).sum()) + a = type( + "A", + (), + {"status": status, "is_perturbed": is_pert, "n_target_guides": 2, "n_other_guides": 3}, + )() + + wanted = np.array([0.85, 0.5, 0.95]) + matrix = effect_size_matrix(a, wanted, guide_sd=0.13, rng=rng) + assert matrix.shape == (3, n_cells) + np.testing.assert_allclose(matrix[:, is_pert].mean(axis=1), wanted) + np.testing.assert_allclose(matrix[:, ~is_pert].mean(axis=1), 1.0) + + +def test_effect_sizes_never_go_negative(): + rng = np.random.default_rng(4) + n_cells = 200 + is_pert = np.zeros(n_cells, dtype=bool) + is_pert[:100] = True + status = np.where(is_pert, 1, 2) + a = type( + "A", + (), + {"status": status, "is_perturbed": is_pert, "n_target_guides": 1, "n_other_guides": 1}, + )() + # A near-zero target with a wide spread is where clamping actually bites. + matrix = effect_size_matrix(a, np.array([0.02]), guide_sd=0.5, rng=rng) + assert (matrix >= 0).all() + + +# --- the draw ---------------------------------------------------------------------------- + + +def test_counts_follow_the_negative_binomial_they_were_asked_for(): + rng = np.random.default_rng(5) + n_cells = 60_000 + baseline = np.full((2, n_cells), 4.0) + theta = np.array([2.0, 50.0]) + counts = draw_counts(baseline, np.ones((2, n_cells)), theta, rng, dtype=None) + + np.testing.assert_allclose(counts.mean(axis=1), 4.0, rtol=0.03) + # Var = mu + mu^2/theta, which is what distinguishes theta from noise. + np.testing.assert_allclose(counts.var(axis=1), 4.0 + 16.0 / theta, rtol=0.06) + + +def test_the_effect_size_scales_the_mean(): + rng = np.random.default_rng(6) + n_cells = 40_000 + counts = draw_counts( + np.full((1, n_cells), 10.0), + np.full((1, n_cells), 0.5), + np.array([20.0]), + rng, + dtype=None, + ) + np.testing.assert_allclose(counts.mean(), 5.0, rtol=0.03) + + +def test_a_fully_knocked_down_gene_draws_zeros_rather_than_failing(): + counts = draw_counts( + np.ones((1, 50)), np.zeros((1, 50)), np.array([3.0]), np.random.default_rng(7) + ) + assert (counts == 0).all() + + +def test_counts_come_back_as_int16_and_overflow_is_refused_not_wrapped(): + counts = draw_counts( + np.full((1, 20), 3.0), np.ones((1, 20)), np.array([5.0]), np.random.default_rng(8) + ) + assert counts.dtype == np.int16 + with pytest.raises(OverflowError, match="int16"): + draw_counts( + np.full((1, 200), 1e5), + np.ones((1, 200)), + np.array([1e6]), + np.random.default_rng(9), + ) + + +def test_the_draw_rejects_inputs_that_do_not_line_up(): + rng = np.random.default_rng(10) + with pytest.raises(ValueError, match="effect_size"): + draw_counts(np.ones((2, 5)), np.ones((2, 6)), np.array([1.0, 1.0]), rng) + with pytest.raises(ValueError, match="theta is"): + draw_counts(np.ones((2, 5)), np.ones((2, 5)), np.array([1.0]), rng) + with pytest.raises(ValueError, match="finite and positive"): + draw_counts(np.ones((2, 5)), np.ones((2, 5)), np.array([1.0, 0.0]), rng) From ca8ad94a1c701d23c7d7634495829f3ea03120be Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Mon, 21 Sep 2026 15:14:53 -0400 Subject: [PATCH 23/83] ADD the phase-3 benchmark: two engine shapes, four cost terms Times one pysceptre call per (target, replicate) against replicates stacked as pseudo-genes, and breaks the cost into drawing counts, the per-gene Poisson fits, the per-target binomial fit and CRT draws, and the per-pair test. Terms 2 and 3 are measured by calling fit_all_genes and fit_all_targets standalone rather than by instrumenting pysceptre, so they attribute the work rather than decomposing one call's clock -- the engine overlaps a chunk's target fit with the previous chunk's gene tests. The stage pass runs outside the windows that time the shapes; inside one, it was being charged to that shape and inflated the first chunk size measurably. Co-Authored-By: Claude Opus 5 (1M context) --- workflow/benchmark_simulation.py | 349 +++++++++++++++++++++++++++++++ 1 file changed, 349 insertions(+) create mode 100755 workflow/benchmark_simulation.py diff --git a/workflow/benchmark_simulation.py b/workflow/benchmark_simulation.py new file mode 100755 index 0000000..c5fe8c1 --- /dev/null +++ b/workflow/benchmark_simulation.py @@ -0,0 +1,349 @@ +#!/usr/bin/env python3 +"""Phase 3's gate: how much does the Python simulation actually cost? + + workflow/benchmark_simulation.py [--reps 100] [--chunk-sizes 10,25,50,100] + +Two ways to drive pysceptre, timed against each other: + +**per-rep** -- one `run_discovery_ntcells_complement` call per (target, +replicate), which is R's structure with a faster engine inside it. Every call +rebuilds the per-run state, including `x_outer_flat(covariate_matrix)`, which at +567,690 cells and p = 11 is not small. + +**stacked** -- one call per (target, replicate chunk), with the replicates +stacked as pseudo-genes `gene@rep` against pseudo-targets `target@es@rep`. The +pseudo-target keys are not cosmetic: pysceptre seeds each target's resampling +stream from its *name*, so one shared key would hand every replicate the same +synthetic treated sets and correlate replicates the Wilson interval assumes are +independent. The effect size is in the key for the same reason -- one `seed` +across a sweep would otherwise share draws between effect sizes. + +**Reported as four terms, not one wall clock**, because which one dominates is +what any later optimisation depends on: + +1. drawing the counts; +2. the per-(gene, replicate) Poisson fit -- unavoidable, and the price of the + faithful null model (`docs/pysceptre-backend.md` section 2.1); +3. the per-target binomial fit and CRT draws; +4. the redundancy in (3) from refitting identical cell sets per replicate, + measured as (3) minus the same work for a single target. + +Terms 2 and 3 are measured by calling `fit_all_genes` and `fit_all_targets` +**standalone**, on the same inputs the engine would hand them, rather than by +instrumenting pysceptre. So they are an attribution of where the work is, not a +decomposition of one call's clock: the engine overlaps a chunk's target fit +with the previous chunk's gene tests, so the parts can sum to more than the +whole. The remainder, total minus (2) minus (3), is the per-pair test. + +The comparison point is R's measured cost model, `1.140s + 0.5561s x pairs` per +(target, simulation), and 634 CPU-h per effect size on day0. +""" + +from __future__ import annotations + +import argparse +import resource +import sys +import time +from contextlib import contextmanager +from pathlib import Path + +import numpy as np +import pandas as pd + +from watteg.baseline import baseline_expression +from watteg.perturbation import effect_size_matrix, guide_assignment, target_cells +from watteg.seeds import SETUP_REP, rng_for +from watteg.sim_input import read_sim_input +from watteg.simulate import draw_counts + +# day0's measured R cost model, for the only comparison that matters. +R_INTERCEPT, R_PER_PAIR = 1.140, 0.5561 +R_CPU_HOURS = 634 +R_TARGETS, R_PAIRS = 3026, 34886 + +TIMINGS: dict[str, float] = {} + + +@contextmanager +def timed(term: str): + start = time.perf_counter() + try: + yield + finally: + TIMINGS[term] = TIMINGS.get(term, 0.0) + time.perf_counter() - start + + +def peak_rss_gb() -> float: + """macOS reports bytes, Linux kilobytes.""" + peak = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss + return peak / 1024**3 if sys.platform == "darwin" else peak / 1024**2 + + +def simulate_replicates(sim, target, genes, target_guides, effect_size, reps, seed, grna_csc): + """Counts for one target across `reps` replicates, stacked as rows. + + The guide assignment and the baseline are drawn once per target, not per + replicate: neither depends on the draw. That hoisting is in the R + implementation too, and it is why the setup seed is `rep = 0`. + """ + is_perturbed = target_cells(sim.cre_perts, sim.target_ids, target) + setup_rng = rng_for(seed, target, SETUP_REP, effect_size) + assignment = guide_assignment(grna_csc, sim.grna_ids, target_guides, is_perturbed, setup_rng) + + rows = [sim.genes.index(g) for g in genes] + baseline = baseline_expression( + "fitted", + fitted_coefs=sim.fitted_coefs[rows], + covariate_matrix=sim.covariate_matrix, + ) + theta = 1.0 / sim.row_data["dispersion"].to_numpy()[rows] + wanted = np.full(len(genes), 1.0 - effect_size) + + blocks = [] + for rep in range(1, reps + 1): + rng = rng_for(seed, target, rep, effect_size) + with timed("1. draw counts"): + es = effect_size_matrix(assignment, wanted, guide_sd=0.13, rng=rng) + blocks.append(draw_counts(baseline, es, theta, rng)) + return np.vstack(blocks), is_perturbed + + +def measure_terms(sim, target, genes, counts, is_perturbed, params, seed, reps, effect_size): + """Where the time goes, by running each stage on its own. + + The redundancy term is the honest way to price the `target@es@rep` keys: + every pseudo-target holds the same cell set, so all but one of their + binomial fits is arithmetic the engine has already done. Fitting one target + and `reps` of them measures exactly that difference. + """ + from pysceptre.pipeline.discovery import fit_all_genes, fit_all_targets + + cells = np.flatnonzero(is_perturbed) + gene_ids = [f"{g}@{r}" for r in range(1, reps + 1) for g in genes] + + with timed("2. per-gene Poisson fits"): + fit_all_genes( + counts, + gene_ids, + sim.covariate_matrix, + gene_rows=list(range(len(gene_ids))), + batch_width=1, + ) + draws = dict(B1=params["B1"], B2=params["B2"], B3=params["B3"], seed=seed) + with timed("3. target fits + CRT draws"): + fit_all_targets( + {f"{target}@{effect_size:g}@{r}": cells for r in range(1, reps + 1)}, + sim.covariate_matrix, + **draws, + ) + with timed("4. of which, one target's share"): + fit_all_targets({target: cells}, sim.covariate_matrix, **draws) + + +def call_engine(counts, gene_ids, covariates, target_cells_map, pairs, params, seed): + from pysceptre.pipeline.discovery import run_discovery_ntcells_complement + + return run_discovery_ntcells_complement( + counts, + gene_ids, + covariates, + target_cells_map, + pairs, + B1=params["B1"], + B2=params["B2"], + B3=params["B3"], + side_code=params["side_code"], + seed=seed, + ) + + +def run_per_rep(sim, target, genes, counts, is_perturbed, params, seed, reps): + """R's shape: one engine call per (target, replicate).""" + cells = np.flatnonzero(is_perturbed) + n_genes = len(genes) + out = [] + for rep in range(reps): + block = counts[rep * n_genes : (rep + 1) * n_genes] + pairs = pd.DataFrame({"response_id": genes, "grna_target": target}) + out.append( + call_engine( + block, genes, sim.covariate_matrix, {target: cells}, pairs, params, seed + rep + ) + ) + return pd.concat(out) + + +def run_stacked(sim, target, genes, counts, is_perturbed, params, seed, reps, effect_size): + """One engine call per (target, replicate chunk).""" + cells = np.flatnonzero(is_perturbed) + gene_ids = [f"{g}@{rep}" for rep in range(1, reps + 1) for g in genes] + target_map = {f"{target}@{effect_size:g}@{rep}": cells for rep in range(1, reps + 1)} + pairs = pd.DataFrame( + { + "response_id": gene_ids, + "grna_target": [ + f"{target}@{effect_size:g}@{rep}" for rep in range(1, reps + 1) for _ in genes + ], + } + ) + return call_engine(counts, gene_ids, sim.covariate_matrix, target_map, pairs, params, seed) + + +def main() -> int: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("prepared", type=Path, help="output of watteg-prepare-sim-input") + parser.add_argument("--reps", type=int, default=100) + parser.add_argument("--chunk-sizes", type=str, default="10,25,50,100") + parser.add_argument("--effect-size", type=float, default=0.15) + parser.add_argument("--seed", type=int, default=1) + parser.add_argument( + "--skip-per-rep", action="store_true", help="skip the per-rep shape, which is the slow one" + ) + args = parser.parse_args() + + sim = read_sim_input(args.prepared / "sim_input.h5") + pairs = pd.read_csv(args.prepared / "pairs.tsv", sep="\t") + design = pd.read_csv(args.prepared / "grna_targets.tsv", sep="\t") + mode = dict( + line.split("\t", 1) + for line in (args.prepared / "analysis_mode.tsv").read_text().splitlines() + if "\t" in line + ) + params = { + "B1": int(mode.get("B1", 499)), + "B2": int(mode.get("B2", 4999)), + "B3": int(mode.get("B3", 0)), + "side_code": {"left": -1, "both": 0, "right": 1}[mode.get("side", "both")], + } + print( + f"{sim.describe()}\n B1/B2/B3 = {params['B1']}/{params['B2']}/{params['B3']}, " + f"side_code = {params['side_code']}" + ) + + # One small, one median, one large target, so the per-pair term is visible. + by_target = pairs.groupby("grna_target")["response_id"].apply(list) + sizes = by_target.apply(len).sort_values() + chosen = [sizes.index[0], sizes.index[len(sizes) // 2], sizes.index[-1]] + print(" targets: " + ", ".join(f"{t} ({len(by_target[t])} pairs)" for t in chosen)) + + guides_of = design.groupby("grna_target")["grna_id"].apply(list) + grna_csc = sim.grna_perts.tocsc() + chunk_sizes = [int(c) for c in args.chunk_sizes.split(",")] + + print( + f"\n{'target':<28} {'pairs':>5} {'shape':<10} {'chunk':>6} {'s/(tgt,rep)':>12} " + f"{'peak GB':>8}" + ) + results = [] + for target in chosen: + genes = by_target[target] + counts, is_perturbed = simulate_replicates( + sim, + target, + genes, + guides_of[target], + args.effect_size, + args.reps, + args.seed, + grna_csc, + ) + n_trt = int(is_perturbed.sum()) + + # Its own pass: run inside a shape's timing window it would be counted + # as that shape's cost, which is what inflated the first chunk size on + # the first attempt. + stage_reps = min(chunk_sizes[0], args.reps) + measure_terms( + sim, + target, + genes, + counts[: stage_reps * len(genes)], + is_perturbed, + params, + args.seed, + stage_reps, + args.effect_size, + ) + + for shape in ("per-rep", "stacked"): + if shape == "per-rep" and args.skip_per_rep: + continue + for chunk in chunk_sizes if shape == "stacked" else [1]: + if chunk > args.reps: + continue + start = time.perf_counter() + done = 0 + while done < args.reps: + take = min(chunk, args.reps - done) + block = counts[done * len(genes) : (done + take) * len(genes)] + with timed("engine"): + if shape == "per-rep": + run_per_rep( + sim, target, genes, block, is_perturbed, params, args.seed, take + ) + else: + run_stacked( + sim, + target, + genes, + block, + is_perturbed, + params, + args.seed, + take, + args.effect_size, + ) + done += take + per_unit = (time.perf_counter() - start) / args.reps + print( + f"{target[:28]:<28} {len(genes):>5} {shape:<10} {chunk:>6} " + f"{per_unit:>12.3f} {peak_rss_gb():>8.2f}" + ) + results.append((target, len(genes), shape, chunk, per_unit, n_trt)) + + print("\n--- where the time goes, measured stage by stage ---") + engine_total = TIMINGS.pop("engine", 0.0) + for term, seconds in sorted(TIMINGS.items()): + print(f" {term:<32} {seconds:8.1f}s") + genes_t = TIMINGS.get("2. per-gene Poisson fits", 0.0) + targets_t = TIMINGS.get("3. target fits + CRT draws", 0.0) + one_t = TIMINGS.get("4. of which, one target's share", 0.0) + print(f" {'5. per-pair test (remainder)':<32} {engine_total - genes_t - targets_t:8.1f}s") + print(f" {'(engine calls, total)':<32} {engine_total:8.1f}s") + if targets_t: + print( + f"\n The redundancy the target@es@rep keys cost: {targets_t - one_t:.1f}s of " + f"{targets_t:.1f}s in term 3 ({100 * (targets_t - one_t) / targets_t:.0f}%), which is " + f"{100 * (targets_t - one_t) / max(engine_total, 1e-9):.1f}% of engine time. That is " + "the number section 5.3 of the plan would have to beat to be worth building." + ) + + print("\n--- against R ---") + frame = pd.DataFrame(results, columns=["target", "pairs", "shape", "chunk", "sec", "n_trt"]) + for shape, group in frame.groupby("shape"): + best = group.loc[group["sec"].idxmin()] + # Fit the same shape of model R's was reported in, if there is more than one size. + fit = np.polyfit(group["pairs"], group["sec"], 1) if group["pairs"].nunique() > 1 else None + model = ( + f"{fit[1]:.3f}s + {fit[0]:.4f}s x pairs" + if fit is not None + else f"{best['sec']:.3f}s at {int(best['pairs'])} pairs" + ) + predicted = ( + (fit[1] * R_TARGETS + fit[0] * R_PAIRS) * 100 / 3600 + if fit is not None + else float("nan") + ) + print(f" {shape:<10} {model}") + if fit is not None: + print( + f"{'':<12} -> {predicted:,.0f} CPU-h per effect size, against R's {R_CPU_HOURS} " + f"({R_CPU_HOURS / max(predicted, 1e-9):.1f}x)" + ) + print(f" R {R_INTERCEPT}s + {R_PER_PAIR}s x pairs -> {R_CPU_HOURS} CPU-h") + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) From 7f798770d6580cf109b5484a65441d908c0df78b Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Mon, 21 Sep 2026 15:28:37 -0400 Subject: [PATCH 24/83] Run the simulation in Python, one call per (target, replicate) The simple shape, chosen on measurement rather than assumption. Stacking replicates as pseudo-genes was built and timed first, because the plan assumed it was the architectural win; it is not. The per-target binomial fit and CRT draws are about 1.5% of the work and the duplicate fits stacking would remove about 0.8%, while the per-pair test is 85% and no call shape touches it. So the pseudo-target keys, the replicate-chunk memory knob and the question of whether replicates share a resampling stream all go away. Each call is its own target, so each replicate draws its own. pysceptre's defaults are left alone -- target_chunk_size and chunk_memory_gb are tuned. What is passed is not a default: B1/B2/B3 and the side are the screen's own analysis parameters, read from analysis_mode.tsv, and the inner entry point is used rather than the public one, which would size B2/B3 from the pair count. Output columns match the R implementation's, so the two can be compared without a conversion step, and a permutations or low-MOI screen is refused by name rather than silently analysed on the wrong path. Co-Authored-By: Claude Opus 5 (1M context) --- pyproject.toml | 1 + src/watteg/cli/run_power_simulation.py | 178 +++++++++++++++++++++++++ src/watteg/engine.py | 152 +++++++++++++++++++++ 3 files changed, 331 insertions(+) create mode 100644 src/watteg/cli/run_power_simulation.py create mode 100644 src/watteg/engine.py diff --git a/pyproject.toml b/pyproject.toml index f0596ea..aa673fe 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -16,6 +16,7 @@ dependencies = [ [project.scripts] watteg-prepare-sim-input = "watteg.cli.prepare_sim_input:main" +watteg-run-power-simulation = "watteg.cli.run_power_simulation:main" [project.optional-dependencies] dev = ["pytest>=7", "ruff>=0.5"] diff --git a/src/watteg/cli/run_power_simulation.py b/src/watteg/cli/run_power_simulation.py new file mode 100644 index 0000000..aced232 --- /dev/null +++ b/src/watteg/cli/run_power_simulation.py @@ -0,0 +1,178 @@ +"""Run the power simulation for one split, one effect size, one chunk of replicates. + + watteg-run-power-simulation --prepared prepared/ --pairs split_01.tsv \ + --effect-size 0.15 --reps 20 --rep-offset 0 --seed 1 --out sim.tsv + +For each target in the split and each replicate, this simulates a count matrix +under the given effect size and asks the screen's own test whether it would have +called the association. The fraction of replicates in which it would is the +power, computed downstream. + +Output columns match `src/run_power_simulation.R`'s, so the two implementations' +results can be compared without a conversion step and so downstream readers do +not care which produced a file. +""" + +from __future__ import annotations + +import argparse +import sys +import time +from pathlib import Path + +import numpy as np +import pandas as pd + +from watteg.engine import DEFAULT_GUIDE_SD, AnalysisParams, simulate_target +from watteg.sim_input import read_sim_input + +# What a reader needs, and nothing more. At 100 replicates x 34,886 pairs x six +# effect sizes the columns nothing reads were 39% of a 3 GB output. +KEEP = [ + "grna_target", + "response_id", + "p_value", + "log_2_fold_change", + "rep", + "effect_size", + "num_pert_cells", + "pass_qc", + "n_nonzero_trt", + "n_nonzero_cntrl", +] + + +def main(argv: list[str] | None = None) -> int: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument( + "--prepared", type=Path, required=True, help="output directory of watteg-prepare-sim-input" + ) + parser.add_argument( + "--pairs", + type=Path, + required=True, + help="one split, with columns grna_target and response_id", + ) + parser.add_argument( + "--effect-size", + type=float, + required=True, + help="a fractional decrease: 0.15 is a 15%% knockdown. 0 is the null arm", + ) + parser.add_argument("--reps", type=int, required=True) + parser.add_argument( + "--rep-offset", + type=int, + default=0, + help="replicates already covered by earlier chunks, so `rep` stays unique across them", + ) + parser.add_argument( + "--seed", + type=int, + required=True, + help="required: results are stochastic and must be reproducible", + ) + parser.add_argument("--out", type=Path, required=True) + parser.add_argument( + "--guide-sd", + type=float, + default=DEFAULT_GUIDE_SD, + help="across-gRNA spread of the effect size [default %(default)s]", + ) + parser.add_argument( + "--n-jobs", type=int, default=8, help="workers for the per-gene tests [default %(default)s]" + ) + parser.add_argument( + "--expression-model", + choices=("fitted", "size_factor"), + default="fitted", + help="'fitted' draws from exp(X.beta), sceptre's own null model. " + "'size_factor' reproduces the pre-2026-09-21 behaviour and exists " + "only for comparison against sweeps produced with it", + ) + args = parser.parse_args(argv) + + # Zero is allowed on purpose: it is the null arm, and running it is the only + # way to measure this pipeline's type-I error from its own output. + if not 0.0 <= args.effect_size < 1.0: + raise SystemExit( + f"--effect-size must be a fractional decrease in [0, 1) (got {args.effect_size}); " + "0 is the null arm" + ) + if args.reps < 1: + raise SystemExit("--reps must be at least 1") + + started = time.perf_counter() + sim = read_sim_input(args.prepared / "sim_input.h5") + params = AnalysisParams.from_analysis_mode(args.prepared / "analysis_mode.tsv") + design = pd.read_csv(args.prepared / "grna_targets.tsv", sep="\t") + split = pd.read_csv(args.pairs, sep="\t") + for column in ("grna_target", "response_id"): + if column not in split.columns: + raise SystemExit(f"{args.pairs} has no {column!r} column") + + # Many-to-many: a guide inside two overlapping elements belongs to both, so + # this is read from the design table and never from a per-unit annotation. + guides_of = design.groupby("grna_target")["grna_id"].apply(list) + genes_of = split.groupby("grna_target")["response_id"].apply(list) + reps = range(args.rep_offset + 1, args.rep_offset + args.reps + 1) + + print( + f"{len(genes_of)} targets / {len(split)} pairs, replicates {reps.start}-{reps.stop - 1}, " + f"effect size {args.effect_size} (relative expression {1 - args.effect_size:g})\n" + f" baseline: {args.expression_model}, n_jobs {args.n_jobs}, " + f"B1/B2/B3 {params.B1}/{params.B2}/{params.B3}, side_code {params.side_code}" + ) + + grna_csc = sim.grna_perts.tocsc() + frames = [] + for target, genes in genes_of.items(): + if target not in guides_of: + raise SystemExit(f"no gRNA maps to target {target!r} in grna_targets.tsv") + at = time.perf_counter() + frames.append( + simulate_target( + sim, + target, + list(genes), + guides_of[target], + effect_size=args.effect_size, + reps=reps, + seed=args.seed, + params=params, + grna_csc=grna_csc, + guide_sd=args.guide_sd, + n_jobs=args.n_jobs, + expression_model=args.expression_model, + ) + ) + elapsed = time.perf_counter() - at + print( + f" {target}: {len(genes)} pairs in {elapsed:.1f}s " + f"({elapsed / args.reps:.2f}s/replicate)" + ) + + combined = pd.concat(frames, ignore_index=True) + for column in ("pass_qc", "n_nonzero_trt", "n_nonzero_cntrl"): + if column not in combined: + combined[column] = np.nan + info = args.prepared / "pairs_with_info.tsv" + if info.exists(): + # Real-data QC counts, constant across replicates, joined rather than + # recomputed per draw -- which is what the R implementation does too. + known = pd.read_csv(info, sep="\t") + cols = [c for c in ("n_nonzero_trt", "n_nonzero_cntrl", "pass_qc") if c in known] + combined = combined.drop(columns=cols).merge( + known[["grna_target", "response_id", *cols]], + on=["grna_target", "response_id"], + how="left", + ) + + args.out.parent.mkdir(parents=True, exist_ok=True) + combined[[c for c in KEEP if c in combined]].to_csv(args.out, sep="\t", index=False) + print(f"\nwrote {len(combined):,} rows to {args.out} in {time.perf_counter() - started:.1f}s") + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/src/watteg/engine.py b/src/watteg/engine.py new file mode 100644 index 0000000..41c9e91 --- /dev/null +++ b/src/watteg/engine.py @@ -0,0 +1,152 @@ +"""Running the screen's own test on simulated counts. + +One call to pysceptre per (target, replicate). The alternative -- stacking +replicates as pseudo-genes against pseudo-targets, to amortise per-call setup -- +was built and measured, and there is almost nothing to amortise: the per-target +binomial fit and CRT draws are about 1.5% of the work and the duplicate fits it +would remove about 0.8%. The per-pair test dominates, and no call shape touches +it. So this is the simple shape, and it stays simple: no pseudo-genes, no +pseudo-target keys, no replicate-chunk memory knob, and no question about +whether replicates share a resampling stream -- each call is its own target, so +each replicate draws its own. + +**pysceptre's defaults are left alone.** `target_chunk_size` and +`chunk_memory_gb` are tuned, and the latter is documented as the fastest and +leanest setting measured. What is passed is not a default: `B1`, `B2`, `B3` and +the side are properties of the screen, read off its own analysis parameters, and +`n_jobs` is the caller's. + +The inner entry point rather than `pipeline.api.run_discovery_analysis`, which +sizes `B2`/`B3` from `len(pairs)` the way R's `run_qc` does. The screen's own +budget is what the real analysis used and is what a simulation of it has to +reuse. +""" + +from __future__ import annotations + +from dataclasses import dataclass +from pathlib import Path + +import numpy as np +import pandas as pd + +from .baseline import baseline_expression +from .perturbation import effect_size_matrix, guide_assignment, target_cells +from .seeds import SETUP_REP, rng_for +from .simulate import draw_counts + +DEFAULT_GUIDE_SD = 0.13 + + +@dataclass(frozen=True) +class AnalysisParams: + """The screen's own test configuration, not this pipeline's choices.""" + + B1: int + B2: int + B3: int + side_code: int + + @classmethod + def from_analysis_mode(cls, path) -> AnalysisParams: + text = Path(path).read_text() + mode = dict(line.split("\t", 1) for line in text.splitlines() if "\t" in line) + if mode.get("resampling_mechanism") != "crt": + raise ValueError( + f"this screen used {mode.get('resampling_mechanism')!r}, and the Python path " + "covers the CRT only. Run it through the R implementation." + ) + if mode.get("moi") == "low": + raise ValueError( + "this screen is low MOI, which pysceptre does not cover. Run it through the R " + "implementation." + ) + return cls( + B1=int(mode["B1"]), + B2=int(mode["B2"]), + B3=int(mode["B3"]), + side_code={"left": -1, "both": 0, "right": 1}[mode["side"]], + ) + + +def simulate_target( + sim, + target: str, + genes: list[str], + target_guides: list[str], + *, + effect_size: float, + reps: range, + seed: int, + params: AnalysisParams, + grna_csc, + guide_sd: float = DEFAULT_GUIDE_SD, + n_jobs: int = 8, + expression_model: str = "fitted", +) -> pd.DataFrame: + """Simulate and test one target, returning one row per (pair, replicate). + + The guide assignment and the baseline are drawn once per target: which + guide a cell carries is a fact about the screen, and the baseline depends + on neither the effect size nor the draw. That is what `rep = 0` is for. + """ + from pysceptre.pipeline.discovery import run_discovery_ntcells_complement + + is_perturbed = target_cells(sim.cre_perts, sim.target_ids, target) + n_perturbed = int(is_perturbed.sum()) + if n_perturbed == 0: + raise ValueError(f"target {target!r} perturbs no cell") + + assignment = guide_assignment( + grna_csc, + sim.grna_ids, + target_guides, + is_perturbed, + rng_for(seed, target, SETUP_REP, effect_size), + ) + + rows = [sim.genes.index(g) for g in genes] + baseline = baseline_expression( + expression_model, + fitted_coefs=sim.fitted_coefs[rows], + covariate_matrix=sim.covariate_matrix, + mean=sim.row_data["mean"].to_numpy()[rows], + size_factors=sim.col_data["size_factors"].to_numpy(), + ) + theta = 1.0 / sim.row_data["dispersion"].to_numpy()[rows] + # The config gives a fractional decrease; the simulation multiplies by a + # relative expression level. + wanted = np.full(len(genes), 1.0 - effect_size) + + treated = np.flatnonzero(is_perturbed) + pairs = pd.DataFrame({"response_id": genes, "grna_target": target}) + + out = [] + for rep in reps: + rng = rng_for(seed, target, rep, effect_size) + counts = draw_counts( + baseline, effect_size_matrix(assignment, wanted, guide_sd, rng), theta, rng + ) + result = run_discovery_ntcells_complement( + counts, + genes, + sim.covariate_matrix, + {target: treated}, + pairs, + B1=params.B1, + B2=params.B2, + B3=params.B3, + side_code=params.side_code, + seed=int(rng.integers(0, 2**31 - 1)), + n_jobs=n_jobs, + ) + result = result.assign(rep=rep, effect_size=effect_size, num_pert_cells=n_perturbed) + out.append(result) + + frame = pd.concat(out, ignore_index=True) + # log2 of the fold change, because that is the column the power step reads + # and the one the R implementation writes. `fold_change < 1` and + # `log_2_fold_change < 0` are the same predicate; carrying the log keeps + # the two implementations' outputs comparable without a conversion step. + frame["log_2_fold_change"] = np.log2(frame["fold_change"]) + return frame From 7bbb8d3d5398c20641cffce84e1495d4a983dd47 Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Mon, 21 Sep 2026 15:36:31 -0400 Subject: [PATCH 25/83] Check the two implementations agree, against noise rather than a tolerance They cannot agree draw for draw -- different generators, so different counts and different CRT draws -- so the bar is the Monte Carlo noise between two correct runs. On 40 replicates of one target: per-pair power agrees for 8 of 8 pairs within 2 sd, with no systematic shift (+0.022 against a 2-se bound of 0.039), and every pair's fold change agrees within its own replicate spread. The fold-change check is the sharpest because it is ABSOLUTE. Both implementations should land on log2(1 - effect_size) whatever their draws, and both do: medians of -0.2325 and -0.2277 against log2(0.85) = -0.2345. A comparison between two implementations cannot catch both being wrong in the same direction; that one can. The first version of this script asked the fold changes to agree to 0.05 in log2, a number picked by eye, and reported a failure at 0.0883 -- against that pair's own 2-se bound of 0.1286, because its per-replicate spread is 0.24 while another pair's is 0.02. No fixed tolerance is right for both. Every bound here is now derived per pair from the data's own spread, which is the rule this repository applies to R's numbers and should apply to its own. Co-Authored-By: Claude Opus 5 (1M context) --- workflow/compare_simulation.py | 164 +++++++++++++++++++++++++++++++++ 1 file changed, 164 insertions(+) create mode 100755 workflow/compare_simulation.py diff --git a/workflow/compare_simulation.py b/workflow/compare_simulation.py new file mode 100755 index 0000000..19e191c --- /dev/null +++ b/workflow/compare_simulation.py @@ -0,0 +1,164 @@ +#!/usr/bin/env python3 +"""Do the two implementations give the same answer? + + workflow/compare_simulation.py --threshold-file discovery_threshold.txt + +**They cannot agree draw for draw and it would be wrong to want them to.** The +two use different random number generators, so a given replicate's counts and +its CRT draws differ. What has to agree is the quantity the pipeline reports: +per-pair power, the fraction of replicates in which the screen's own test would +have called the association. + +So this is a statistical comparison, and the bar is the Monte Carlo noise +between two correct runs rather than a tolerance. Two independent estimates of +the same binomial `p` over `n` replicates differ with standard deviation +`sqrt(2p(1-p)/n)`: at `p = 0.5` and `n = 40` that is 0.11, which is what +"agreement" can mean at this replicate count. Raising the replicate count is +the only thing that tightens it. + +Reported alongside power, because power at the extremes hides everything: the +per-pair median p-value on a log scale, which is sensitive where power is +saturated at 0 or 1, and the fold change, which is the sharpest signal +available because it is an absolute quantity -- both implementations should +land on `log2(1 - effect_size)`, whatever their draws. + +**Every bound here is derived from the data's own spread, per pair.** A first +version of this script asked the fold changes to agree to 0.05 in log2, a +number picked by eye, and one pair missed it at 0.0883 -- against its own +two-standard-error bound of 0.1286, because that gene's per-replicate spread +is 0.24 while another's is 0.02. A fixed tolerance is either too tight for the +noisy genes or too loose for the quiet ones, and there is no single value that +is neither. +""" + +from __future__ import annotations + +import argparse +import sys +from pathlib import Path + +import numpy as np +import pandas as pd + +FAILURES: list[str] = [] + + +def check(label: str, ok: bool, detail: str = "") -> None: + print(f"{'PASS' if ok else 'FAIL'} {label}{' ' + detail if detail else ''}") + if not ok: + FAILURES.append(label) + + +def power_per_pair(frame: pd.DataFrame, threshold: float) -> pd.DataFrame: + """The pipeline's definition: a replicate counts only if the test called it + AND the simulated perturbation reduced expression.""" + called = (frame["p_value"] < threshold) & (frame["log_2_fold_change"] < 0) + return ( + frame.assign(called=called) + .groupby(["grna_target", "response_id"]) + .agg( + power=("called", "mean"), + reps=("called", "size"), + median_log10_p=("p_value", lambda p: float(np.median(np.log10(np.maximum(p, 1e-300))))), + median_log2fc=("log_2_fold_change", "median"), + sd_log2fc=("log_2_fold_change", "std"), + ) + .reset_index() + ) + + +def main() -> int: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("python_tsv", type=Path) + parser.add_argument("r_tsv", type=Path) + parser.add_argument("--threshold-file", type=Path, required=True) + args = parser.parse_args() + + threshold = float(args.threshold_file.read_text().strip()) + py = pd.read_csv(args.python_tsv, sep="\t") + r = pd.read_csv(args.r_tsv, sep="\t") + print(f"threshold {threshold:.6g}") + print(f"python: {len(py):,} rows, {py['rep'].nunique()} replicates") + print(f"r: {len(r):,} rows, {r['rep'].nunique()} replicates\n") + + a, b = power_per_pair(py, threshold), power_per_pair(r, threshold) + key = ["grna_target", "response_id"] + check( + "the same pairs were simulated", + set(map(tuple, a[key].to_numpy())) == set(map(tuple, b[key].to_numpy())), + f"{len(a)} pairs", + ) + merged = a.merge(b, on=key, suffixes=("_py", "_r")) + n = int(merged["reps_py"].iloc[0]) + check("the same replicate count", (merged["reps_py"] == merged["reps_r"]).all(), f"{n}") + + d = merged["power_py"] - merged["power_r"] + # The noise floor for the comparison, per pair, from R's estimate. + p = merged["power_r"].clip(1e-9, 1 - 1e-9) + sd = np.sqrt(2 * p * (1 - p) / n) + within = np.abs(d) <= 2 * sd + 1e-12 + print(f"\npower, {len(merged)} pairs at {n} replicates") + print(f" mean difference {d.mean():+.4f} (expected ~0)") + print(f" mean |difference| {np.abs(d).mean():.4f}") + print(f" largest |difference| {np.abs(d).max():.4f}") + print(f" noise floor, 2 sd {2 * sd.mean():.4f} on average") + check( + "every pair within 2 sd of Monte Carlo noise", + bool(within.all()), + f"{int(within.sum())}/{len(within)}", + ) + # A mean shift is what a systematic difference looks like; individual pairs + # crossing is what noise looks like. + se_mean = float(np.sqrt(np.sum(sd**2)) / len(sd)) + check( + "no systematic shift in power", + abs(d.mean()) <= 2 * se_mean + 1e-12, + f"{d.mean():+.4f} against a 2-se bound of {2 * se_mean:.4f}", + ) + + print("\nfold change in log2 -- an absolute quantity, so it is the sharpest check") + fc = merged["median_log2fc_py"] - merged["median_log2fc_r"] + # The standard error of a median is about 1.2533 sd / sqrt(n) for a normal, + # and these two are independent, so their difference has the pooled one. + fc_se = 1.2533 * np.sqrt((merged["sd_log2fc_py"] ** 2 + merged["sd_log2fc_r"] ** 2) / n) + print(f" mean difference {fc.mean():+.4f}") + print(f" largest |difference| {np.abs(fc).max():.4f}") + print(f" its own 2-se bound {2 * fc_se[np.abs(fc).idxmax()]:.4f}") + check( + "every pair's fold change within 2 se of its own replicate spread", + bool((np.abs(fc) <= 2 * fc_se + 1e-12).all()), + f"{int((np.abs(fc) <= 2 * fc_se + 1e-12).sum())}/{len(fc)}", + ) + + # The absolute check: whatever the draws, a 15% knockdown should land on + # log2(0.85). This is what would catch both implementations being wrong in + # the same direction, which no comparison between them can. + effect_size = float(py["effect_size"].iloc[0]) + wanted = np.log2(1.0 - effect_size) + for label, column in (("python", "median_log2fc_py"), ("r", "median_log2fc_r")): + off = merged[column] - wanted + check( + f"{label} hits the requested knockdown (log2 {wanted:.4f})", + bool((np.abs(off) <= 2 * fc_se + 1e-12).all()), + f"median {merged[column].median():+.4f}, largest miss {np.abs(off).max():.4f}", + ) + + print("\nmedian p-value per pair, log10 -- sensitive where power saturates") + lp = merged["median_log10_p_py"] - merged["median_log10_p_r"] + print(f" mean difference {lp.mean():+.3f} decades") + print(f" largest |difference| {np.abs(lp).max():.3f} decades") + + saturated = int(((merged["power_r"] == 0) | (merged["power_r"] == 1)).sum()) + if saturated: + print( + f"\n NOTE: {saturated} of {len(merged)} pairs have power exactly 0 or 1 in R, " + "where power cannot distinguish the two implementations. Read the fold change and " + "the p-value rows for those." + ) + + print("\n" + ("ALL PASS" if not FAILURES else f"{len(FAILURES)} FAILED: {FAILURES}")) + return 1 if FAILURES else 0 + + +if __name__ == "__main__": + sys.exit(main()) From 6c8966fc7a284f2b0637a69efe6f6ab8b58a577b Mon Sep 17 00:00:00 2001 From: Eugenio Mattei Date: Mon, 21 Sep 2026 15:44:02 -0400 Subject: [PATCH 26/83] Port the remaining four pipeline steps split_pairs, consolidate_replicates, compute_power and summarize_power, which is everything between the simulation and the published tables. Checked against R on identical input, which makes these exact rather than statistical comparisons. compute_power agrees on every column to machine epsilon -- 6.7e-16 at worst, which is summation order. summarize_power agrees on every column and in the same order, with one exception that is not one: `dispersion` differs by 7.3e-12, because it is pysceptre's fitted theta against sceptre's cache, and that has its own gate at 1e-6. Two things the comparison caught: max_effect_size_tested was missing altogether, and the per-gene columns sat after the power columns rather than before them. The subtleties that have gone wrong before are tested for the specific failure. Wilson rather than the normal approximation, because power_ci_low is what establishes a negative and 0 of 100 must not read as [0, 0]. effect_label not trimming 0.2 into "2". The suffix rule on min_detectable_effect_size, which exists because taking the first effect size that clears lets Monte Carlo noise win in one direction only. And the interval on power inverting into the interval on effect size, where getting it backwards would report the conservative edge as the optimistic one. Two of those tests failed first and were wrong themselves: they asserted a Wilson bound was exactly 0, where center - halfwidth cancels to a denormal. R's does the same; its printed 0.000000 is rounding. Co-Authored-By: Claude Opus 5 (1M context) --- pyproject.toml | 4 + src/watteg/cli/compute_power.py | 87 ++++++++++++ src/watteg/cli/consolidate_replicates.py | 59 ++++++++ src/watteg/cli/split_pairs.py | 103 ++++++++++++++ src/watteg/cli/summarize_power.py | 158 ++++++++++++++++++++ src/watteg/power.py | 168 ++++++++++++++++++++++ tests/test_power.py | 174 +++++++++++++++++++++++ 7 files changed, 753 insertions(+) create mode 100644 src/watteg/cli/compute_power.py create mode 100644 src/watteg/cli/consolidate_replicates.py create mode 100644 src/watteg/cli/split_pairs.py create mode 100644 src/watteg/cli/summarize_power.py create mode 100644 src/watteg/power.py create mode 100644 tests/test_power.py diff --git a/pyproject.toml b/pyproject.toml index aa673fe..fe0fbd0 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -17,6 +17,10 @@ dependencies = [ [project.scripts] watteg-prepare-sim-input = "watteg.cli.prepare_sim_input:main" watteg-run-power-simulation = "watteg.cli.run_power_simulation:main" +watteg-split-pairs = "watteg.cli.split_pairs:main" +watteg-consolidate-replicates = "watteg.cli.consolidate_replicates:main" +watteg-compute-power = "watteg.cli.compute_power:main" +watteg-summarize-power = "watteg.cli.summarize_power:main" [project.optional-dependencies] dev = ["pytest>=7", "ruff>=0.5"] diff --git a/src/watteg/cli/compute_power.py b/src/watteg/cli/compute_power.py new file mode 100644 index 0000000..6a6ebab --- /dev/null +++ b/src/watteg/cli/compute_power.py @@ -0,0 +1,87 @@ +"""Turn per-replicate simulation results into a power estimate per pair. + + watteg-compute-power --simulations sims/ --threshold-file discovery_threshold.txt \ + --out power_es0.15.tsv + +`--simulations` takes files or directories; a directory expands to the TSV and +Parquet files inside it, sorted, so row order does not depend on the filesystem. +See `watteg/power.py` for what power means here and why the interval is Wilson's. +""" + +from __future__ import annotations + +import argparse +import sys +from pathlib import Path + +import pandas as pd + +from watteg.power import compute_power + +SUFFIXES = (".parquet", ".tsv", ".tsv.gz") + + +def expand(paths: list[Path]) -> list[Path]: + out: list[Path] = [] + for path in paths: + if path.is_dir(): + out.extend(sorted(p for p in path.iterdir() if p.name.endswith(SUFFIXES))) + else: + out.append(path) + return out + + +def read_one(path: Path) -> pd.DataFrame: + if path.suffix == ".parquet": + return pd.read_parquet(path) + return pd.read_csv(path, sep="\t") + + +def main(argv: list[str] | None = None) -> int: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--simulations", nargs="+", type=Path, required=True) + parser.add_argument("--threshold-file", type=Path) + parser.add_argument( + "--alpha", type=float, help="set the threshold explicitly instead of reading it" + ) + parser.add_argument("--conf-level", type=float, default=0.95) + parser.add_argument("--out", type=Path, required=True) + args = parser.parse_args(argv) + + if (args.threshold_file is None) == (args.alpha is None): + raise SystemExit("pass exactly one of --threshold-file and --alpha") + threshold = ( + args.alpha if args.alpha is not None else float(args.threshold_file.read_text().split()[0]) + ) + + paths = expand(args.simulations) + if not paths: + raise SystemExit(f"--simulations matched no files: {args.simulations}") + frame = pd.concat([read_one(p) for p in paths], ignore_index=True) + print( + f"read {len(frame):,} simulation rows from {len(paths)} file(s), threshold {threshold:.6g}" + ) + + power = compute_power(frame, threshold, conf_level=args.conf_level) + args.out.parent.mkdir(parents=True, exist_ok=True) + power.to_csv(args.out, sep="\t", index=False) + + width = power["power_ci_high"] - power["power_ci_low"] + print(f"wrote {len(power):,} pairs to {args.out}") + print(f" replicates per pair: {power['n_reps'].min()}-{power['n_reps'].max()}") + print( + f" mean power {power['power'].mean():.3f} | at 0: {(power['power'] == 0).sum():,} " + f"| at 1: {(power['power'] == 1).sum():,} " + f"| in (0.1, 0.9): {power['power'].between(0.1, 0.9, 'neither').sum():,}" + ) + print(f" median 95% CI width {width.median():.3f} (widest {width.max():.3f})") + if power["n_reps"].min() < 100: + print( + " note: below ~100 replicates a per-pair estimate is coarse; see " + "docs/choosing-num-replicates.md" + ) + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/src/watteg/cli/consolidate_replicates.py b/src/watteg/cli/consolidate_replicates.py new file mode 100644 index 0000000..1a65ece --- /dev/null +++ b/src/watteg/cli/consolidate_replicates.py @@ -0,0 +1,59 @@ +"""Concatenate a sweep's per-split simulation output into one Parquet file. + + watteg-consolidate-replicates --simulations sims/ --out per_replicate/es0.15.parquet + +The per-replicate output is the one thing power can be re-derived from, so it is +read repeatedly -- and as 1,000 gzipped TSVs that costs 30-90 s of parsing every +time. One Parquet file per effect size is seconds. + +The split files stay in the work directory, so a failed consolidation loses +nothing. +""" + +from __future__ import annotations + +import argparse +import sys +from pathlib import Path + +import pandas as pd + +from watteg.cli.compute_power import expand, read_one + +# Stored as categories: a sweep repeats a few thousand target and gene names +# across millions of rows, and dictionary encoding is what makes the file small. +CATEGORICAL = ("grna_target", "response_id") + + +def main(argv: list[str] | None = None) -> int: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--simulations", nargs="+", type=Path, required=True) + parser.add_argument("--out", type=Path, required=True) + parser.add_argument("--compression", default="zstd") + args = parser.parse_args(argv) + + paths = expand(args.simulations) + if not paths: + raise SystemExit(f"--simulations matched no files: {args.simulations}") + frame = pd.concat([read_one(p) for p in paths], ignore_index=True) + + if "effect_size" in frame.columns: + sizes = frame["effect_size"].unique() + if len(sizes) > 1: + raise SystemExit( + f"the input mixes effect sizes ({', '.join(map(str, sizes))}); consolidate " + "one effect size per file, or power would be averaged across knockdown levels" + ) + for column in CATEGORICAL: + if column in frame.columns: + frame[column] = frame[column].astype("category") + + args.out.parent.mkdir(parents=True, exist_ok=True) + frame.to_parquet(args.out, compression=args.compression, index=False) + size_mb = args.out.stat().st_size / 1e6 + print(f"wrote {len(frame):,} rows from {len(paths)} file(s) to {args.out} ({size_mb:.1f} MB)") + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/src/watteg/cli/split_pairs.py b/src/watteg/cli/split_pairs.py new file mode 100644 index 0000000..c97fdcc --- /dev/null +++ b/src/watteg/cli/split_pairs.py @@ -0,0 +1,103 @@ +"""Split the discovery pairs into per-task chunks, balanced by cost. + + watteg-split-pairs --pairs pairs.tsv --outdir splits/ --n-splits 480 + +A target's pairs cannot be separated: the simulation draws one count matrix per +(target, replicate) and tests every one of that target's genes against it, so +splitting a target would simulate it twice. Targets are therefore the unit, and +the task is bin packing them. + +**Weighted by pairs plus a per-target overhead**, because a task's cost is not +proportional to its pairs. The measured model is `intercept + slope x pairs`, +and on day0 the intercept is a real share of a small target's cost, so weighting +by pairs alone over-fills the splits that hold many small targets. + +Longest-processing-time-first: sort by weight descending, then put each target +in whichever split is currently lightest. Deterministic, and within a small +constant factor of optimal for this shape. +""" + +from __future__ import annotations + +import argparse +import sys +from pathlib import Path + +import pandas as pd + + +def assign_splits(weights: pd.Series, n_splits: int) -> pd.Series: + """Target -> split number (1-based), by longest-processing-time-first.""" + load = [0.0] * n_splits + assignment = {} + for target, weight in weights.sort_values(ascending=False).items(): + lightest = min(range(n_splits), key=lambda k: load[k]) + assignment[target] = lightest + 1 + load[lightest] += float(weight) + return pd.Series(assignment, name="split") + + +def main(argv: list[str] | None = None) -> int: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--pairs", type=Path, required=True) + parser.add_argument("--outdir", type=Path, required=True) + parser.add_argument("--n-splits", type=int, required=True) + parser.add_argument("--prefix", default="split_") + parser.add_argument( + "--target-overhead", + type=float, + default=2.0, + help="pairs-equivalent fixed cost of a target, from the measured cost model " + "[default %(default)s]", + ) + args = parser.parse_args(argv) + + if args.target_overhead < 0: + raise SystemExit("--target-overhead must not be negative") + + pairs = pd.read_csv(args.pairs, sep="\t") + for column in ("grna_target", "response_id"): + if column not in pairs.columns: + raise SystemExit(f"{args.pairs} has no {column!r} column") + if pairs.empty: + raise SystemExit(f"{args.pairs} contains no pairs") + + per_target = pairs.groupby("grna_target").size() + if args.n_splits > len(per_target): + raise SystemExit( + f"--n-splits ({args.n_splits}) exceeds the number of targets ({len(per_target)}). " + "Every split must hold at least one target; lower --n-splits." + ) + print( + f"{len(pairs):,} pairs over {len(per_target):,} targets " + f"({per_target.min()}-{per_target.max()} each, median {per_target.median():g})" + ) + + split_of = assign_splits(per_target + args.target_overhead, args.n_splits) + pairs = pairs.assign(split=pairs["grna_target"].map(split_of)) + + args.outdir.mkdir(parents=True, exist_ok=True) + # Zero-padded so lexicographic order matches numeric order, which keeps a + # workflow engine's channel ordering and a manual listing predictable. + width = max(2, len(str(args.n_splits))) + written = 0 + sizes = [] + for k in range(1, args.n_splits + 1): + chunk = pairs.loc[pairs["split"] == k, ["grna_target", "response_id"]] + chunk = chunk.sort_values(["grna_target", "response_id"]) + chunk.to_csv(args.outdir / f"{args.prefix}{k:0{width}d}.tsv", sep="\t", index=False) + written += len(chunk) + sizes.append(len(chunk)) + + if written != len(pairs): + raise SystemExit(f"wrote {written} pairs but read {len(pairs)}; the split lost rows") + print(f"wrote {args.n_splits} splits to {args.outdir}") + print( + f" pairs per split: {min(sizes)}-{max(sizes)}, " + f"imbalance (max/min) {max(sizes) / max(min(sizes), 1):.3f}" + ) + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/src/watteg/cli/summarize_power.py b/src/watteg/cli/summarize_power.py new file mode 100644 index 0000000..97ae41f --- /dev/null +++ b/src/watteg/cli/summarize_power.py @@ -0,0 +1,158 @@ +"""One row per pair, across every effect size in a sweep. + + watteg-summarize-power --power power_es0.05.tsv power_es0.15.tsv ... \ + --sim-input prepared/sim_input.h5 --out power_summary.tsv + +Each effect size contributes `power_at_effect_size_