Usage

nextflow run . -profile sherlock -params-file config/config.yml

That is the whole pipeline. Every step is also a console script taking explicit arguments and printing --help, so any one of them can be run, re-run or debugged on its own; nothing reads a global config object.

Getting a dataset in

The pipeline’s input is a .h5mu, not a sceptre object. Producing one is a one-off step per dataset, and it is the only place R appears:

Rscript pysceptre/scripts/export_sceptre_dataset.R \
    --sceptre-object results/sample1/sceptre_object.rds \
    --out-dir export/ --all-genes --all-cells
python  pysceptre/scripts/make_h5mu.py export/

The object must have been through assign_grnas() and run_qc() with grna_integration_strategy = "union".

Both flags are required, and for different reasons.

--all-genes because the poscounts size factors are a per-cell reduction over the whole gene set.

--all-cells because they are also computed against a per-gene geometric mean taken over every cell in the object, including the ones QC removed. An export restricted to the QC-passing cells gives different size factors for the cells that remain, and different normalised means for every gene — measured on day0, a 0.46 % shift in size factor and 2.2 % in the raw gene mean. Nothing downstream would report it. watteg-prepare-sim-input warns if it is handed an export without them.

The export handles odm-backed (out-of-core) matrices itself, which is why nothing in the pipeline takes an --response-odm flag any more.

The five steps

dataset.h5mu
      |
      |  watteg-prepare-sim-input      once per sample
      +--> sim_input.h5                per-gene and per-cell statistics, the model
      |                                coefficients, and the perturbation matrices
      +--> pairs.tsv                   QC-passing discovery pairs
      +--> pairs_with_info.tsv         every pair with its real-data QC counts
      +--> grna_targets.tsv            gRNA -> target mapping (many-to-many)
      +--> discovery_threshold.txt     the p-value a simulation must beat
      +--> analysis_mode.tsv           which test the screen ran, and its resampling budget
      |
      |  watteg-split-pairs            once per sample
      +--> split_001.tsv ...           balanced chunks of targets
      |
      |  watteg-fit-null-models        once per sample, under --null-fits reuse
      +--> null_fits.h5                each gene's null-model fit, per simulation
      |
      |  watteg-run-power-simulation   once per (split, effect size, simulation chunk)
      +--> *.partial.tsv.gz            per-pair counts (--partials-out), by default
      +--> *.tsv.gz                    one row per (pair, simulation) (--out), optional
      |
      |  watteg-consolidate-replicates once per effect size, only if the rows were written
      +--> replicates_es*.parquet      the same rows, as one file
      |
      |  watteg-compute-power          once per effect size, from the counts or the rows
      +--> power_es*.tsv               one row per pair, with confidence intervals
      |
      |  watteg-summarize-power        once per sample
      +--> power_summary.tsv           one row per pair, one column per effect size

There is no sceptre_template.rds. It was an R object carrying the covariate matrix and the analysis parameters, which now live in sim_input.h5 and analysis_mode.tsv.

watteg-fit-null-models is the Python version of R’s FIT_NULL_MODELS: each gene’s null model is fitted once per simulation, on an independent draw with no knockdown, and every target the gene is tested with reuses that fit (--null-fits reuse). Without it (--null-fits refit) each simulation refits every gene once per target, the exact configuration, at about 1.7× the cost of a simulated pair-test. See Methods.

1. watteg-prepare-sim-input

watteg-prepare-sim-input --dataset export/dataset.h5mu --outdir prepared/

Fits each gene’s Poisson GLM and negative-binomial theta — the same model sceptre’s own perform_response_precomputation fits, reproduced to about ten significant digits — and stores the coefficients. A gene whose theta hits the estimator’s bounds is refused rather than simulated from: a clamped theta is not an estimate, and drawing counts from it would state a noise level the data never supported. (sceptre clamps and carries on, which is reasonable for an analysis and not for a simulation.)

2. watteg-split-pairs

watteg-split-pairs --pairs prepared/pairs.tsv --n-splits 280 --outdir splits/

A target’s pairs cannot be separated — the simulation draws one count matrix per (target, simulation) and tests every one of that target’s genes against it — so targets are the unit and the job is bin packing. Weighted by pairs plus an overhead, because a task’s cost is intercept + slope × pairs and weighting by pairs alone over-fills the splits that hold many small targets.

3. watteg-run-power-simulation

watteg-run-power-simulation \
    --prepared prepared/ --pairs splits/split_001.tsv \
    --effect-size 0.15 --reps 100 --seed 20250812 \
    --null-fits-file prepared/null_fits.h5 \
    --out sim/split_001_es0.15.tsv

The defaults are the fast configuration, each part measured before it was adopted (Methods). The previous configuration is --driver engine --null-fits refit, plus --permutations per-replicate --nulls scan for the exact earlier behaviour. A CRT screen needs --driver engine --null-fits refit: the fast driver runs the permutation test only, and says so rather than testing a CRT screen by permutations.

--effect-size 0 is allowed on purpose. It simulates no effect: perturbed cells are drawn from the unperturbed mean and the screen’s own test is run on them, so the p-values should be uniform. It is the only way to measure this pipeline’s type-I error from its own output. Run it as its own sweep, not as another point on a curve. The two estimands are the same simulation there.

--expression-model. fitted draws from exp(X·β), the expected count sceptre’s own null model gives that cell, so the simulation and the test that judges it are on one scale. size_factor reproduces the pre-2026-09-21 behaviour — a size-factor-normalised gene mean scaled by the cell’s poscounts factor — and exists only to compare against sweeps already run with it. It mixes two models, runs about 4 % low, and reproduces 86.5 % of the observed count variance against the fitted model’s 99.5 %.

3b. watteg-fit-null-models

watteg-fit-null-models --prepared prepared/ --reps 100 --seed 20250812 --n-jobs 8 \
    --out prepared/null_fits.h5

Pass the file to watteg-run-power-simulation --null-fits reuse --null-fits-file. A task given no file fits its own genes first with the same keyed draws, so the output is the same bytes either way; the file only saves refitting a gene in every task that tests it. Refused for a CRT screen, whose simulations the fast driver (the file’s only reader) does not run.

4. watteg-compute-power and watteg-summarize-power

watteg-compute-power \
    --simulations sim/ --threshold-file prepared/discovery_threshold.txt \
    --out power_es0.15.tsv

watteg-summarize-power \
    --power power_es0.05.tsv power_es0.15.tsv \
    --sim-input prepared/sim_input.h5 --out power_summary.tsv

watteg-compute-power takes either --partials (the simulation tasks’ per-pair counts) or --simulations (one row per pair and simulation), and gives the same table from both when each pair’s simulations came from one task. Both take files or directories; a directory expands to the .tsv, .tsv.gz and .parquet files inside it, sorted. In the pipeline the rows are written only under keep_per_simulation: true, or when reps_per_chunk < num_replicates. --power-threshold (default 0.8) sets what “detectable” means for the minimum-detectable-effect-size columns.

Reproducibility

A pair’s power does not depend on how the work was split up. Every draw is keyed on (seed, target, simulation, effect size) and nothing else, so running simulations 1–100 in one task and in five tasks of twenty gives byte-identical results, and so does moving a target from one split to another. A reused null-model fit is keyed on (seed, gene, simulation), so it is the same wherever it was made. That is checked end to end:

WATTEG_PREPARED=prepared/ pytest -m realdata

The streams differ from the R implementation’s and always would — different generators — so the two agree statistically rather than draw for draw. See Plan - pysceptre backend for what has been measured.

Sizing a cluster run

The unit of work is (split × effect size × simulation chunk), and essentially all of the compute is in step 3. A task’s cost follows intercept + slope × pairs, so:

• halving the pairs saves less than half, because the intercept does not move;
• halving the simulations saves exactly half;
• effect sizes cost the same as each other.

--n-jobs buys wall clock rather than CPU time: measured on a 10-performance-core machine, 8 workers gave 2.5× the wall-clock speed for 1.5× the CPU, and 14 was slower than 8.

A task’s wall time is bounded by its largest target, whose simulations are one unit of work. On one machine, several narrower tasks at once therefore beat one wide one: on a 14-core laptop, three 330-pair cis tasks at 4 workers each finished in 727 s, against 296 s each at 8 workers one after another. For a local sweep, -profile local --power_simulation_cpus 4 with Nextflow’s executor.cpus = 12 runs three at once; each took ~4.2 GB on macOS.