Usage

Every step is a standalone executable in src/ that takes explicit command-line arguments and prints --help. Nothing reads a global config object, so any step can be run, re-run or debugged on its own.

Orchestration is in progress. The steps below are the complete, working pipeline as run by hand or from a shell loop. A Nextflow workflow that wires them together with a SLURM profile is the next piece of work and is not part of this repository yet. config/config.yml already holds the parameters it will consume.

The five steps

sceptre_object.rds
      |
      |  prepare_sim_input.R          once per sample
      +--> sim_input.rds              per-gene / per-cell statistics + perturbation matrices
      +--> sceptre_template.rds       the sceptre object with both count matrices emptied
      +--> pairs.tsv                  QC-passing discovery pairs
      +--> grna_targets.tsv           gRNA -> target mapping
      +--> discovery_threshold.txt    the p-value a replicate must beat
      |
      |  split_pairs.R                once per sample
      +--> split_001.tsv ...          balanced chunks of targets
      |
      |  run_power_simulation.R       once per (split, effect size, replicate chunk)
      +--> sim_*.tsv                  one row per (pair, replicate)
      |
      |  compute_power.R              once per effect size
      +--> power_es*.tsv              one row per pair, with confidence intervals
      |
      |  summarize_power.R            once per sample
      +--> power_summary.tsv          one row per pair, one column per effect size

1. prepare_sim_input.R

Derives everything the simulation needs from the sceptre object, once, so the parallel tasks read small files instead of a multi-gigabyte object.

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

It logs which categorical covariates are available for --cell-batches, which is worth reading: the column is dataset-specific (e.g. batch_factor, replicate_factor) and is often not called batch.

2. split_pairs.R

Splits targets into balanced chunks. Whole targets stay together, because the simulation builds per-target state once and then loops over replicates.

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

Balancing is on pairs, not on target count, and it reports the achieved imbalance so you can see it. It also refuses to produce empty splits.

3. run_power_simulation.R

The actual simulation. One invocation handles one split × one effect size × one chunk of replicates.

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 100 \
    --seed 20250812 \
    --out sim/split_001_es0.15.tsv

Running the full grid by hand

for es in 0.15 0.2; do
    for split in splits/split_*.tsv; do
        name=$(basename "$split" .tsv)
        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/${name}_es${es}.tsv"
    done
done

4. compute_power.R

Turns per-replicate results into power per pair, with a Wilson interval.

Rscript src/compute_power.R \
    --simulations "$(ls sim/*_es0.15.tsv | paste -sd, -)" \
    --threshold-file prepared/discovery_threshold.txt \
    --out power_es0.15.tsv

Run it once per effect size — it refuses input that mixes effect sizes, and refuses duplicated (target, gene, replicate) rows, which is what an overlapping --rep-offset would produce.

5. summarize_power.R

One row per pair, one power_at_effect_size_<N> column per effect size, plus the smallest tested effect size at which each pair reaches a target power.

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

Reproducibility

--seed is required. Seeds are derived per (seed, target, replicate, effect_size), not per task, which has two consequences worth relying on:

• Results do not depend on how work is divided. The same --seed gives identical numbers whether you use 1 split or 280, and whether replicates run in one chunk or ten. This is verified in the test suite.
• Runs are extensible. Going from 100 to 400 replicates leaves replicates 1–100 byte-identical, so you can add replicates with --rep-offset 100 instead of recomputing.

Sizing a cluster run

Measured on a 292-gene × 586,309-cell dataset with 2,798 targets and 32,386 QC-passing pairs, using all control cells:

Prefer more splits over more replicate chunks. An extra split costs one extra ~1.1s input load; an extra replicate chunk re-pays the per-target setup (control selection, guide status, cell permutation) once per chunk per target. Only chunk replicates when you need more parallelism than one-target-per-task, or when a single task would exceed your queue’s time limit.

For the dataset above, ~10 targets per task is about an hour of work, so --n-splits 280 with no replicate chunking gives 280 tasks per effect size — comfortably under the MaxArraySize of 1000 that most SLURM installations use.