Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
84 commits
Select commit Hold shift + click to select a range
f93fe24
Plan the pysceptre backend for the high-MOI power analysis
emattei Sep 21, 2026
308b40e
Record what phase 1 found
emattei Sep 21, 2026
c530cc2
Call the size of the size-factor question an estimate, not a bound
emattei Sep 21, 2026
9fe30fc
Note that the real-data regression passed too
emattei Sep 21, 2026
deee55f
Point the plan at the squashed commit
emattei Sep 21, 2026
3ddbfc4
Port the expression statistics to Python, and gate them against R
emattei Sep 21, 2026
5d63834
Fit the dispersions in Python instead of reading sceptre's cache
emattei Sep 21, 2026
61152ac
Port prepare_sim_input to Python, and check it against R's own output
emattei Sep 21, 2026
3a9225c
Name the degenerate fits, gate the dispersions, record phase 2
emattei Sep 21, 2026
9b2e9a4
Record phase 2's results in the plan
emattei Sep 21, 2026
8e25a5b
Validate in Python, not across the language boundary
emattei Sep 21, 2026
8d8dfa2
Decide against porting the control-sampling flags, and say why
emattei Sep 21, 2026
db520a4
Record that the simulated mean is biased low, and by how much
emattei Sep 21, 2026
04c7646
Say plainly that sceptre is not involved in the mean bias
emattei Sep 21, 2026
f8ab01b
Connect the mean question to the analytical side
emattei Sep 21, 2026
90308fb
Record what PerturbPlan's post-hoc source actually requires
emattei Sep 21, 2026
8aa55aa
Gate the fitted coefficients, not just theta
emattei Sep 21, 2026
d4fe805
Simulate from sceptre's own model, not from a normalised mean
emattei Sep 21, 2026
5fbb839
Stop claiming the sweeps are conservative
emattei Sep 21, 2026
859df71
Draw from sceptre's own model in R too
emattei Sep 21, 2026
87cfc6a
Carry the per-pair QC counts through to the simulation
emattei Sep 21, 2026
c25f56c
Port the simulation core: seeds, guide assignment, the draw
emattei Sep 21, 2026
ca8ad94
ADD the phase-3 benchmark: two engine shapes, four cost terms
emattei Sep 21, 2026
7f79877
Run the simulation in Python, one call per (target, replicate)
emattei Sep 21, 2026
7bbb8d3
Check the two implementations agree, against noise rather than a tole…
emattei Sep 21, 2026
6c8966f
Port the remaining four pipeline steps
emattei Sep 21, 2026
115c10e
Rewire the pipeline onto the Python backend, and drop R from the envi…
emattei Sep 21, 2026
ed303c3
Check end to end that the chunk layout reaches no number
emattei Sep 21, 2026
b76ec8c
Point the docs at the Python pipeline
emattei Sep 21, 2026
f48a73b
Record where the port stands, and what would be wrong to skip
emattei Sep 21, 2026
3071abb
Generate the synthetic fixture rather than committing one
emattei Sep 21, 2026
24f3da1
Rebuild the container for the Python environment
emattei Sep 21, 2026
c85cc34
Rewrite Usage for the Python pipeline
emattei Sep 21, 2026
54298a5
Describe the model the simulation actually uses
emattei Sep 21, 2026
dd1315f
Bring the docs landing page to the Python pipeline
emattei Sep 21, 2026
6d55f19
Update Output for the Python intermediates
emattei Sep 21, 2026
5ab9d9f
Say that the control-sampling flags are gone, and why that is safe
emattei Sep 21, 2026
c7d3cea
Lint Python in the hook, and test the pipeline in CI
emattei Sep 21, 2026
b6290bb
Give Development the commands that exist
emattei Sep 21, 2026
011d3dc
ADD the Stage B comparison: Python power against a published R sweep
emattei Sep 21, 2026
d10fd5d
Point at the r-implementation branch, not legacy
emattei Sep 21, 2026
c2faa4f
Write down how to run Stage B, and what each replicate count buys
emattei Sep 21, 2026
09105ac
Measure the benchmark's stage breakdown outside the shapes' timing
emattei Sep 21, 2026
2c976ce
Record what re-running the sweeps at the new scale would cost
emattei Sep 21, 2026
5aeb3f3
Centre the effect-size matrix after reordering it, not before
emattei Sep 21, 2026
c1dde55
Stage B needs a regenerated R reference, not the old sweep
emattei Sep 21, 2026
e4469a3
Stop calling the old sweeps published
emattei Sep 21, 2026
bbc8c95
Lock pysceptre at the v0.2.0 tag
emattei Sep 21, 2026
3da0fdd
Finish the sweep: the R help strings and the pysceptre pin
emattei Sep 21, 2026
fb3262a
Add the fixture the R and Python suites share
emattei Sep 24, 2026
c815caf
Simulate power at a fixed element effect, with control cells at exact…
emattei Sep 24, 2026
1549c1c
Stage analysis_mode.tsv for the simulation, and give each sample its …
emattei Sep 24, 2026
48b41d4
Write down what simulated power means, as on the R branch
emattei Sep 24, 2026
8d00177
Solve the pin exactly instead of iterating
emattei Sep 24, 2026
593dd3e
Bring the reference R simulation on this branch to the same spec
emattei Sep 24, 2026
9a2242a
Say plainly what a skipped target does, and retire a stale rationale
emattei Sep 24, 2026
ddefd5c
Keep the pin exact on huge targets; bring the R reference fully in line
emattei Sep 24, 2026
2918aee
Ship pysceptre's scripts/ in the cloud image, and fit the recipe unde…
emattei Sep 25, 2026
900f4d5
Add a samplesheet for moi5 cis on the Python pipeline
emattei Sep 25, 2026
f485883
Run the screen's own permutation test, and keep clamped-theta genes a…
emattei Sep 25, 2026
6570024
Give guides a spread that vanishes at no effect: Beta, sd = c * es * …
emattei Sep 25, 2026
f2d684c
Build the paired-gene set once, not once per gene
emattei Sep 25, 2026
efd0f05
List COMPUTE_POWER's inputs with nullglob, not ls
emattei Sep 25, 2026
fb03d30
Record Stage B at full scale: Python and R agree on moi5 cis
emattei Sep 25, 2026
a85188a
Run a task's (target, replicate) units in parallel, on 8 CPUs and 8 GB
emattei Sep 25, 2026
7aa7eaa
Plan the next step: a random-effect estimand, shared permutations, a …
emattei Sep 25, 2026
30c826e
Plan a re-check of the guide spread's value and form, noise removed
emattei Sep 25, 2026
e6c6db6
Add --permutations per-target: one permutation set per target, drawn …
emattei Sep 25, 2026
40e5371
Record what the data audit settled about the guide spread
emattei Sep 25, 2026
de0eaf3
Make a broken worker pool retryable, and give each worker one thread
emattei Sep 25, 2026
157c4b9
Plan the sparse-route fix first: 2.1-2.3x, bit-identical
emattei Sep 25, 2026
3ea5baa
Add --nulls sparse: pysceptre's draw-matrix route for the permutation…
emattei Sep 25, 2026
36fe149
Plan fit reuse after the fast driver, and retiring --nulls after pysc…
emattei Sep 25, 2026
288cb0f
Add --driver fast: one sparse permutation matrix per target, all gene…
emattei Sep 25, 2026
61c62c6
Add fit reuse: each gene's null model fitted once per simulation, sha…
emattei Sep 25, 2026
596364e
Add --estimand fixed|random: power at a fixed element effect, or aver…
emattei Sep 25, 2026
44e3057
Make the fast configuration the default: per-target permutations, spa…
emattei Sep 25, 2026
9edd600
Compute per-pair power inside each simulation task; keep the per-simu…
emattei Sep 25, 2026
298c464
Size POWER_SIMULATION from a measurement of the new defaults: 8 GB st…
emattei Sep 25, 2026
0651653
Fit null models on the input passed in, not one already in the worker…
emattei Sep 25, 2026
e093248
Record item 5 of the plan as measured: trans 1 h 52 min and $67.57
emattei Sep 26, 2026
2b533fd
Rewrap one line of the plan's status
emattei Sep 26, 2026
0315035
Record Stage 0: Python vs R now equals Python vs Python at another seed
emattei Sep 26, 2026
8e24c18
Merge main (the legacy R pipeline) into the Python pipeline, keeping …
emattei Sep 26, 2026
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
22 changes: 22 additions & 0 deletions .githooks/pre-commit
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@
# * any file larger than 512 KB
# * filenames that are not snake_case
# * camelCase / PascalCase identifiers in R code
# * ruff check and ruff format on Python, when ruff is on PATH
#
# Install: pixi run install-hooks (sets core.hooksPath to .githooks)
# Lint all: pixi run lint (same checks across the whole repo, read-only)
Expand Down Expand Up @@ -140,6 +141,27 @@ while IFS= read -r file; do
rm -f "$scratch"
fi

# --- Python lint and format ------------------------------------------------------------------
# ruff, not a bespoke regex. The R check below exists because nothing else was going to
# enforce R style; Python has a linter that already knows the rules, and the same one the
# project's `pixi run format` task uses, so the hook and the task cannot disagree.
#
# Read-only here even under --fix: ruff format rewrites whole files, and a hook that silently
# reformats what you are committing makes `git diff --cached` stop describing your change.
case "$file" in
*.py)
if command -v ruff >/dev/null 2>&1; then
if ! ruff check --quiet "$file"; then
red "REJECT $file fails ruff check"
failures=$((failures + 1))
elif ! ruff format --check --quiet "$file" >/dev/null 2>&1; then
red "REJECT $file is not ruff-formatted; run 'pixi run -e dev format'"
failures=$((failures + 1))
fi
fi
;;
esac

# --- R identifier convention ---------------------------------------------------------------
# Flags camelCase and PascalCase assignments. Deliberately does not demand strict snake_case:
# SCREAMING_SNAKE constants (PERT_LEVELS) and dotted S3 methods (print.sim_input) are correct
Expand Down
48 changes: 48 additions & 0 deletions .github/workflows/test.yml
Original file line number Diff line number Diff line change
@@ -0,0 +1,48 @@
# The tests that need neither R nor a real screen, which is all of them by default: the two that
# do are behind `-m realdata` and are skipped here.
#
# Not a pixi environment. pixi resolves pysceptre from git and would rebuild the whole conda
# environment on every push; uv installs the same dependencies from pyproject.toml in seconds. What
# that gives up is checking that pixi.lock still resolves, which is what the container build does
# and which does not need to happen on every commit.

name: tests

on:
push:
branches: [main, 'feat/**']
pull_request:

jobs:
test:
runs-on: ubuntu-latest
steps:
- uses: actions/checkout@v4

- uses: astral-sh/setup-uv@v5
with:
enable-cache: true

- name: Install
run: |
uv venv --python 3.12
# pysceptre is not on PyPI; the pin matches pixi.toml's and has to move with it.
uv pip install "pysceptre @ git+https://github.com/broadinstitute/pysceptre.git@v0.2.0"
uv pip install -e ".[dev]"

- name: Lint
run: |
uv run ruff check src tests workflow
uv run ruff format --check src tests workflow

- name: Test
run: uv run pytest -q

- name: Nextflow DAG
# Catches a rewiring mistake without running anything. The stub needs the synthetic
# dataset, which is generated rather than committed.
run: |
curl -fsSL https://get.nextflow.io | bash
uv run watteg-make-test-data --out tests/data/synthetic.h5mu
./nextflow run . -profile standard -params-file config/test.yml -stub-run \
--samplesheet assets/samplesheet_synthetic.csv --outdir /tmp/stub
4 changes: 4 additions & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -88,3 +88,7 @@ tests/data/

# Created by tests/config/run.sbatch -- sparse files whose SIZE is the fixture, 3.3 GB of nothing.
tests/config/inputs/
.venv/
__pycache__/
*.egg-info/
tests/data/*.h5mu
84 changes: 51 additions & 33 deletions README.md
Original file line number Diff line number Diff line change
@@ -1,7 +1,8 @@
# WattEG

**W**atts for **E**lement–**G**ene pairs: power analysis for element–gene pairs in single-cell
CRISPR screens, built on [sceptre](https://katsevich-lab.github.io/sceptre/).
CRISPR screens, built on [pysceptre](https://github.com/broadinstitute/pysceptre) — a Python port
of [sceptre](https://katsevich-lab.github.io/sceptre/)'s statistical engine.

**📖 [Documentation](https://engreitzlab.github.io/WattEG/)**

Expand All @@ -22,67 +23,84 @@ written for one specific analysis. This repository generalises it.

```sh
pixi install
pixi run setup # installs sceptre from a pinned commit
pixi run check-api # verifies that pin against the internals the pipeline uses
```

sceptre is not on conda-forge or bioconda, so it is installed from a pinned commit rather than
captured in `pixi.lock`. The pipeline reads several of its unexported S4 slots, which is why the
version is pinned and checked.
That is the whole of it. The environment is one Python package plus Nextflow: no R, no pinned
sceptre commit, no patch applied to it, and no separate install step. pysceptre enters as a git
dependency pinned to a commit, because a sweep has to be re-runnable against the engine that
produced it.

**The R implementation has not been deleted.** It is the reference the Python path was validated
against and what the methods paper describes, and it still lives in `src/*.R` and `lib/*.R`. What
is gone is the environment that ran it; to run it, use the **`r-implementation`** branch.

That branch is **maintained, not frozen**. Two bugs found during the port change results, and both
were fixed in R as well as in Python rather than being quarantined on a snapshot -- nothing was
published, so there were no numbers owed a bit-for-bit reproduction, and freezing the branch would
only have preserved the bugs. Sweeps produced before 2026-09-21 predate both fixes.
(`legacy` is something else again -- the Snakemake implementation that preceded both.)

## Input

**One file: a sceptre object** (`.rds`) on which `assign_grnas()` and `run_qc()` have been called,
using `grna_integration_strategy = "union"`.
**One file: a `.h5mu` dataset** exported from a sceptre object on which `assign_grnas()` and
`run_qc()` have been called, using `grna_integration_strategy = "union"`.

Everything else is derived from it — the discovery pairs, the gRNA-to-target mapping and the
significance threshold all already live inside the object.
significance threshold all travel with the export.

Making one is a one-off step per dataset, run wherever R and sceptre are available. It is the only
place R appears at all, and it is not part of the pipeline:

```sh
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/
```

If the object's response matrix is odm-backed (out-of-core), pass `--response-odm` on
`prepare_sim_input.R` with the path to the backing `.odm` file — see [Usage](https://engreitzlab.github.io/WattEG/usage/).
**`--all-cells` is required, not optional.** DESeq2 "poscounts" size factors are a per-cell
reduction against a per-gene geometric mean taken over every cell in the object, so an export
restricted to the QC-passing cells gives different size factors for the cells that remain. The
export handles odm-backed (out-of-core) matrices itself, which is why there is no longer an
`--response-odm` flag anywhere in the pipeline.

List your samples in a CSV (see `assets/samplesheet.csv`):

```csv
sample,sceptre_object
sample1,results/sample1/sceptre_object.rds
sample,dataset
sample1,export/dataset.h5mu
```

Add an optional `response_odm` column (see `assets/samplesheet_seqera_test.csv`) for any sample whose
sceptre object's response matrix is odm-backed — same file `prepare_sim_input.R --response-odm` takes
standalone. Leave it blank, or omit the column entirely, for in-memory-backed objects.

## Quickstart

Each step is a standalone executable in `src/` with `--help`.
```sh
nextflow run . -profile sherlock -params-file config/config.yml
```

Each step is also a console script with `--help`, so a sweep can be driven by hand or by another
runner:

```sh
# derive the simulation inputs (once per sample)
Rscript src/prepare_sim_input.R \
--sceptre-object results/sample1/sceptre_object.rds \
--outdir prepared/
watteg-prepare-sim-input --dataset export/dataset.h5mu --outdir prepared/

# split targets into per-task chunks
Rscript src/split_pairs.R --pairs prepared/pairs.tsv --n-splits 280 --outdir splits/
watteg-split-pairs --pairs prepared/pairs.tsv --n-splits 280 --outdir splits/

# simulate (once per split x effect size)
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 \
watteg-run-power-simulation \
--prepared prepared/ --pairs splits/split_001.tsv \
--effect-size 0.15 --reps 100 --seed 20250812 \
--out sim/split_001_es0.15.tsv

# power per pair, then one table across effect sizes
Rscript src/compute_power.R \
--simulations "$(ls sim/*_es0.15.tsv | paste -sd, -)" \
--threshold-file prepared/discovery_threshold.txt \
watteg-compute-power \
--simulations sim/ --threshold-file prepared/discovery_threshold.txt \
--out power_es0.15.tsv

Rscript src/summarize_power.R \
--power power_es0.15.tsv,power_es0.2.tsv \
--out power_summary.tsv
watteg-summarize-power \
--power power_es0.15.tsv power_es0.2.tsv \
--sim-input prepared/sim_input.h5 --out power_summary.tsv
```

Parameters live in `config/config.yml`.
Expand Down
4 changes: 2 additions & 2 deletions assets/samplesheet.csv
Original file line number Diff line number Diff line change
@@ -1,2 +1,2 @@
sample,sceptre_object
day0_grna20_no_shuffle,results/day0_grna20_no_shuffle/sceptre_object.rds
sample,dataset
day0_grna20_no_shuffle,results/day0_grna20_no_shuffle/dataset.h5mu
10 changes: 5 additions & 5 deletions assets/samplesheet_all_samples.csv
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
sample,sceptre_object
day0_grna20,results/day0_grna20/sceptre_object.rds
day0_grna20_no_shuffle,results/day0_grna20_no_shuffle/sceptre_object.rds
day2_grna20,results/day2_grna20/sceptre_object.rds
day4_grna20,results/day4_grna20/sceptre_object.rds
sample,dataset
day0_grna20,results/day0_grna20/dataset.h5mu
day0_grna20_no_shuffle,results/day0_grna20_no_shuffle/dataset.h5mu
day2_grna20,results/day2_grna20/dataset.h5mu
day4_grna20,results/day4_grna20/dataset.h5mu
4 changes: 2 additions & 2 deletions assets/samplesheet_day0_grna20.csv
Original file line number Diff line number Diff line change
@@ -1,2 +1,2 @@
sample,sceptre_object
day0_grna20,results/day0_grna20/sceptre_object.rds
sample,dataset
day0_grna20,results/day0_grna20/dataset.h5mu
4 changes: 2 additions & 2 deletions assets/samplesheet_moi5_cis.csv
Original file line number Diff line number Diff line change
@@ -1,2 +1,2 @@
sample,sceptre_object,response_odm
moi5,gs://hgrm-seqera/element-gene-power-analysis-union/moi5/sceptre_object.rds,gs://hgrm-seqera/element-gene-power-analysis-union/moi5/response.odm
sample,dataset
moi5,gs://hgrm-seqera/element-gene-power-analysis-union/moi5/dataset.h5mu
4 changes: 2 additions & 2 deletions assets/samplesheet_moi5_trans.csv
Original file line number Diff line number Diff line change
@@ -1,2 +1,2 @@
sample,sceptre_object,response_odm
moi5_trans_discovery,gs://hgrm-seqera/sceptre-nf-ondisc-test/moi5-trans-output/sceptre_object.rds,gs://hgrm-seqera/20260508-wtc11-sceptre/grna-single-target-moi5-ondisc/gene.odm
sample,dataset
moi5_trans_discovery,gs://hgrm-seqera/sceptre-nf-ondisc-test/moi5-trans-output/dataset.h5mu
4 changes: 2 additions & 2 deletions assets/samplesheet_seqera_test.csv
Original file line number Diff line number Diff line change
@@ -1,2 +1,2 @@
sample,sceptre_object,response_odm
wtc11_union_odm,gs://hgrm-seqera/element-gene-power-analysis-test/sceptre_object.rds,gs://hgrm-seqera/element-gene-power-analysis-test/response.odm
sample,dataset
wtc11_union_odm,gs://hgrm-seqera/element-gene-power-analysis-test/dataset.h5mu
4 changes: 2 additions & 2 deletions assets/samplesheet_synthetic.csv
Original file line number Diff line number Diff line change
@@ -1,2 +1,2 @@
sample,sceptre_object
synthetic,tests/data/sceptre_object_high_default.rds
sample,dataset
synthetic,tests/data/synthetic.h5mu
Loading
Loading