Output

power_summary.tsv — the main result

One row per element–gene pair, one set of columns per effect size.

NA in any of the three is a statement about the effect sizes you ran, not evidence that a pair is undetectable. If you tested 0.15 and 0.2 and a pair needs 0.4, it will be NA. That is why max_effect_size_tested sits next to it.

Two things about how these are derived are easy to get wrong, so they are stated explicitly.

The direction of the interval inverts. Power rises with effect size, so a lower bound on power gives a larger minimum detectable effect size. min_detectable_effect_size_ci_high is therefore computed from power_ci_low, not from power_ci_high. It is the pessimistic reading — the smallest knockdown the data can establish you would have caught — which is why it, not the point estimate, is what a false-negative argument rests on. The bracket is built by thresholding each effect size’s own 95 % interval, so read it as a conservative-to-optimistic range rather than an exact 95 % interval for the effect size itself.

A pair must clear the threshold at its effect size and at every larger one tested. Reporting the first effect size that clears, in isolation, lets Monte-Carlo noise decide: at 100 simulations a pair whose true power is 0.75 clears 0.8 about a third of the time, so across six effect sizes a spurious early clear is likely, and the error only ever runs one way — reporting the pair as more detectable than it is. Effect sizes a pair was not tested at count as unknown, not as failures, so they do not block the run.

Why gene_mean and dispersion are here

Power for a pair is governed by

SE^2 ~ (1 / n_pert_cells) * (1 / mu + 1 / theta)

so anything that models power from covariates needs the gene’s dispersion alongside its expression and the perturbed-cell count. Until these columns existed, only average_expression_all_cells reached this table and every such analysis had to load a separate sim_input.h5 to find the rest.

With the convention above, the bracket is 1/gene_mean + dispersion — no conversion, because dispersion already is 1/theta.

They are joined by watteg-summarize-power at summary time rather than emitted per simulation, because they are per-gene constants: putting them in the simulation output would only help sweeps run after the change, while re-summarising an existing sweep from its stored power tables takes seconds. Both columns are absent if --sim-input is not passed, so older outputs and hand-run summaries are unaffected.

Interpreting a negative result

The usual reason to run this pipeline is to decide what a non-significant pair means: did CRISPRi fail to show an effect because there is no regulatory link, or because the experiment could not have seen one? Power is what separates those, but only if it is read the right way.

Per pair — threshold the lower bound

The claim “this negative is biological” is a claim that power was at least 0.8, not that the point estimate landed above it. So threshold power_ci_low (or take min_detectable_effect_size_ci_high, which does it for you across effect sizes). That is a stricter test than it looks, and the strictness is the reason simulation count matters here more than anywhere else:

Measured on 34,886 pairs at effect size 0.15 and 100 simulations: 37.7 % of pairs have measured power ≥ 0.8 (power_ci_low), 12.4 % are ambiguous, and 49.9 % are clearly underpowered (power_ci_high < 0.8). That last half is not a precision problem — no number of simulations fixes it; those pairs needed more perturbed cells or a larger effect.

The sentence the table supports, per pair, is then:

min_detectable_effect_size_ci_high = 0.25 → we would have detected a knockdown of 25 % or more, so the absence of a call rules out effects that large, and says nothing about smaller ones.

Per element — power is mostly a property of the gene, not the element

It is tempting to summarise power per grna_target and say “this element is well powered whatever gene you pair it with”. The data do not support that framing. Decomposing per-pair power at effect size 0.15:

The gene matters about four times more than the element, and a single number per element hides a ±0.34 spread. A mean over the pairs an element happens to have tested is also not comparable between elements: an element surrounded by highly expressed genes scores better than one surrounded by lowly expressed genes at identical perturbation efficiency.

What is well defined is the conditional version. Per-pair power answers “could we have detected this link?”; per-element power answers “how well did we perturb this element?” — a question whose sufficient statistic is the perturbed-cell count, not anything about genes. So report it at a reference gene: “for a gene at median expression, element E has 0.7 power to detect a 15 % knockdown.” That is comparable across elements. It also costs nothing extra: power_summary.tsv carries both mean_pert_cells and average_expression_all_cells per pair, so fit power ~ f(mean_pert_cells, average_expression_all_cells) and evaluate at a reference expression. Simulating a reference-gene panel instead would cost ~450 CPU-hours, because the per-target term in the cost model is paid whether a target carries one gene or thirty — see Choosing num_replicates.

The element-level statement that does hold up is a floor: 21.9 % of elements (663 of 3,026) have no tested pair that could have reached power 0.8 at all, and those elements cannot support a “regulates nothing” claim under any reading. In the other direction, only 5 of 3,026 elements have measured power >= 0.8 for every tested pair, so element-wide negative claims are essentially never assertable at 100 simulations and a 15 % knockdown. Worth knowing before building a figure around them.

power_es<effect_size>.tsv — per effect size

By default this table is built inside the simulation: each task holds all simulations of its pairs and writes per-pair counts (simulations called, simulations used, and the sums behind the means), which COMPUTE_POWER adds up. That is the same table the per-simulation rows give, byte for byte: the means are sums over counts either way, which is exactly how pandas computes a grouped mean, and the rows are read back at full precision. The counts are not published.

What power actually counts

A simulation counts as a success when both hold:

p_value < threshold   AND   log_2_fold_change < 0

The second condition makes this one-sided: a simulation only counts if the simulated perturbation reduced expression. threshold is the largest nominal p-value that survived multiple-testing correction in the real discovery analysis, so “detected” means the same thing here as it did in your actual results.

Reading the confidence intervals

Always read power together with its interval, especially at the boundaries. power = 0 does not mean the pair is undetectable — at 100 simulations the Wilson upper bound is 0.037, so the data are consistent with a true power of a few percent. Likewise power = 1 has a lower bound of 0.963.

The interval is Wilson rather than p̂ ± 1.96·SE precisely because the normal approximation returns [0, 0] for 0 successes, asserting certainty the data do not support. See Choosing num_replicates.

per_replicate/replicates_es<effect_size>.parquet — per-simulation detail (optional)

One row per (pair, simulation). Written only with keep_per_simulation: true, or when simulations are chunked across tasks (reps_per_chunk < num_replicates), where power is computed from these rows. It is the only output from which power can be re-derived — subsampling simulations to study a reduced design, bootstrapping, or re-thresholding — so keep it for a sweep you will want to re-analyse. It dominates the output volume: 74 million rows for moi5 trans.

Five columns sceptre returns are dropped before writing, because at 100 simulations × 34,886 pairs × 6 effect sizes they were 39 % of a 3 GB output and nothing read them:

Together with gzip that takes the per-simulation output from ~3 GB to a few hundred MB for a six-point sweep at 100 simulations. Nothing needs a decompression step: read.delim sniffs the magic number and handles .tsv.gz and .tsv alike, and watteg-compute-power accepts either.

Intermediates

There is no sceptre_template.rds: it was an R object carrying the covariate matrix and the analysis parameters, which are in sim_input.h5 and analysis_mode.tsv. null_fits.h5 plays the part R’s null_precomputations.rds did, made in Python by pysceptre’s own fit (see Methods); under null_fits: refit there is none, and every gene is refitted on each simulation’s own counts.

split_*.tsv is not published. The splits are parallelisation bookkeeping — 1,000 files per sample — and they are regenerable: the bin packing is deterministic given pairs.tsv and n_splits. They live in the Nextflow work directory and go when it is cleaned.

Two means, deliberately

sim_input.h5 carries two per-gene means and they are not interchangeable:

• mean — size-factor normalised. The simulation drew counts from this until 2026-09-21; it now draws from exp(X·β) instead, and mean is kept because --expression-model size_factor reproduces the sweeps that were run on it.
• average_expression_all_cells — raw. Reported in the output so power can be related to expression level.

Only the raw one appears in the output tables. The two differ by more than a rounding: on day0 the normalised mean sits 16 % below the raw one, which is what a normalisation does and is also why feeding the wrong one to anything expecting counts is a scale error rather than a small one.