Status and handoff

State of the refactor, written so the work can be picked up by someone else.

Short version: the five pipeline steps are complete, run standalone, and have now run on Sherlock as well as on a laptop. workflow/slurm_executor/ holds one plain sbatch script per step — deliberately no orchestration layer, so a failure is attributable to the step rather than the runner. The full simulation array at effect size 0.15 ran on 2026-08-13 as job 38849611, 1,000 tasks at 100 replicates, 999 of them completing; it is the source of every measured number below.

Since then, and not yet reflected in the older sections of this document, three things landed: null_fit is implemented (Step 6a, below), the sceptre pin moved to 0.99.0, and sceptre now carries a local patch that reuses the gRNA precomputation across replicates. All three are in flight as of 2026-08-13 afternoon: the null-model array 38879935 is running, its merge 38879939 and the null_fit rerun of the simulation 38882207 are queued behind it.

The old-vs-new comparison has now run and passed (job 38916341) — the refactored pipeline agrees with the old one on every check, including a zero difference in perturbed cells per pair across all 34,886 pairs. See Comparison 1 below. Two equivalence checks had already removed the other confounders: the version bump and the gRNA patch each provably change no number.

Still missing: a test suite (Step 8). Nextflow (Step 7) is in progress on the nextflow branch — four of seven processes run, and PREPARE_SIM_INPUT’s output is byte-identical to the sbatch runner’s. Step 9 has been retired: the gRNA precomputation reuse already collected most of what it was for.

Done

The pipeline

Five standalone executables in src/, each with --help, none referencing snakemake@:

prepare_sim_input.R -> split_pairs.R -> run_power_simulation.R -> compute_power.R -> summarize_power.R

Shared code lives in lib/ at the repo root (cli.R, sim_input.R, pert_input.R, simulate.R, sceptre_io.R, apply_patch.R), located from the calling script’s own path — each executable resolves dirname(dirname(<its own path>))/lib — so it works both standalone and staged onto PATH by a workflow engine.

Layout changed on 2026-08-13, ahead of the Nextflow work: bin/ became src/, and bin/lib/ was lifted to a top-level lib/. The old Snakemake implementation was removed at the same time and lives on the legacy branch. Anything referring to bin/ predates this.

Measured improvements

Correctness fixes

Two of my own claims that measurement overturned

• Control-cell sampling is not a usable speed lever. It costs 21–60 % of power. Off by default. Table in Methods.
• The speedup from smaller matrices was ~7×, not the ~109× the size reduction suggested. Per-replicate cost is sub-linear in cell count; ~0.3s per replicate is fixed overhead.

Infrastructure

• pixi.toml / pixi.lock — three environments (runtime / build / dev), sceptre pinned to 0.99.0 (3ba046b8, an untagged main commit — v0.10.3 is still the newest tag, so SCEPTRE_REF names a DESCRIPTION version, not a tag) plus a local patch, both guarded by pixi run check-api. See the sceptre section below
• .githooks/pre-commit — auto-fixes line endings, tabs (2 spaces for R/YAML, 4 otherwise), trailing whitespace, final newline; rejects files >512 KB, non-snake_case filenames, camelCase R identifiers. pixi run lint is the read-only repo-wide version
• Documentation published to GitHub Pages
• config/config.yml — parameters only; dead reps key removed

Reference numbers

Needed to size cluster runs. Measured, not estimated. The right-hand column is Sherlock, on day0_grna20_no_shuffle, and supersedes the laptop figures wherever they differ.

The ~295 CPU-hour figure was a laptop measurement and is roughly 3× optimistic. Sherlock’s cores are slower per-thread than an M-series laptop, and sceptre’s per-call cost is untouched by the refactor, so it dominates. Sized against the old pipeline’s own sacct record for its 2026-05-14 run, the refactor is still ahead on both axes:

Both ran 100 replicates on the same dataset, so this is like for like. The CPU win is 1.6×; the wall-clock win is mostly the unthrottled array rather than the code.

48 % of the cost is the per-target term — one run_discovery_analysis() call carries the full 586,309-cell bookkeeping however few gene pairs ride along. Targets can only share a call when their gene sets are disjoint, which is Step 9 below. Consequences worth internalising before optimising:

• Halving the pairs saves ~26 %, not 50 %. Halving the replicates saves exactly 50 %.
• Filtering pairs for a second pass is much dearer than it looks: the 4,322 pairs whose interval straddles power 0.8 (12.4 % of pairs) are spread across 1,820 of 3,026 targets, so a top-up run pays 60 % of the per-target overhead to redo 12 % of the pairs — 82 % of its cost is overhead.
• Effect sizes cost the same as each other. At the current 635 CPU-h a six-point sweep is ~3,800 CPU-h (it was ~5,100 before the gRNA reuse), and the 2,900-task QOS budget still allows one wave: 480 splits × 6 effect sizes = 2,880 tasks, ~75 min each, so ~1.5 h wall clock.

Settled — which per-gene null model does the simulation test against?

Verdict: fit the null model on a null simulation of each replicate and reuse it (null_fit). Today’s inherited-cache behaviour understates power, and the fix costs 1.7 %. Measured by job 38854980 and the threshold check 38861277; the numbers are below.

Before testing any pair, sceptre fits a Poisson GLM of the gene’s counts on the cell covariates (perform_response_precomputation, precomputation_functions.R:13) — the null model the perturbed cells are compared against. It skips that fit when @response_precomputations already holds an entry for the response_id (medium_level_functs_v2.R:63, reached because run_outer_regression <- calibration_check || control_group_complement and our control_group_complement is TRUE).

Two measured facts follow:

• glm.fit is 27 % self / 59 % total under Rprof, against 7 % for our own rnbinom draws, so GLM fitting is the single largest cost in a call. That 59 % is not this precomputation. The profile was taken over the current path, with the inherited cache populated and the precomputation therefore skipped — so the 59 % is GLM work sceptre does anyway. Clearing the cache adds the precomputation on top of it, and the table below measures that: 4.3×, all of it in the per-pair term. What the resident 59 % actually is remains unidentified, and it bounds what any pooling optimisation such as Step 9 can recover.
• sceptre_template.rds inherits 272 precomputations from the real discovery analysis, and they cover all 237 genes that have QC-passing pairs — none are refitted. So every simulated count is tested against coefficients fitted to real counts. Measured directly by job 38854980: 237 of 237 present, 0 absent. The bias is therefore uniform rather than mixed across genes, which is easier to reason about but no more correct — a faithful emulation fits its null model on the data it is given, and the real analysis did exactly that.

The three configurations, measured

Job 38854980 ran all three on 3 targets × 5 replicates, rotating the order within each replicate so no configuration is always first:

cleared is 4.3×, and the extra cost sits in the per-pair slope, not the intercept — each pair in a call is a distinct gene, so refitting scales with the genes tested. It is a reference, never a production candidate. null_fit buys its fidelity for +1.7 % over as_is.

Treat the ratios as sound and the absolute CPU-hours as indicative only: this job timed as_is at 11.5s on the 5-pair target where the profiling run measured 6.15s, so its clock runs ~1.9× slow — probably because it holds three caches and a null simulation in memory at once. That is also why the 1,305 CPU-h here disagrees with the ~858 in the table above. The 100-replicate array’s own sacct record supersedes both.

Why as_is is not good enough

Aggregate p-value agreement looks reassuring and is misleading:

cleared vs null_fit   n=265  max|diff|=0.012  median=1.2e-10  spearman=1.0000
cleared vs as_is      n=265  max|diff|=0.08   median=4.6e-05  spearman=0.9990

Power is not aggregate. It is the fraction of replicates with p < 7.2629e-4, so only crossings of that threshold matter, and a median difference of 4.6e-05 is ~6 % of the threshold itself. Counting the calls instead of correlating the p-values (job 38861277, threshold_check.R):

cleared vs null_fit   flips=0 / 265 (0.00 %)   per-pair power identical
cleared vs as_is      flips=7 / 265 (2.64 %)   [cleared-only=7, as_is-only=0]

• null_fit is exactly equivalent to the faithful reference — zero disagreements in 265 calls, and per-pair power identical to three decimal places.
• as_is is biased, one-directionally. All 7 flips go the same way: cleared calls the pair significant and as_is does not, never the reverse (137 of 265 significant against 144). Random noise would split both ways — P ≈ 0.016 for 7–0 — so as_is understates power, by a mean |Δpower| of 0.026 on this sample.
• It concentrates where it matters. Only ~37 of 265 p-values sit within a factor of 10 of the threshold, so roughly 19 % of near-threshold calls flip.

Caveat on the magnitude: 3 targets, 53 pairs, 5 replicates. The direction is solid and the equivalence of null_fit is unambiguous; the 2.64 % figure is thin and should not be quoted as precise.

Superseded on the magnitude. The full 100-replicate run under both configurations (see Consequences below) puts the mean shift at +0.0063, not the ~0.026 this sample gave — the direction held, the size did not. The 2.64 % flip rate was measured on 5 replicates, where a single replicate crossing the threshold moves power by 0.2; at 100 replicates it is 0.8 percentage points of pairs crossing the certification line.

Note the coefficients themselves differ substantially — inherited theta median 20.0 against 18.8 from the null fit, median worst-coefficient difference 7.77 — while the p-values barely move. sceptre’s conditional-resampling test is far more robust to the null model than the null model is stable.

Consequences

• The effect size 0.15 array (38849611) ran as_is, and the bias has now been measured directly rather than extrapolated. Both configurations exist at 100 replicates over all 34,886 pairs (power_as_is/ and power_null_fit/), and they are exactly paired — same simulated counts, so every difference is the null model with zero Monte-Carlo noise between them.

  	as_is	null_fit
  Mean power	0.5933	0.5996 (+0.0063)
  Certified (power_ci_low ≥ 0.8)	13,094 (37.5 %)	13,374 (38.3 %)
  Ambiguous	4,382 (12.6 %)	4,408 (12.6 %)
  Underpowered	17,410 (49.9 %)	17,104 (49.0 %)

  The direction is confirmed and the magnitude is ~4× smaller than predicted. as_is understates power: 11,649 pairs rise against 3,631 that fall, a 3.2:1 ratio that noise would not produce. But the shift is +0.0063 mean, 0.0100 mean absolute, against the ~0.03 this document extrapolated from 3 targets × 53 pairs × 5 replicates. That extrapolation was too thin, as it said at the time; treat 0.0063 as the number and the earlier 0.03 as withdrawn.

  It concentrates exactly where the mechanism says it should — in the transition band, not at the saturated ends:

  as_is power band	pairs	mean Δ	% moved
  0.00–0.02 (floor)	4,411	+0.0010	9.7 %
  0.02–0.20	4,409	+0.0094	52.7 %
  0.20–0.50	4,411	+0.0121	73.4 %
  0.50–0.80	5,894	+0.0130	73.5 %
  0.80–0.98	8,937	+0.0047	48.6 %
  0.98–1.00 (ceiling)	6,824	+0.0003	9.0 %

  Practical consequence: null_fit is the correct configuration and should stay the default, but nothing built on the as_is numbers was materially wrong — it moved 280 of 34,886 pairs across the certification line, 0.8 percentage points. Keep both outputs.

• The pre-refactor pipeline inherited the same cache, so the old-vs-new comparison is unaffected: both sides share the bias. Compare them as_is against as_is.
• Clearing the slot without supplying null fits gives cleared, at 4.3×. The two changes must land together.

An earlier version of this section put the cost of cleared at ~2.4× per call and cited a p-value divergence of 0.942 with Spearman 0.224 between as_is and null_fit. Both came from profile_call.R, whose first configuration was labelled “refit, as the pipeline runs today” but inherited the template’s cache and refitted nothing (line 100 only overwrites the slot when precomp is non-NULL). The 2.4× was 1 / (1 - 0.59) applied to a glm.fit share that is not the precomputation. The proper experiment gives as_is vs null_fit max 0.088 and Spearman 0.9989, so both numbers are retracted. Why that script reported 0.942 at all is unexplained; it is superseded by cache_experiment.R and should not be used.

Step 6a — null_fit: implemented, running

Not a one-line change, because the null fit has to be shared across array tasks. A gene’s null model is target-independent, so one fit per (gene, replicate) serves every target — but targets are spread over 1,000 splits, so fitting inside each task would pay for it 1,000 times instead of once. It is also effect-size independent: a null simulation has no knockdown, so one set of fits serves the whole sweep. 100 replicates, not 100 × 6.

What shipped:

1. src/fit_null_models.R — for each replicate, simulates a null matrix (no knockdown) and runs one run_discovery_analysis() with @response_precomputations empty, keeping only the resulting slot. An entry is 11 named coefficients plus a theta scalar, so 100 replicates × 272 genes is a few hundred KB. Seeded from (seed, rep) only — deliberately not (seed, target, rep, effect_size), since the point is that it depends on neither target nor effect size.
2. run_power_simulation.R --null-precomputations assigns that replicate’s entry to @response_precomputations inside the replicate loop (line 317), and validates that the bundle’s seed and replicate range match the run’s.
3. The inherited slot is cleared in slim_sceptre_object() (lib/sceptre_io.R:199), after build_dispersion_vector has read it (prepare_sim_input.R:259) — dispersions must keep coming from the real data, since they set the noise the simulation is meant to reproduce. Only the null model moves. Clearing is what makes an accidental as_is impossible to reintroduce.
4. Steps 02b (array, one task per replicate), 02b_submit.sh and 02c (merge) in workflow/slurm_executor/, with --null-precomputations threaded through 04_run_power_simulation.sbatch.

The guard is in. run_power_simulation.R:166 errors if the template arrives with a non-empty @response_precomputations and no --null-precomputations — rather than silently running as_is. That is the bug that went unnoticed for the entire refactor, and it can no longer recur silently.

Done as of 2026-08-13. The null models are built (38879935, merged by 38879939 into prepared/null_precomputations.rds), and the null_fit rerun of effect size 0.15 completed as 38887744: 1,000/1,000 tasks, every output at the exact expected row count, every output carrying a matching .provenance sidecar. Results are in results/refactor/sim_null_fit/es0.15/. The first attempt, 38882207, silently skipped all 1,000 tasks and was cancelled — see the resume-check entry under Gotchas.

And it is 1.58× cheaper than the as_is run, not 1.7 % dearer.

The null model itself costs ~1.7 % more, so the saving is not from it — it is the gRNA precomputation reuse, which shipped in the same commit. That retires the “the saving is small” caveat recorded under the equivalence check below: 1.02× was measured on the single worst-case split (1 target, 36 pairs, in the cleared configuration), and at the median 9 pairs per target the fixed per-target gRNA fit is a far larger share of the call. 1.58× at production scale is the number to use; the step-10 figure measured what it measured, on a split chosen to be pessimal.

Note this also invalidates the cost model in the reference table for any future run: 5.076s + 0.5923s × pairs was fitted on as_is output with no gRNA reuse. Refit it on 38887744 before sizing the next sweep. 02b_submit.sh had never run successfully before commit 08b0d04: PASSTHROUGH[@]: unbound variable under set -u, the bash 4.2 empty-array bug already fixed in 04_submit.sh but never back-ported. That is why null_models/ was empty.

Keep the as_is output rather than deleting it — it is the only direct measurement of how much the bias moved real numbers at 100 replicates, and the threshold check that found the bias used 5.

sceptre: pinned to 0.99.0, and patched

Two changes to the dependency, both verified not to move any number.

The version bump

The pin moved from v0.10.3 (9faa4373) to 0.99.0 (3ba046b8, 2026-08-07). There is no v0.99.0 tag — v0.10.3 is still the newest tag — so SCEPTRE_REF names an untagged commit’s DESCRIPTION version. The slot set is byte-for-byte the same 55 slots (reformatted only), so sceptre_object.rds and sceptre_template.rds serialised under 0.10.3 still deserialize, and 0.99.0’s NEWS is Bioconductor prep with no statistical change.

Two things the bump forced, both easy to miss:

• cat() became message(). capture.output() only catches stdout, so the “consider parallel = TRUE” note would have leaked ~279,800 lines per effect size. suppressMessages() now wraps the run_discovery_analysis() calls in run_power_simulation.R and fit_null_models.R; both wrappers are kept so the pin can move either way.
• New Imports parallelly and withr, both on our code path. R CMD INSTALL enforces Imports, so they had to enter pixi.toml, forcing a pixi.lock re-solve — which must run inside the ubuntu container, since the host is CentOS 7 (glibc 2.17) and pixi solves against the host’s virtual __glibc. That is 00b_relock_env.sbatch (job 38877284); the only package added was r-parallelly 1.48.0. Commit pixi.toml and pixi.lock together.

Verified not to be a problem: the one RNG change on the association path (set.seed(4) → withr::with_seed(4, ...) in partition_response_ids()) sits inside if (parallel), and the pipeline always calls parallel = FALSE.

The patch — reusing the gRNA precomputation

sceptre’s CRT path regresses perturbation status on the covariates and draws synthetic treatment assignments from the fitted probabilities. That fit depends on the gRNA-to-cell assignments and the covariate matrix — never on the response counts. The simulation calls run_discovery_analysis() once per (target, replicate) with only @response_matrix changing, so every replicate after the first refit an identical model.

Unlike @response_precomputations there is no slot to hand it back through: fitted_probabilities is a local variable inside the two CRT workhorses. So this ships as patches/0001-reuse-grna-precomputation.patch, adding an optional argument with a NULL default — unpatched behaviour is unchanged. Reuse is on by default; --no-grna-precomp-reuse turns it off.

Three design points are load-bearing — do not “simplify” them away:

• The cache holds coefficients, not fitted values. Fitted values are one double per cell, which at this scale is ~2.3 GB held in the caller’s session. Reconstruction via binomial()$linkinv(drop(X %*% coefs)) is bit-identical — glm.fit computes the fitted.values it returns exactly that way, and the offset is zero here. Verified: identical() TRUE, max abs diff 0.
• Each entry records n_cells, and the guard checks it. The two workhorses fit on different cell sets but the same covariates, so the coefficient vector is the same length in both — length alone cannot tell a cache built for one path from a cache built for the other.
• Never cache crt_index_sampler_fast()’s output. The fitted probabilities are deterministic; the synthetic index sets drawn from them must be redrawn every call. Caching the indices would freeze the resampling distribution across replicates and destroy the CRT null — producing a power estimate that is garbage but looks plausible.

A target that is rank deficient within its own cells is omitted from the cache with a warning: glm.fit returns NA for an aliased coefficient while keeping fitted.values valid, so a naive X %*% coefs would yield all-NA.

Supporting machinery: lib/apply_patch.R is a pure-R unified-diff applier, needed because the cluster runs everything in a stock ubuntu:22.04 container that has neither patch nor git. It is stricter than patch(1) — no fuzz, no offset search, exact context match — so a moved pin fails loudly. src/install_sceptre.R stamps RemoteSha and a LocalPatches md5 into DESCRIPTION and short-circuits only when both match; src/check_sceptre_api.R asserts the patched API separately, because an unapplied patch would otherwise surface only at the first run_discovery_analysis() call, long after steps 01 and 02.

Equivalence check 10_grna_precomp_equivalence — does the patch change results? No. And the saving is small.

Job 38879261: outputs byte-identical between reuse and refit. That was the point of the job and it is confirmed.

The cost result is weak, and much weaker than the retracted reasoning implied:

refit every replicate   620s   (124.0s per call)
reuse per target        608s   (121.6s per call)
speedup                 1.02x  (-1.9% wall clock)

split_0001 is 1 target with 36 pairs, so 5 replicates = 5 calls and reuse removes 4 of 5 gRNA fits: 12s / 4 ⇒ ~3s per gRNA fit. Two caveats pointing opposite ways: 1.9 % on 620s is within plausible node-to-node variation, so 3s is an upper bound from a noisy difference, not a measurement — and this is the least favourable target in the dataset, since 36 pairs is the maximum against a median of 9, and the run was in the cleared configuration, which inflates the per-pair term ~5.75×. Against the null_fit model a fixed ~3s is ~8 % at 36 pairs but ~22 % at the median 9.

The “up to 2.4×” claim stays retracted. The replacement is 1.58×, measured at production scale — see the Step 6a table above, where the full null_fit array came in at 634 CPU-h against the as_is array’s 999. Both caveats above pointed that way and the second one was right: this split is the least favourable in the dataset. Do not quote the 1.02× as the saving; quote it as what the equivalence check happened to cost on one worst-case split, which is not what it was measuring.

Equivalence check 11_sceptre_version_equivalence — does the version bump change results? No.

Job 38879265: stock v0.10.3 against patched 0.99.0, reuse off both sides, identical — the bump does not change any number, so the old-vs-new comparison (Comparison 1 below, 07_compare_old_new.sbatch) is not confounded by it.

Both checks are exactly paired: seeds derive from (seed, target, rep, effect_size), which involves neither the cache nor the sceptre version, so the simulated counts are bit-identical across each pair and any difference would be attributable. Between them they isolate the two confounders the old-vs-new comparison would otherwise carry.

Gotchas found along the way

Recorded because each cost real time and none is discoverable from the code.

• vendor/sceptre.tar.gz is a snapshot of main, not the pin, and is untracked (vendor/ is gitignored). Do not use it as an install source.
• The login node caps open file descriptors at 256, hard, and pixi lock exhausts it mid-solve (No file descriptors available (os error 24)). The solve must run in a job — hence 00b_relock_env.sbatch. More generally, nothing heavy belongs on the login node.
• pixi cannot install R packages from GitHub at all — conda channels only, and its git support is PyPI-only. Hence the separate pixi run setup task plus a SHA pin.
• .githooks/pre-commit strips trailing whitespace from source extensions, which would corrupt a .patch file’s blank context lines. .patch is not in its is_source_file list, so it is safe today — and apply_patch.R also treats a zero-length line as blank context. Do not add .patch to that list.
• A patch’s context only guards the lines it touches. The API assertions in check_sceptre_api.R are the layer that actually catches a silently-unapplied patch.
• 02c’s dependency is afterok:38879935,afterany:<step 10>:<step 11>, and the afterany is deliberate: 02c writes prepared/null_precomputations.rds, and the two equivalence checks test for that file at runtime, so an unguarded merge would silently flip them from cleared to null_fit depending on scheduling.

• Step 04’s resume check skipped all 1,000 tasks of the null_fit rerun. Job 38882207 reported 1,000 successes in ~15 s each and recomputed nothing: out was ${SIM_DIR}/es${ES}/${SPLIT_NAME}.tsv with no reference to the null-model configuration, and “complete” was judged by row count — which is identical under every configuration, since the null model does not change how many rows come out. So the as_is output from 38849611 satisfied the check. The R-level guard at run_power_simulation.R:166 cannot catch this, because R never starts. Worse, the skip is the only thing that prevented a second bug: step 04 now passes --null-precomputations unconditionally, so had it run it would have overwritten the as_is baseline in place — the same directory 09_compare_null_fit.sbatch requires as the as_is side.

  Fixed two ways. NULL_MODEL_CONFIG is now part of the output path (SIM_DIR=${OUTDIR}/sim_${NULL_MODEL_CONFIG}, POWER_DIR likewise), so two configurations can no longer address the same file; the as_is baseline moved to results/refactor/sim_as_is/ and is named explicitly by AS_IS_SIM_DIR. And each output now carries a .provenance sidecar recording the configuration, the null-bundle md5, reps, seed, guide-sd and effect size — the resume check requires it to match, so an output produced differently is redone rather than trusted. A missing sidecar means the file predates the check and cannot be vouched for, so it is also redone. Resubmitted as 38887744.

• Step 03’s replicate-count check compared against stale output. It failed with “971 rows, 376 pairs, 109 with != 2 reps” while the three files the run actually wrote were each exactly 36 pairs × 2 reps. The directory held output from an earlier 280-split run, and because split_pairs.R pads the index to the width of --n-splits, those files are named smoke_split_001_* against smoke_split_0001_* — the check’s glob could not tell them apart. Fixed in 628d5c9: step 03 now records the outputs it writes and checks only those.

Left to do

Step 7 — Nextflow + SLURM profile — in progress, on the nextflow branch

Four of seven processes exist and run: PREPARE_SIM_INPUT, SPLIT_PAIRS, FIT_NULL_MODELS, MERGE_NULL_MODELS. Still to write: POWER_SIMULATION, COMPUTE_POWER, SUMMARIZE_POWER.

Minimal DSL2, not nf-core scaffolding — same principle as workflow/slurm_executor/: no framework layer, so a failure is attributable to the step rather than the runner. The sbatch executor stays until Nextflow reproduces a run; it produced every measurement in this document.

main.nf                            the DAG
nextflow.config                    params and profiles
conf/base.config                   measured resources, per process
conf/sherlock.config               SLURM + apptainer + pixi
config/config.yml                  production parameters
config/test.yml                    3 splits x 2 replicates
modules/local/*.nf                 one file per process
workflow/nextflow/run.sbatch       the driver, which runs as a job

The DAG, which is not the shape planned above — that predates null_fit:

samplesheet -> PREPARE_SIM_INPUT ---+-> SPLIT_PAIRS -----------------+
                                    |                               |
                                    +-> FIT_NULL_MODELS (x reps)    |
                                        -> MERGE_NULL_MODELS -------+
                                                                    v
                                    POWER_SIMULATION (split x effect size)
                                         -> collectFile by (sample, effect size)
                                         -> COMPUTE_POWER -> SUMMARIZE_POWER

FIT_NULL_MODELS hangs off PREPARE_SIM_INPUT rather than off SPLIT_PAIRS because a gene’s null model is fitted on a null simulation and so depends on neither the target nor the effect size. One set of fits serves every split and every effect size in a sweep — 100 replicates, not 100 × splits × effect sizes. Making it a sibling rather than a descendant is what expresses that.

The design question that had to be settled first

The sbatch runner does cd $REPO_ROOT && pixi run. Nextflow cannot: a task must run in its own work directory or its outputs are invisible to the engine. Two facts, both verified before building on them, make it work anyway:

• pixi run --manifest-path <repo>/pixi.toml preserves the working directory (tested from /tmp), so the task’s outputs land in the task directory.
• R_LIBS_USER is derived from $PIXI_PROJECT_ROOT, not the cwd, so the environment follows the manifest.

Every process invocation is built on that pair.

Validated at production scale

The Nextflow runner reproduces the sbatch runner exactly, on a full 1,000-split run. Job 38978285, 2h19m wall clock, all seven processes, 34,886 pairs at 100 replicates:

That is the check that decides whether the two runners can be considered equivalent, and it passes on the deliverable rather than on an intermediate. Note the per-replicate p-values do still differ in the far tail between the two runs — see “Reproducibility is bitwise on one node” — but not one of those differences reaches the power estimate.

Two incidental findings from that run:

• results/day0_grna20/sceptre_object.rds and results/day0_grna20_no_shuffle/sceptre_object.rds are byte-identical. The _no_shuffle suffix describes how the old pipeline processed the sample — with the size-factor permutation removed — not a different input. The dataset inventory under Step 11 lists them as separate datasets, which is true of their old-pipeline outputs and false of their inputs. Running the refactored pipeline on both is redundant.
• Preemption on owners is handled. Three tasks were killed with exit 143 and retried automatically by conf/base.config’s errorStrategy, without intervention. The sbatch runner needs a resubmission plus its resume check to do the same.

Verified so far

PREPARE_SIM_INPUT’s output is byte-identical to the sbatch runner’s — all five files, both .rds objects included (sim_input.rds 16,057,759 bytes; sceptre_template.rds 46,520,066). That is a much stronger result than “it ran”: a single difference anywhere in the container, pixi environment, R version or patched sceptre would perturb the serialised objects. The environment integration is therefore the same one the measurements came from.

Four traps, each of which reported something other than its cause

Recorded because every one cost time and none is guessable from the error text.

1. .ifEmpty { error ... } aborts every run. Nextflow invokes the closure while building the DAG, not when the channel proves empty, so a guard for a malformed samplesheet fired unconditionally — and, being the first error raised, printed “samplesheet has no rows” in place of the real failure further down, masking trap 2 completely. Validate the file eagerly instead.
2. take rejects a def local or an as int cast. Both arrive wrapped in a PojoWrapper and fail with Missing process or function take([DataflowStream, PojoWrapper]), which reads like a missing operator. A literal or a params.* value works; nothing else tried did. Bound the range when you build the channel rather than trimming it afterwards where you can.
3. -profile test was silently doing almost nothing. Nextflow gives -params-file precedence over profile params, and the driver always passes one, so every test override sharing a name with a production parameter was discarded — the run used 2 null chunks of 1 replicate where the profile asked for 1 of 2, and num_replicates stayed at 100. A “test” run would have been a full-size run. Test settings now live in config/test.yml, selected with PARAMS_FILE. Profiles choose where a run executes; params files choose what it runs.
4. A local def scratch is shadowed by the scratch process directive. conf/sherlock.config sets process.scratch = false, so "--bind ${scratch}:${scratch}" expanded to --bind false:false and apptainer died with “unable to add false to mount list: destination must be an absolute path”. It had also silently reached PIXI_CACHE_DIR and friends. The variables are named scratchDir/groupHomeDir for that reason — do not tidy them back.

Two smaller ones: the driver must resolve REPO_ROOT from scontrol show job ... Command= rather than $BASH_SOURCE, because Slurm copies the script into its spool directory; and Nextflow’s file() resolves a relative path against launchDir and returns an absolute path, so file(p).isAbsolute() is always true and cannot be used to test whether the user gave a relative path.

Notes still standing

• The combine step needs no script: collectFile(keepHeader: true, skip: 1) replaces the deleted combine_sceptre_power_analysis.R.
• Resources come from the measured table above, not the old Snakemake guesses (which asked 8 GB for simulation and 64 GB for prepare). conf/base.config names the source of each.
• cpus 1 on the simulation process — sceptre runs with parallel = FALSE and parallelism comes from the fan-out.
• The driver runs as a job, never on a login node: it is a JVM that lives as long as the pipeline. NXF_HOME and NXF_WORK are redirected to $GROUP_HOME and $SCRATCH, since the defaults are ~/.nextflow and ./work and $HOME is 15 GB.
• nextflow run . -profile stub -params-file config/test.yml -stub-run checks the wiring in seconds without touching the cluster. The stub profile uses resourceLimits so the measured requests do not have to be rewritten to fit a login node.
• Removing Snakefile, rules/ and R/ belongs to this step. Done ahead of it, on 2026-08-13, once Comparison 1 had passed and the old implementation had nothing left to prove. Snakefile, rules/, R/ and envs/ are gone from this branch and preserved on legacy. pixi run lint is now green — those files were the only things it flagged.

Why not sceptre’s own Nextflow pipeline

sceptre now ships one (timothy-barry/sceptre-pipeline). It is the right tool for a real discovery analysis and the wrong one for a power simulation. Evaluated against a clone of it:

• It requires ondisc matrices. Every script calls read_ondisc_backed_sceptre_object(sceptre_object_fp, response_odm_fp, grna_odm_fp), with no in-memory path. Using it would mean serialising each simulated matrix to ODM.
• The simulation needs one object per gene-disjoint batch, not one per replicate — the same constraint as Step 9, since a gene is paired with ~147 targets. That is 100 × 300 = 30,000 objects for one effect size, each needing its own ODM.
• It re-runs everything upstream per object. The workflow is set_analysis_parameters -> prepare_assign_grnas -> assign_grnas -> combine -> run_qc -> prepare -> analysis, unconditionally; pipeline_stop truncates the tail but cannot skip the head. gRNA assignment alone is 43,432 gRNAs at their own 2s/gRNA heuristic, ~24 CPU-hours per object, none of it dependent on the simulated response matrix.
• Its 0.05s per pair is a walltime heuristic, not a benchmark — run_association_analysis_time_per_pair, used to size --time. It is not comparable to our measured per-pair cost.

What is worth taking from it: prepare_association_analyses.R sorts pairs by response_id and pods on that key, so all ~147 pairs of a gene share one @response_precomputations entry — one GLM fit serving 147 tests. That is the same insight as Step 9 (pool more pairs per call) and it is what led to the null-model question above.

Reproducibility is bitwise on one node, not across the cluster

Found on 2026-08-14 while checking the Nextflow runner against the sbatch one on the same three splits, at the same seed. Two of the three splits came out byte-identical. The third did not:

This is not a different RNG stream — that would perturb every row, not four of seventy-two. It is last-bit floating-point difference between CPU generations (this cluster spans six; see the node table under Open questions), amplified by sceptre’s skew-normal tail approximation, which is very sensitive in the far tail and nowhere else. A p-value of 10⁻⁸¹ is zero for every purpose the pipeline has.

What the seeding guarantees is the random draws, not the arithmetic. So:

• The same seed on the same node type reproduces a run bit for bit.
• Across node types, expect agreement to ~15 significant figures on estimates, and agreement on every significance call — which is what power is computed from, and therefore what matters.
• Do not write a test, a comparison or an acceptance criterion that demands byte-identity of simulation output. workflow/nextflow/verify_against_sbatch.sbatch compares numerically and asserts zero significance flips, which is the correct bar; an earlier version demanded byte-identity and failed on hardware rather than on the pipeline.

Step 8 — tests and comparisons

1. Unit tests — done: 62 tests in tests/testthat/, run with pixi run test or sbatch workflow/slurm_executor/12_unit_tests.sbatch (R does not belong on a login node). They cover the pure functions only — no sceptre, no lab data, no cluster — and every one of them corresponds to an entry under “Correctness fixes” above:

   File	Covers
   test-stats.R	wilson_interval against published values and at both boundaries, where the normal approximation collapses; effect_label pinning the 0.2 → "2" bug
   test-simulate.R	center_effect_size_matrix putting each gene’s perturbed mean on its target; create_effect_size_matrix orientation and clamping; build_dispersion_vector erroring rather than recycling
   test-seeding.R	derive_seed separating every key component, and being invariant to how replicates are chunked — the property that makes --n-splits a purely computational knob

   Two of these tests failed on first run because the tests were wrong, not the code: the effect-size matrix is genes × cells, not cells × genes (run_power_simulation.R:314 reorders columns by cell), and init_seed logs its seed deliberately. Both were corrected against the actual contract rather than by loosening the assertion — a test encoding a wrong assumption is worse than no test, and the transposed one would have passed a genuinely broken implementation.

   wilson_interval and effect_label moved to lib/stats.R to make this possible; they were defined inside compute_power.R and summarize_power.R, where testing them meant invoking a command line.

   Not covered yet: split binning conservation, which is enforced in SPLIT_PAIRS and 02_split_pairs.sbatch at runtime rather than by a unit test.

2. Synthetic test data generator — done: src/make_test_data.R builds a sceptre object the whole pipeline runs on, in 43 seconds, with no lab data. sbatch workflow/slurm_executor/13_make_test_data.sbatch builds it and then runs all six steps on it end to end. Latest: 360 QC-passing pairs, 27 significant discoveries, threshold 7.3 × 10⁻³, 30 precomputations, 648 KB.

   Built through sceptre’s own API — import_data → set_analysis_parameters → assign_grnas → run_qc → run_discovery_analysis — rather than by fabricating slots, so it cannot drift from the pinned library. Two constraints are not obvious and both were found by the object being rejected:

   • It has to plant real knockdowns. discovery_threshold() errors when nothing is significant, deliberately, because the old code returned -Inf there and reported zero power for every pair. Pure noise produces an unusable object, not a boring one.
   • Cells must carry a variable number of gRNAs. With exactly one each, grna_n_nonzero is constant, collinear with the intercept, and sceptre rejects the formula as containing “redundant information”. moi = "high" for the same faithfulness reason: it is what makes the control group the complement, which is the CRT path production actually uses.

   The object is not committed — tests/data/ is gitignored. It is 648 KB against the pre-commit hook’s 512 KB ceiling, and regenerating it is the point of having a generator.

   sceptredata was considered instead. It is not on conda, not installed, and not even in sceptre’s Suggests, so it would need the same SHA-pin-and-install machinery as sceptre itself; and since it ships raw matrices rather than a sceptre_object, the API calls above would still have to be made by hand. Reasonable, but it removes less than it adds.

3. Old-vs-new comparison — done and passed, job 38916341; see Comparison 1 below.
4. Two-stage replicate allocation — documented in Choosing num_replicates but not orchestrated. The scripts already support it through --rep-offset.

Step 9 — batch gene-disjoint targets into one sceptre call — RETIRED, do not build this

The gRNA precomputation reuse already collected most of what this was for, by a different mechanism. Step 9 amortises the per-call setup by putting many targets in one call; the patch eliminated the dominant component of that setup outright, by caching the fit instead. Refitting the cost model on 38887744 shows what is left:

Batching 302,600 calls into ~30,000 would leave ~10 CPU-h of per-target overhead, so it now saves about 86 CPU-h per effect size out of 635 — against the ~385 it promised. That is not worth a refactor that touches seeding, splitting and the target loop, and whose two open risks were never resolved: a 194-gene batch needs ~0.9 GB for the count matrix and as much again for the effect-size matrix against ~54 MB per target today, and sceptre’s per-call cost was only ever shown linear over 5–36 pairs, not the 194 a batch would carry.

The per-pair term is now 85 % of the bill, and Step 9 explicitly does not move it — a gene is still simulated once for every target it pairs with. Any further optimisation has to attack per-pair work, not per-call setup.

The reasoning below is kept because the analysis is still correct and the batching argument is reusable if the per-target term ever grows back.

Originally: the largest single saving left on the table, up to 1.8×, worth about 370 CPU-hours per effect size. Costed below because the reasoning is easy to get wrong in both directions.

Where the cost is. The unit of work is one run_discovery_analysis() call per (target, replicate) — 3,026 × 100 = 302,600 calls per effect size — and each one carries the full 586,309-cell setup regardless of how few gene pairs it covers. That per-call term is 4.91s of the 4.91s + 0.459s × pairs model, so it is 413 of the 858 CPU-hours: half the bill is setup paid and over.

Why the obvious fix does not work. The tempting version is one simulated object per replicate covering all 34,886 pairs, so the setup is paid 100 times instead of 302,600. It is not valid. Each gene is paired with 147 targets on average, up to 299 — only 10 of 237 genes belong to a single target — and pair (A, g) needs gene g knocked down in A’s perturbed cells while pair (B, g) needs it knocked down in B’s. One response matrix cannot hold both. Applying every knockdown at once puts B-perturbed cells, also knocked down, into A’s control group, which shrinks the contrast and understates power. The pre-refactor pipeline was per-target for the same reason; this is not an artifact of the refactor.

What does work. Two targets can share a call whenever their gene sets are disjoint, because a knockdown on gene gA never touches gene gB’s row. Cell overlap is fine: a cell with gRNAs for both A and B is in each one’s treatment group for its own gene, and unmodified for the other’s. Greedy first-fit over the real pair table gives 300 batches, against a lower bound of 299 forced by the gene that appears in 299 targets — so greedy is essentially optimal here.

Batches hold at most 31 targets, 194 pairs and 194 of the 237 genes.

The per-pair term does not move, and that bounds the win. Gene g is still simulated once for every target it pairs with, so the total rnbinom volume is unchanged. This is 1.8×, not the 10× the call-count ratio suggests.

parallel = TRUE is a different question and probably not the lever. sceptre’s internal parallelism forks across pairs within a call; the array fan-out already provides that parallelism across nodes, without fork overhead and without pinning a replicate to one node’s cores. It does not reduce the number of calls, which is where the waste is. Batching and who-schedules-the-parallelism are orthogonal.

Three things to check before building it.

1. Settled by the array: the per-call overhead is 427 CPU-h and Step 9’s ceiling is 1.63×. Fitting time = reps × (a × targets + b × pairs) over all 999 completed tasks gives 5.076s + 0.5923s × pairs per (target, replicate) — 1,001 CPU-h predicted against 999 measured, so the model is sound. The per-target term is 427 CPU-h, 42.6 % of the bill.

   Batching 302,600 calls into ~30,000 leaves ~42 CPU-h of overhead, so 427 + 574 → 42 + 574 ≈ 616 CPU-h, a 1.63× saving worth ~385 CPU-h per effect size.

   This refutes the previous revision of this point, which used a three-target profiled fit (3.66s + 0.512s × pairs) to argue the intercept was mostly amortizable and the ceiling only ~1.5×. The real intercept is 5.076s, close to the original 36-target estimate of 4.91s. That three-target fit was thin and this document said so; it was still wrong to revise the headline number on it. The slope was the part the smoke test got wrong: 0.5923s against 0.459s, 29 % low.

2. Memory. A 194-gene batch is 194 × 586,309 dense doubles ≈ 0.9 GB for the count matrix and as much again for the effect-size matrix, against ~54 MB per target today. The 8 GB request may hold but has not been tested.
3. Linearity. The 0.459s slope was fitted over 5–36 pairs per call. Batches carry up to 194, and sceptre’s per-call cost has not been shown to stay linear that far out.

Shape of the change. A batching function (first-fit over gene sets, ordered deterministically so layout is reproducible), the target loop in run_power_simulation.R iterating batches instead of single targets, and split_pairs.R splitting on batches, not targets. Seeds must stay derived from (seed, target, rep, effect_size) so results remain invariant to the batch layout, exactly as they are invariant to the split layout today — that invariance is also what makes the two designs directly comparable pair by pair.

Step 10 — re-derive the environment — DONE

Audited 2026-08-14. The list was in better shape than expected: one package was in the wrong environment, and nothing else needed to change.

Our own code uses only Matrix, sceptre and optparse via library(), plus stats, utils, methods, tools and testthat via ::. Everything else in pixi.toml is there for sceptre’s benefit, which is where the real question was.

r-bh moved to the build environment. sceptre declares LinkingTo: Rcpp, BH — header-only Boost, used to compile the C++ and never loaded at runtime. It had been sitting in the runtime environment that every Nextflow task activates, under a comment describing it as a “NAMESPACE import”, which it is not: sceptre’s NAMESPACE imports only Matrix and Rcpp. That wrong comment is presumably how it got there.

src/audit_dependencies.R is the tool that settled it, and it separates three things that get conflated:

BH was the only declared dependency with no reference anywhere in the package’s code. Everything else earns its place, so the rest of the runtime list stands as documented in pixi.toml.

What the audit cannot tell you is which of those 15 :: references lie on our path specifically — that needs coverage instrumentation. It narrows the question to a list short enough to reason about, and the synthetic object now makes a candidate removal cheap to test by running the whole pipeline without the package.

The re-solve was a no-op apart from the intended change (job 39034832): four lines removed from pixi.lock, no package added, removed or moved version. r-base stayed at 4.4.3-h35b0bb1_11. Verified afterwards by job 39037197, which rebuilt the synthetic object and ran all six pipeline steps against the new environment — if BH were needed at runtime, sceptre would not load at all.

A caveat that outlived this step. Nearly every spec in [dependencies] is "*", so the lockfile is the only thing pinning versions and any re-solve is a gamble on what conda-forge has published since. This one happened to come back clean; the next one may not. Pinning the specs to their locked versions would make future re-solves honest, and is worth doing before the next dependency change rather than during one.

Step 11 — fit the power curve and run three effect sizes instead of six

Prototyped in src/fit_power_curve.R, validated on one dataset, not yet adopted. Halves the cost of a sweep (~2,575 CPU-hours at 100 replicates on this dataset) and makes the minimum detectable effect size continuous rather than snapped to whichever effect sizes were run.

The model is derived, not curve-fitted. For a Wald-type test at a fixed threshold,

power(effect_size) = Phi(beta / SE - z),   beta = -log(1 - effect_size),  z = qnorm(1 - alpha)

beta is the effect on the scale the test works on, SE collects everything pair-specific (perturbed cells, expression, dispersion), and z is fixed by the discovery threshold and shared by every pair. On the probit scale that is a straight line through -z with slope 1/SE: one free parameter per pair, so three effect sizes leave two degrees of freedom to check the fit rather than just enough to force it. src/fit_power_curve.R fits it as a binomial GLM with a probit link, no intercept and -z as an offset — a GLM rather than least squares on qnorm(power) because 0/100 and 100/100 still carry information about the slope, whereas qnorm() would be infinite there.

Validation. Held-out test on dc_tap_paper_wtc11_no_shuf, which has a complete six-point sweep (0.05 / 0.1 / 0.15 / 0.2 / 0.25 / 0.5) over 6,574 pairs at 100 replicates. Fit on three of them, predict the other three, compare against what was measured:

“Noise floor” is the mean binomial standard error of the measured value being compared against, so the best grid predicts held-out power to within about 1.6× the noise of simply measuring it, and reproduces the six-point sweep’s minimum detectable effect size exactly for 89 % of pairs.

Two supporting checks on the same data. The functional form holds: fitting per pair on all six points gives a probit-scale residual sd of 0.132 (median), against ~0.13 expected from 100-replicate binomial noise alone at power 0.5 — the straight line fits as well as the data can distinguish. And monotonicity holds: of 32,870 consecutive-effect-size comparisons only 220 (0.67 %) decrease, largest decrease 0.060, which is about one standard error of a difference.

Grid placement is dataset-specific and matters more than the number of points. wtc11 is well powered and transitions between 5 % and 15 %; day0_grna20_no_shuffle is not, and 84 % of its pairs transition between 5 % and 25 %. Projecting each day0 pair’s curve from its measured 0.15 power:

A one-parameter sigmoid is only well determined by points away from 0 and 1, so a grid that brackets the transition rather than sampling it wastes replicates. Pick the grid from a cheap pilot — one effect size at 30 replicates locates the transition distribution — rather than reusing another dataset’s.

Open points before adopting it.

• Validated on one dataset, and this is now the critical gap. dc_tap_paper_wtc11_no_shuf is the only dataset with a full six-point sweep. Inventory of what exists today:

  Dataset	Effect sizes with power output
  dc_tap_paper_wtc11_no_shuf	0.05 / 0.1 / 0.15 / 0.2 / 0.25 / 0.5
  day0_grna20	0.1 / 0.15 / 0.2
  dc_tap_paper_k562 / dc_tap_paper_wtc11	0.05 / 0.1 (and both carry the size-factor shuffle)
  day0_grna20_no_shuffle, dc_tap_paper_k562_no_shuf, day2, day4	0.15 only

  So the held-out test cannot be repeated on a second dataset without new compute. Sweeps are needed on at least two more datasets, and the two worth doing are dc_tap_paper_k562_no_shuf — same protocol as wtc11, different cell type, so it isolates cell type from method — and Gasperini et al., which is a different lab, protocol and scale entirely and is therefore the real test of generalization. The DC-TAP data is Ray et al.; wtc11 and k562 are two of its cell types, so validating only within it is close to validating within one experiment.

  Note the shuffled variants are not usable as references for absolute power, since the size-factor permutation paired each cell’s library size with another cell’s perturbation status. They may still be usable for checking the functional form of the curve, which is a claim about shape rather than level — worth deciding deliberately rather than by accident.

• z should come from the threshold, but did not match. Fitting z per pair on wtc11 gives a median of 2.045, and the script’s --threshold-file route would use qnorm(1 - threshold). Where those disagree, --fit-z profiles a single shared z by total deviance. Worth understanding why they differ — sceptre’s test is a conditional-resampling test with a skew-normal approximation, not a Wald test, so some deviation is expected.
• Keep the raw per-effect-size power values. The fitted curve is monotone by construction, so the monotonicity check above is only meaningful on unfitted numbers.
• Low-count genes may plateau below power 1, because the resampling p-value has a granularity floor. A two-parameter version (gamma * Phi(...)) would cover that, but then three points leave no slack for checking. deviance / df per pair, which the script reports, is the diagnostic.
• The deviance says the model is imperfect, even though it predicts well. Running the prototype on wtc11 (fit 0.05 / 0.1 / 0.25, predict the rest) reproduces the held-out accuracy — MAE 0.0175 with --fit-z, 0.0156 with z fixed at 2.045 — but reports deviance / df of 2.4–2.7 (median, 90th percentile 8.2). The earlier probit-residual check missed this because it discarded saturated points; the GLM includes them, and binomial deviance is very sensitive there. So the straight line is not the true curve at the saturated ends, while still interpolating the transition to within about 1.5× the Monte-Carlo noise. Use it for interpolation, not for extrapolation past the fitted range, and treat a high deviance / df as “check this pair” rather than as a verdict on the pair.
• Profiling z did slightly worse than fixing it (0.0175 against 0.0156), so --fit-z is a diagnostic for whether the threshold and the curves agree, not the recommended default.
• The MDES columns in the table above treat the measured six-point grid as truth. It is not: each point carries ±0.05, so pairs near the 0.8 boundary flip grid steps easily. Some of the 11 % disagreement is the measurement being wrong rather than the fit, which pools 300 draws.

Predicting the curve from covariates instead of measuring it

The obvious extension is to skip simulation altogether: the theory says SE^2 ~ (1 / n_pert_cells) * (1 / mu + 1 / theta), so k = 1 / SE should be predictable from the perturbed-cell count and the gene’s expression, both of which the pipeline already reports per pair. Regressing the fitted k on those two, over the 6,030 wtc11 pairs with a usable fit:

log(k) = -0.605 + 0.502 * log(perturbed cells) + 0.351 * log(expression)
                       theory: 0.500              theory: 0 to 0.5
R^2 = 0.867      residual sd of log k = 0.243  ->  k predicted to within x1.28

The perturbed-cell exponent lands on 0.502 against a theoretical 0.5, so the sqrt(n) law is confirmed on real data. It is still not a substitute for measurement where per-pair claims are concerned. A x1.28 error in k is about ±0.20 in power near the middle of a pair’s transition, against 0.017 for the per-pair curve fit above — and that error is model error, not sampling noise, so it does not shrink with more simulation and is likely concentrated in particular genes rather spread evenly. Certifying a negative on a prediction would mis-certify a non-random subset of pairs.

Where it is the right tool:

• Design counterfactuals — the power gained from twice the cells or more gRNAs per element. Measurement cannot answer these at all.
• Aggregate statements, where ±0.2 per pair averages down over thousands of pairs.
• Choosing the effect-size grid without a pilot run, and deciding which pairs are worth simulating.
• Pairs never tested — the 1,564 that failed pairwise QC on day0, or a planned experiment.

Note one trap in the numbers above: perturbed-cell count alone explains only 2.3 % of the variance while expression alone explains 83.7 %. That is not evidence that cell count is unimportant — its exponent is confirmed — but that it barely varies across pairs in this dataset. Causally important, empirically near-constant.

The defensible per-pair version, if wanted, is to calibrate the residual quantiles on a measured subset and widen each pair’s effect-size interval by them. That gives honest intervals, merely wider than measured ones — and quantifies what a measured sweep is buying.

What this needs before it can carry a paper

This subsection was written as a caveat on a side result. It is no longer a side result: predicting power for pairs that were never simulated is the intended headline, because it is the only route to power for the trans pairs a bootstrap cannot reach (an element against genes on another chromosome — sceptre supports the test, but simulating it pair by pair is infeasible). The claim therefore rests on the least-validated result in this document, and these are the gaps, recorded so they are not discovered late:

• Fitted on one dataset. Cross-dataset transfer is the claim, and it is untested. Needs dc_tap_paper_k562_no_shuf (same protocol, different cell type) and Gasperini et al. (different lab, protocol and scale) — the same two sweeps the power-curve section above needs, so one round of compute serves both. Refit log(k) ~ log(pert_cells) + log(expression) per dataset and ask whether the coefficients transfer, not just whether each dataset fits well on its own. A model that has to be refitted per dataset is still useful but is a different, weaker claim.
• A ×1.28 error in k is ≈ ±0.20 in power mid-transition, and it is model error, not sampling noise, so it does not shrink with more replicates. It is also unlikely to be spread evenly across genes. For training labels in scE2G / ENCODE-rE2G, that means mislabelled negatives concentrated in a non-random subset of pairs — which is worse for a classifier than uniform noise. Quantify which pairs the residuals concentrate in (expression decile? dispersion? cell count?) rather than reporting the residual sd alone.
• The covariate model is close to a one-covariate model on this data. Expression alone explains 83.7 % of the variance in k, perturbed cells alone 2.3 %. The sqrt(n) exponent is confirmed at 0.502 against a theoretical 0.500, so cell count is causally right and empirically near-constant here — but a reviewer will ask, and a dataset with real variation in perturbed-cell count is what answers it. Gasperini et al. differs enough in scale to provide that.
• The trans case needs a demonstration, not an argument. Simulate a tractable subset of trans pairs directly, then compare against what the covariate model predicts for them. Without that, the headline is an extrapolation from cis pairs to a regime where nothing has been measured — and trans pairs may differ systematically, since the gene sets are not matched to the element’s neighbourhood.
• Calibrated intervals should be a contribution, not a footnote. Prediction plus honest residual-calibrated intervals is a defensible per-pair claim; a point prediction with ±0.20 is not.

Prior art to check before writing any of this up: scPower (Schmid et al.) is the work reviewers will name, and the claim that no satisfying method exists needs an actual literature search rather than an assumption.

Running the comparisons on the cluster

Setup

git clone https://github.com/EngreitzLab/element-gene-power-analysis.git
cd element-gene-power-analysis
pixi install
pixi run setup        # compiles sceptre; needs a compiler and network
pixi run check-api    # must print "All checks passed"

pixi run setup installs sceptre into .pixi/rlibs inside the checkout, so the whole environment is self-contained and needs no module system. It does require outbound network access to GitHub — if the compute nodes are offline, run setup on a login node first.

Smoke test before committing CPU-days

Three targets, 2 replicates — should finish in a couple of minutes and proves the environment works.

mkdir -p prepared splits sim
Rscript src/prepare_sim_input.R \
  --sceptre-object results/sample1/sceptre_object.rds --outdir prepared/

Rscript src/split_pairs.R --pairs prepared/pairs.tsv --n-splits 280 --outdir splits/

Rscript src/run_power_simulation.R \
  --sim-input prepared/sim_input.rds --sceptre-template prepared/sceptre_template.rds \
  --pairs splits/split_001.tsv --grna-targets prepared/grna_targets.tsv \
  --effect-size 0.15 --reps 2 --seed 20250812 --out sim/smoke.tsv

Full run as a SLURM array

prepare_sim_input.R once, then one array task per (split, effect size).

These snippets are illustrative and predate workflow/slurm_executor/, which is what actually runs — prefer it. In particular they omit --null-precomputations, so as written they would run the cleared configuration at ~4.3× the cost: the template’s inherited slot is now emptied by slim_sceptre_object(), so sceptre refits every gene’s null model inside every call. Add steps 02b and 02c and pass the bundle.

#!/usr/bin/env bash
#SBATCH --job-name=power-prepare
#SBATCH --mem=12G
#SBATCH --time=1:00:00
#SBATCH --cpus-per-task=1
pixi run Rscript src/prepare_sim_input.R \
    --sceptre-object results/sample1/sceptre_object.rds \
    --outdir prepared/
pixi run Rscript src/split_pairs.R \
    --pairs prepared/pairs.tsv --n-splits 280 --outdir splits/

#!/usr/bin/env bash
#SBATCH --job-name=power-sim
#SBATCH --array=1-280
#SBATCH --mem=4G
#SBATCH --time=3:00:00
#SBATCH --cpus-per-task=1
#SBATCH --output=logs/sim_%A_%a.out

# Selected by position rather than by formatting the index: split_pairs.R pads the filename to
# the width of --n-splits, so 280 splits give split_001.tsv but 1000 would give split_0001.tsv.
SPLIT=$(ls splits/split_*.tsv | sed -n "${SLURM_ARRAY_TASK_ID}p")
for ES in 0.15 0.2; do
    pixi run Rscript src/run_power_simulation.R \
        --sim-input prepared/sim_input.rds \
        --sceptre-template prepared/sceptre_template.rds \
        --pairs "$SPLIT" \
        --grna-targets prepared/grna_targets.tsv \
        --effect-size "$ES" --reps 100 --seed 20250812 \
        --out "sim/split_${SLURM_ARRAY_TASK_ID}_es${ES}.tsv"
done

Then per effect size:

for ES in 0.15 0.2; do
    pixi run Rscript src/compute_power.R \
        --simulations "$(ls sim/*_es${ES}.tsv | paste -sd, -)" \
        --threshold-file prepared/discovery_threshold.txt \
        --out "power_es${ES}.tsv"
done

pixi run Rscript src/summarize_power.R \
    --power power_es0.15.tsv,power_es0.2.tsv --out power_summary.tsv

Expect roughly an hour per array task at 10 targets × 100 replicates × 2 effect sizes. Time limit is set to 3h for headroom — per-target cost varies about 1.5× with pair count.

Comparison 1 — old vs new pipeline — DONE, PASSED (job 38916341, 2026-08-13)

The refactor is validated. The old pipeline’s 2026-05-14 output was still on disk (results/day0_grna20_no_shuffle/power_analysis/power_analysis_results_es_0.15.tsv, 34,886 pairs), so the old environment never had to be rebuilt — which retires the CentOS 7 / pinned-linux-64-solve risk described below. old_vs_as_is is the gating comparison, because the old pipeline inherits the same real-data cache as_is does; comparing it against null_fit would fold in a deliberate change of null model.

Check 2 is the strongest of these and is easy to skim past. Perturbed cells per pair is RNG-independent, so a zero difference across all 34,886 pairs proves the two implementations select identical cells for identical pairs — the entire upstream path is bit-identical, and only the simulation’s RNG stream differs. Checks 3–5 are then about Monte-Carlo agreement alone.

Check 5’s mean absolute difference of 0.0304 is not a discrepancy: at 100 replicates a single power estimate carries a standard error up to 0.05, and both sides carry one, so per-pair scatter of this size is what agreement looks like. What would have been alarming is a directional offset, and check 4 finds none — 13,823 against 13,642 is as close to a coin flip as this many pairs allow.

The old_vs_null_fit comparison “fails” check 3, and that is the correct behaviour. Mean power differs by 0.0066 = 28.1 σ. That is not a refactor problem — it is the deliberate null-model change, re-measured independently: null_fit − as_is was +0.0063 measured directly, null_fit − old is +0.0066, and as_is − old is +0.0003. The three are additive and consistent, which is a useful triangulation. Check 3 tests “is this the same estimator?” against the Monte-Carlo SE of a mean over 34,886 pairs (0.00024), so any real systematic shift blows through it however small. It is reported and does not gate.

Reports: results/refactor/comparison/old_vs_as_is/ and .../old_vs_null_fit/.

Why it was outstanding: the upstream statistics are proven bit-identical (size factors, normalised means, raw means: max absolute difference exactly 0), and dropping @grna_matrix was proven to leave discovery results unchanged. What has not been checked is the simulation itself, because seeding moved inside the replicate loop, which changes the RNG stream. So the comparison has to be distributional rather than exact.

The old code needs the old environment, which is why the cluster is the right place: the old envs/power_analysis.yml is a fully-pinned linux-64 solve. It cannot be solved on macOS, but should solve on a Linux compute node.

# 1. check out the pre-refactor tree alongside the current one.
#    b3cf7ab, not c7448a7: it is the tip of the old `low-mem` branch and includes the
#    "Improve memory usage in power simulations" commit, so it is the version actually
#    being run in practice.
#
#    Sherlock's default git is 1.8.3.1, which has no `worktree` subcommand -- load a
#    modern one first. Do NOT `git checkout` the old commit in the working tree instead:
#    the sbatch scripts run `src/` from the repository root, so swapping branches under a
#    running array silently changes the code mid-run.
ml system git/2.45.1
git worktree add ../egpa-old b3cf7ab
cd ../egpa-old

# 2. the old pipeline expects four inputs under results/<sample>/ with these exact names
mkdir -p results/sample1
cp /path/to/sceptre_object.rds        results/sample1/final_sceptre_object.rds
cp /path/to/grna_groups_table.rds     results/sample1/
cp /path/to/gene_grna_group_pairs.rds results/sample1/

# the fourth input has to be an .rds, so convert the discovery-results CSV:
Rscript -e 'saveRDS(read.csv("/path/to/sceptre_discovery_results.csv"),
                    "results/sample1/results_run_discovery_analysis.rds")'

# 3. run the old pipeline for the power results, NOT `all`
#    (`all` stops at the concatenated per-replicate file)
snakemake --use-conda --conda-frontend conda --cores 8 \
    results/sample1/power_analysis/power_analysis_results_es_0.15.tsv

Set num_replicates in that tree’s config/config.yml to whatever you use for the new run so the two are comparable, and remember the old config/config.yml also has the dead reps: 20 key — it is num_replicates that takes effect.

Then compare against the new power_es0.15.tsv. Suggested acceptance criteria:

A systematic offset in either direction is a real discrepancy and worth chasing. Scatter across individual pairs is expected — at 100 replicates the standard error of a single estimate is up to 0.05, so differences of ±0.1 on individual pairs are normal.

Use the same --reps in both, and note the old pipeline cannot be seeded, so it produces a different draw every run. Running it two or three times gives a sense of its own run-to-run spread, which is the yardstick for judging the old-vs-new difference.

Comparison 2 — monotonicity at full replicate count

Power must not decrease as effect size increases. At 12 replicates on 33 pairs this held (mean power 0.356 → 0.604 → 0.838 for 0.15 / 0.25 / 0.5) with one violation of 0.08 that a two-proportion test could not distinguish from noise. Worth repeating at 100 replicates across all 32,386 pairs, where noise-driven violations should nearly vanish. A violation that survives is a bug.

Comparison 3 — timing and memory at scale

Run with --with-report-style accounting (or just sacct) and record actual per-task wall time and MaxRSS, so the resource requests above can be tightened. The 4 GB figure comes from a laptop run of three targets; confirm it holds for splits containing the largest targets (up to 36 pairs).

Decisions already made

Recorded so they are not relitigated:

The analysis question this is all for

Recorded because it determines which columns matter, and therefore how many replicates to buy.

The goal is false-negative triage: for a pair CRISPRi did not call, is the link absent, or could the experiment not have seen one? Three consequences follow, and Output explains each in full.

1. Threshold power_ci_low, never power. “This negative is biological” asserts that power was at least 0.8. Certifying that needs 88/100 successes — or 29/30, which is why 30 replicates is unusable for this question regardless of what the precision tables say. At 0.15 and 100 replicates, under null_fit (power_null_fit/power_es0.15.tsv): 38.3 % of pairs certified, 12.6 % ambiguous, 49.0 % clearly underpowered (that last half is not a replicate problem — it needs more perturbed cells or a larger effect). The as_is run gave 37.5 / 12.6 / 49.9.

2. Minimum detectable effect size is the deliverable, not six power columns. summarize_power.R now emits min_detectable_effect_size plus a _ci_low / _ci_high bracket, where _ci_high is derived from power_ci_low: the interval inverts, as power rises with effect size. It turns the sweep into one sentence per pair: “we would have caught a knockdown of ≥ 25 %, so the absence of a call rules out effects that large.” Qualifying requires clearing the threshold at that effect size and every larger one tested; taking the first effect size that clears lets noise bias every pair towards looking more detectable than it is.

3. Per-element power needs reframing before it can be reported. The tempting version — “this element is well powered whatever gene you pair it with” — is not what the data say. At 0.15 only 20.9 % of the variance in per-pair power is between elements; 79 % is gene to gene within an element, the mean within-element SD is 0.34, and element mean power correlates with mean_pert_cells at only 0.38. A mean over the pairs an element happens to have tested is also incomparable between elements: it inherits the expression levels of whichever genes sit nearby.

The defensible framing is conditional: per-pair power answers “could we have detected this link?”, per-element power answers “how well did we perturb this element?” — whose sufficient statistic is the perturbed-cell count, not anything about genes. Report it at a reference gene: “for a gene at median expression, element E has 0.7 power at a 15 % knockdown.”

Not implemented. The fit is cheap and needs no new simulation — power_summary.tsv already carries mean_pert_cells and average_expression_all_cells per pair, so fitting power ~ f(pert_cells, expression) and evaluating at reference expression is post-processing. Simulating a reference-gene panel instead would cost ~450 CPU-h, because the per-target term is paid whether a target carries one gene or thirty. Whoever picks this up should decide where it lives: a new src/ script, or a section of summarize_power.R.

What does hold at element level is a floor: 21.9 % of elements (663/3,026) have no tested pair that could have reached power 0.8, so they cannot support a “regulates nothing” claim under any reading, and they should be excluded from biological interpretation rather than reported as null results. Conversely only 5 of 3,026 elements have every tested pair certified — element-wide negative claims are essentially never assertable at 100 replicates and a 15 % knockdown.

Open questions

• Two-stage replicate allocation: worth wiring in? Costed properly it is ~1,768 CPU-h per effect size against 3,430 for a uniform 400 and 858 for a uniform 100 — a 1.9× saving over uniform 400, not the 3.7× an earlier version of this document claimed. That claim came from counting replicate equivalents, which ignores the per-target term; see Choosing num_replicates. It only pays off if per-pair cutoff decisions are being made — which, per the section above, they are.
• Is the 5 % effect size worth 100 replicates? Nearly every pair is underpowered there, so the replicates buy a precise estimate of a number that is 0. Sweeping 5 % at 30 reps and spending the difference on 20–25 % would carry more information. Untested.

• --target-overhead in split_pairs.R: closed, and the answer is “leave it”. It defaults to 0, so every split is balanced on pair count alone. The as_is array fixed the right value at the time (5.076s per target against 0.5923s per pair ⇒ ≈ 8.6 pair-equivalents), but setting it bought almost nothing at N_SPLITS=1000, and the reason this document previously gave for setting it was wrong. The gRNA reuse has since made it moot: at 1.140s per target against 0.5561s per pair the right value is ≈ 2.0 pair-equivalents, on splits that already hold 34–36 pairs. There is nothing left to correct for.

  Measured on the real splits, and on job 38849611’s own sacct record:

  	Predicted tail task	Predicted floor	max/mean
  --target-overhead 0 (what ran)	68.4 min	43.0 min	1.139×
  --target-overhead 8.6	67.4 min	58.0 min	1.123×

  Simulating the same LPT bin packing split_pairs.R uses reproduces the on-disk splits exactly at overhead 0, so the comparison is like for like. The overhead raises the floor a lot (43 → 58 min) and moves the tail by one minute — and the tail is what determines when the array finishes. At 1,000 splits each task holds only 1–4 targets, so there is little imbalance left to reclaim.

  This refutes the previous revision of this point, which blamed the observed 35–110 minute spread on pair-only balancing. It is not that: the cost model predicts a 43–68 minute range, and regressing actual task time on predicted gives r = 0.373, r² = 0.139 — split composition explains 14 % of the variance. The model is unbiased on the mean (36.02s predicted against 35.97s actual per replicate), so it is not simply wrong. What explains the spread is node heterogeneity: across the 229 nodes the array landed on, the slowest node’s mean is 2.73× the fastest (55.8s against 20.5s per replicate). There was no preemption to blame (999 COMPLETED, 1 FAILED). Set --target-overhead 8.6 if N_SPLITS ever drops far enough that splits hold ten or more targets; do not re-split an array for it.

• ODM / out-of-core support is designed for but not implemented. The seam is get_response_matrix() in lib/sceptre_io.R, which currently errors explicitly on an odm. The hard part is that poscounts size factors need a per-cell median over genes — a column-wise reduction — while an odm is row-accessible.