Switch main to the Python (pysceptre) pipeline - #2
Merged
Merged
Conversation
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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
--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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
…rance 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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
…ronment Six processes instead of eight. FIT_NULL_MODELS and MERGE_NULL_MODELS are gone: they existed because R refitting a gene's null model inside every call cost 4.3x, and the Python path does that refit as a matter of course. With them go reps_per_null_chunk, test_max_null_reps and the divisibility check between them and the replicate count. The samplesheet takes a `dataset` column -- a .h5mu from pysceptre's export -- rather than `sceptre_object`, and the response_odm column is gone because the export handles odm-backed matrices itself. Nothing in the pipeline reads R, so pixi.toml holds no R: not r-base, but also not the pinned sceptre commit, its local patch, the pinned ondisc commit, the two source-install tasks, the check-api task that verified the pin against unexported S4 slots, or the HDF5/curl/openssl link-line dependencies ondisc's static build needed. pysceptre enters as a git dependency pinned to a commit, because a branch moves and a sweep has to be re-runnable against the engine it was produced with. Four things had to change that predate this work: conf/base.config has never parsed under the Nextflow 26.04 that pixi.lock already pinned. Its strict syntax rejects a closure defined inside a config closure (the asMem helper, now inlined -- MemoryUnit's toString round-trips, so one expression still accepts both the MemoryUnit the config declares and the String a params file sends), reads `(long) (x)` as a call to `long`, rejects top-level statements like workflow.onComplete, and does not see a `def` closure from inside a workflow body, so resolve() is a function now. Dynamic publishDir paths need the closure form or `meta` is not yet in scope. Verified: a stub run completes all eight tasks across all six processes. The POWER_SIMULATION resource numbers are still R's, and they are now labelled as such rather than left to look calibrated. A task uses task.cpus workers where R's used one, and the footprint is a different function of the cell and gene counts, so they err high until a real sweep's trace replaces them. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
How many splits a cluster wanted and how many replicates fit in a task are operational choices. Four replicates in one task against two tasks of two come out BYTE-identical, not merely similar: the draws are keyed on (seed, target, replicate, effect size) and nothing else, so a chunk boundary has nowhere to enter. The seed-level half of that contract was already tested and runs everywhere. This is the half that would catch it being broken downstream of the seed -- by a per-task generator, a shared one, or anything consuming draws in an order that depends on the chunking -- so it runs the real engine and is opt-in behind -m realdata. Its companion asserts a different seed DOES change the draws, because invariance that came from ignoring the seed would pass the first test and be worthless. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
README's installation, input and quickstart all described the R path. The environment is now `pixi install` and nothing else; the input is a .h5mu exported once per dataset, which is the only place R appears and is not part of the pipeline; and the quickstart drives the watteg-* console scripts. --all-cells is called out as required rather than listed as a flag, because an export without it silently gives different size factors for the cells that remain, and nothing downstream would say so. status.md keeps its history -- the measured numbers and correctness fixes are still true of the R implementation -- behind a note saying which of its conclusions the port has overtaken, including the null-model section it spends the most space settling. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Phases 1-6 are done. Stage D of the validation is done and is the only one that did not need a cluster. The plan now carries a section on what is left rather than leaving it to be inferred from unticked boxes: the resource numbers are still R's and are not a calibration until a trace replaces them, no container exists for the gcb profile, and Stage B at full scale is the thing that decides whether this replaces the R path. Also says plainly what the equivalence evidence is and is not. Everything measured says the two paths agree -- on one target, eight pairs, forty replicates. That is agreement where they have been compared, which is a different claim from agreement. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
From the second adversarial review, each finding upheld 3/3: - _pin_to_mean() carried the rounding of a long cumulative sum. Perturbed cells take one of a few guide values, so the error adds up instead of cancelling, and at ~5e4 perturbed cells it passed the 1e-12 check and raised on a valid pin. R and Python could even disagree on identical input, because numpy's mean is pairwise and R's rowMeans is not. One correction step on the linear piece takes the worst miss to 5.6e-16 up to 3e5 cells. The check's tolerance floor now grows with the cell count, as summation error does. New test at 5e4, 1e5 and 3e5 cells. The docstring says the no-clamp case is the plain shift "to a few ulps", which is what it is. - The R reference copy on this branch now matches r-implementation HEAD (6a09df0) file for file: lib/simulate.R, src/prepare_sim_input.R and tests/testthat/test-simulate.R take the same fixes. src/fit_null_models.R, which was already stale here and called draw_counts() without the baseline this branch's lib/ requires, is brought across too. 1,238 R assertions pass. 48 Python unit tests and the 2 realdata tests pass. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
…r Wave's limit watteg-prepare-sim-input (and make_test_data) import sceptre_io from pysceptre's scripts/, found relative to the installed package. pysceptre is a git dependency that pixi installs as a wheel, and the wheel does not ship scripts/. In the image, every run would therefore have died at PREPARE_SIM_INPUT with "No module named sceptre_io". It never showed, because no real cloud run of the Python pipeline has happened; locally it needs PYTHONPATH pointed at a pysceptre checkout. The image now clones pysceptre at the exact commit pinned in pixi.lock (read from the lock, so the two cannot drift), copies scripts/ to /opt/pysceptre-scripts and puts it on PYTHONPATH. The smoke test imports sceptre_io. The recipe is cut from 8.5 KB to 2.8 KB (3.7 KB base64). Wave rejects a containerfile over 10,000 bytes base64, which the old one exceeded (11.3 KB), the same limit the R image hit. The directives are unchanged apart from the pysceptre step; the rationale that stays is what a maintainer needs. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Points at the dataset.h5mu exported on 2026-09-25 from the same moi5 sceptre object the R samplesheet uses (38,606 genes x 134,790 cells, 2,974 targets, 33,066 QC-passing pairs), so the two implementations run on identical input for the Stage B comparison. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
…s R does Two things that would have stopped the Python pipeline on moi5 at the first real task, found by running it locally on the moi5 cis export before any cloud launch. - Both moi5 screens used sceptre's permutation test, and AnalysisParams refused anything but the CRT, so Python could not simulate moi5 at all. pysceptre implements permutations (run_discovery_ntcells_complement takes resampling_mechanism), so the engine now passes the screen's own mechanism through with its B1/B2/B3, and the simulation re-runs the test the screen ran. No change to pysceptre. Its docs call permutations exercised rather than validated against R, so Stage B (R vs Python on moi5 cis) is also that validation. On real moi5 cis targets it runs at ~0.1 s per pair per replicate on one CPU, peak RSS 1.6 GB. - prepare refused any gene whose theta sits at the estimator's bound. That stopped moi5 cis on one gene (ENSG00000186409, mean 0.0085 counts per cell, theta at 0.01) and would have dropped its pairs while R kept them. sceptre's cache holds the same clamped value, R simulates from it, and sceptre's test uses it in its null model. The gene is now kept at the bound with a warning. With it, Python prepare on moi5 cis gives 33,066 pairs and a discovery threshold byte-identical to R's. Also wraps a long comment in baseline.py that failed ruff. New test for the resampling-mechanism parsing; 49 unit tests pass. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
…(1 - es) Same change as r-implementation 1f44a42, in the Python engine. Each guide's knockdown is now Beta with mean es and sd c * es * (1 - es), c = 0.65, via draw_guide_effects(). It replaces the absolute N(1 - es, 0.13), whose spread did not vanish at es = 0. Around the pinned mean that inflated the false-call rate of highly expressed genes: 3,032 moi5 trans pairs above twice the nominal rate at the discovery threshold, the worst at 41x. The form and value come from a research workflow on 2026-09-25. Per-guide data from DC-TAP K562, DC-TAP WTC11 and moi5 show spread ~0 at es = 0, rising in proportion and levelling off. The form fits at chi-square 6.6 on 4 df, against 128 for the absolute one. Power moves by no more than 0.02 for any moi5 pair. Draws lie in (0, 1), so nothing clamps; the pin is unchanged. Renamed guide_sd -> guide_spread_c (--guide-spread-c, DEFAULT_GUIDE_SPREAD_C 0.65, valid in [0, 2)). The CLI and main.nf refuse the old name, because an old 0.13 read as c would shrink the spread fivefold. The R reference files and SLURM scripts on this branch are brought to r-implementation HEAD (all identical at base). methods.md takes the same sections as r-implementation; usage.md and benchmark_simulation.py are updated. Tests: es = 0 gives exactly 1; the draw's moments match c * es * (1 - es) at es 0.05 / 0.15 / 0.5; c = 0 gives exactly 1 - es; the fake generators draw via beta(). 50 unit tests, 2 realdata tests, 1,247 R assertions and the two-sample stub DAG pass. A real moi5 cis run through the new draw writes 225 rows at ~0.09 s per pair per replicate. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
`g in set(pairs["response_id"])` inside the comprehension rebuilt the set for every gene: 38,606 genes x 742,525 pairs on moi5 trans, measured at 167 ms per build, 108 min in all -- past PREPARE_SIM_INPUT's 1 h limit. cis hid it at 33,066 pairs. The output is unchanged; only the time. Also let compare_stage_b.py read the pipeline's published Parquet as well as a CLI run's TSV. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
`sim_list=$(ls sim/*.parquet sim/*.tsv.gz sim/*.tsv)`: only one pattern ever matches, so ls exits 2 for the others, and under Nextflow's `bash -e` the task died on that line with empty stdout and stderr. Every real Python run hit it; moi5 cis (intergalactic_liskov) failed there after all 1,000 simulations and the consolidation had succeeded. The stub never runs the script, and R's module escaped only because it piped ls into paste. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
33,066 pairs, 100 replicates each: 96.9 % within 2 se, 98.1 % agree at the 0.8 bar with the crossings symmetric (320 / 318). The mean shift, -0.00055, is formally outside its 2-se bound; it is flat across expression and sits in mid-power pairs, the shape of R's hoisted null model against Python's refit. Also records the two bugs the first real Python runs found, and the measured cloud cost that sets trans at cpus = 4. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
pysceptre's n_jobs parallelises over one call's genes, and a cis target has a median of 6, so extra cores sat idle. The task now maps (target, replicate) units over --n-jobs worker processes: forked on Linux, spawned elsewhere (pysceptre keeps discovery state in a module global, so threads collide). Each unit draws from its own seeded stream and results come back in submission order: 1 and 8 workers write byte-identical files, and 8 ran 5.6x faster than 1 on a laptop. New realdata test pins the invariance. POWER_SIMULATION now asks for params.power_simulation_cpus = 8 and power_simulation_memory = 8 GB, replacing the R-calibrated floor + slope per MiB of split. 8 GB is measured (pysceptre-paper: day0 at 8 cores peaks at 5.13 GB with permutations, 7.24-10.3 GB with the CRT), and 4 GB would save nothing: no predefined 8-vCPU machine has less than 8 GB. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
…faster driver Section 13 records what was decided on 2026-09-25: an --estimand fixed|random option in both implementations (skip the pin; the Beta spread is zero at es = 0, so the null arm is identical under both), one permutation set per target drawn from the seed, and the faster driver, task shape and resources aimed at cis + trans in under an hour. The simulation keeps the full design. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Item 6 of section 13: confirm the enhancer bins that fixed c = 0.65 had sampling noise removed, and refit c and the form if not, before the random estimand and the PerturbPlan comparison lean on 0.0829. Not run yet. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
…from the seed Every replicate of a target gets the same pysceptre seed, from a stream spawned off the target's setup seed, so all its replicates are tested against one permutation set -- close to sceptre, whose sampler reseeds mt19937(4) on every call. The default stays per-replicate for now so the check against it is clean. Counts are identical in both modes (equal fold changes); the p-values differ; chunking and worker-count invariance hold. Wired through params.permutations and the POWER_SIMULATION module. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
The enhancer bins were noise-corrected with the same calibrated correction as the null pairs (methods.md did not say so), and the curve runs high at large effects: it peaks at es 0.5 while the bins peak near 0.37, and the bin-optimal c is 0.593 against 0.65. Item 6 now asks whether the calibration is adequate and whether the form should fall faster. Still not run. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
A worker killed mid-task (out of memory, or a crash as the pool starts) broke the pool and Python exited 1, which conf/base.config reads as a code error and stops the run on. It now exits 137, so the task is retried with double the memory. Seen on 1 of 200 moi5 cis tasks, 26 s in (2026-09-25). POWER_SIMULATION exports OMP/OPENBLAS/MKL_NUM_THREADS=1: the task's parallelism is its worker processes, and thread pools in the parent at fork time oversubscribe and can crash children. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Item 3a: in one-target calls pysceptre computes the permutation nulls by a scan that gathers and cumsums a 227 MB array to read one column; its own sparse route gives identical p-values (1,200 cis and 2,550 trans pair-tests, max |dp| = 0) at 17.7 ms instead of 119.5 ms per escalated pair. Steer the current engine onto it now, then build the matrix once per target in the faster driver. Also records the measured cost breakdown and fit reuse. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
… nulls In one-target calls pysceptre's default prefix scan builds a running sum over a (B, n_trt, 14) array to read one column. Returning None from PermutationPrefixSums.statistics -- its documented "use the draw matrix" signal -- selects its own sparse route for the duration of each call, without changing pysceptre. Measured on the real engine: cis target 53.6 -> 27.4 s (1.96x), trans target 97.4 -> 46.3 s (2.10x), output byte-identical; realdata tests pin the identity in both permutation modes. Default stays 'scan' until the running speed benchmarks finish. Upstream: broadinstitute/pysceptre#2. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
…eptre's release Item 3c: a Python FIT_NULL_MODELS step (one fit per gene per simulation on an independent es = 0 draw, reused across targets and effect sizes), default --null-fits reuse, refit kept exact. Measured 1.24-1.29x, call rate unchanged. Item 3d: pysceptre dev 0de1bff picks the sparse route itself for one target; bump the pin and drop --nulls once it is released. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
…s and simulations at once Per target, the permutation draws are made once (the per-target seed, so they equal what each engine call draws) and each stage's sparse matrix is built once and applied to every gene x simulation in one product; later stages only for the columns that escalate. Genes are fitted and tested by pysceptre's own functions (fit_all_genes, run_low_level_test_full with the precomputed nulls), so output is byte-identical to the engine run with --permutations per-target --nulls sparse: checked on the cis and trans reference targets, the largest cis target (607 cells) with a clamped-theta gene, at 1 and 4 workers, plus unit and realdata tests. cis reference target, 100 simulations: 125.7 -> 84.3 s at 1 worker, 40.7 -> 27.6 s at 4. Opt-in; requires --permutations per-target. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
…red across targets watteg-fit-null-models draws each gene's counts with no knockdown from a stream keyed rng_for(seed, "__null_fit__|" + gene, simulation, 0) and fits pysceptre's own null model to them (fit_all_genes, batch width 1, shared x_outer_flat). The fast driver takes those fits in place of refitting under --null-fits reuse, from --null-fits-file or, without one, from the same keyed fits made in the task; --null-fits refit keeps the exact per-target refit. A FIT_NULL_MODELS process between PREPARE_SIM_INPUT and POWER_SIMULATION feeds every simulation task; the file records the seed, the baseline and a digest of sim_input.h5, and a run that does not match is refused. The worker pool moves to watteg.workers so both steps share it. Checked: refit is byte-identical to the previous fast driver and engine output (cis reference target at es 0 and 0.15, trans reference target); a fit depends on (gene, simulation) only (a file for more genes and simulations, and fits made in the task at 1 and 2 workers, give the same bytes). Reuse against refit, same counts and permutations: cis 808 vs 805 calls of 1,200 (McNemar p 0.45), trans 1,280 vs 1,282 of 2,550 (p 0.73), no pair beyond its Wilson half-width; worker time 1.6-1.7x lower. Default unchanged (refit) until the new defaults land. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
…aged over it effect_size_matrix(estimand="random") skips the pin, leaving the element's realised mean over its perturbed cells where the guide draws put it; control cells stay exactly 1 and both estimands draw the same guide effects from the same stream. The option runs through the engine, the fast driver, the CLI and Nextflow (estimand), and is written as the last column of every per-simulation row, into the power table and the summary; the power steps refuse input that mixes estimands. Checked: es = 0 gives identical rows under both (fixture and CLI), fixed output is the previous file plus the new column byte for byte, and under random the realised mean's sd matches c*es*(1-es)*sqrt(sum n^2)/sum n on the shared fixture. cis reference target, es 0.15: 808 vs 781 calls of 1,200. Python only; the R side is still to do. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
…rse nulls, fast driver, fit reuse --permutations per-target, --nulls sparse, --driver fast and --null-fits reuse are now the defaults of watteg-run-power-simulation, nextflow.config and config/config.yml; the engine, refit, scan and per-replicate stay selectable, and an impossible combination is refused up front by the CLI and by main.nf. The fast driver runs the permutation test only, so a CRT screen needs --driver engine --null-fits refit; config/test.yml (the synthetic CRT object) says so. methods.md gains a section on how each simulation's test is laid out (one permutation set per target, the sparse route, the fast driver, fit reuse); usage.md documents every option. A realdata test pins that no flags at all gives the fast configuration's bytes; the stub DAG passes with both drivers. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
…lation table optional When a task holds all simulations of its pairs (reps_per_chunk == num_replicates, the default), POWER_SIMULATION writes per-pair counts (--partials-out: simulations called and used, the fold-change and cell-count sums, the simulation range) and COMPUTE_POWER adds them up (--partials). watteg.power now builds every table from counts (power_counts, merge_power_counts, power_from_counts), so the counts' table and the rows' table are the same bytes: a mean is a sum over a count either way, as pandas' grouped mean is. The per-simulation rows and CONSOLIDATE_REPLICATES run only under keep_per_simulation (default false) or when simulations are chunked across tasks; merging refuses overlapping simulation ranges or mixed thresholds, and the task checks it wrote every pair and every simulation. Fixes the rows' route: pandas' default CSV parser returned a neighbouring float for about half of a TSV's p-values and fold changes, so simulation files are now read with float_precision="round_trip". Checked in unit and realdata tests, and on a real local Nextflow run on moi5 cis (power from counts == power from the published Parquet); the stub DAG passes with and without rows, chunked, and with FIT_NULL_MODELS; day0 (CRT) passes with --driver engine --null_fits refit. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
…ands A 330-pair moi5 cis task (29 targets, 100 simulations, reused fits, per-pair counts) at 8 workers peaked at 3.38 GB in a Linux container (cgroup memory.peak, forked workers), 278 s. 8 GB is 2.4x that and the smallest predefined 8-vCPU machine, so power_simulation_memory stays. FIT_NULL_MODELS peaked at 3.04 GB on Linux at 8 workers (6.9 GB on macOS, 208 s for 244 genes x 100 simulations). On macOS, where workers are spawned and each loads its own inputs, the same task's tree peaked at 9.7 GB; three 4-worker tasks at once (4.2 GB each) beat one 8-worker task by 1.22x, because a task's wall time is bounded by its largest target. trans is not measured. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
… state
compute_null_fits used SHARED.setdefault("sim", sim), so a process whose
shared state already held another prepared input would fit that one while
recording the passed input's fingerprint. Assign it instead; forked workers
inherit it and spawned ones reload from --prepared, as before.
Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
cis ~1 h and $2.48, trans 1 h 52 min and $67.57 (previously 5.6 h and $411), power unchanged against the previous runs at the same seed. Item 1 is done in both implementations. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
On the fast configuration, cis: 97.0 % within 2 se and 98.0 % same 0.8 call in all three comparisons (two seeds, and each against R), and the mean shift is inside its 2-se bound at both seeds. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
…the Python tree main becomes the Python (pysceptre) pipeline. The R implementation lives on r-implementation, which already contains every commit main had. Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
What this does
Replaces the R pipeline on
mainwith the Python (pysceptre) pipeline fromfeat/pysceptre-backend. The R implementation stays onr-implementation.Evidence
docs/pysceptre-backend.md, "Stage B at full scale".Merge conflicts: 14 files, both sides changed
maincarries the R pipeline, and this branch the Python one. The resolution is to take this branch's version of each conflicting file, and to drop the R-only modules (fit_null_models.nf,merge_null_models.nf), which the Python pipeline replaces. That choice belongs to the owner, so this PR does not resolve the conflicts itself.Still open, not blocking
--nullsafter pysceptre's next release (Permutation nulls: prefix-scan route is 6-7x slower than the sparse route for single-target calls (identical results) broadinstitute/pysceptre#2).🤖 Generated with Claude Code