|
| 1 | +# Off-target CRISPR edit finder — a 5-minute demo |
| 2 | + |
| 3 | +A hands-on tour you can run on a login node in about a minute. Nothing here needs cohort |
| 4 | +data, a cluster allocation, or the DRAGEN reference: every command below runs against a |
| 5 | +**5.5 KB synthetic genome** committed to the repo. |
| 6 | + |
| 7 | +For the full scientific writeup — what the workflow is and is not, how the shape model was |
| 8 | +trained, what the AAVS1 cohort actually showed — see |
| 9 | +[`OFFTARGET_WORKFLOW.md`](OFFTARGET_WORKFLOW.md). This page is the "show me it working" one. |
| 10 | + |
| 11 | +--- |
| 12 | + |
| 13 | +## What the finder does, in three sentences |
| 14 | + |
| 15 | +CRISPR cuts where you aimed it, and sometimes where you didn't. This workflow takes deep |
| 16 | +error-corrected sequencing (ECS) at a panel of candidate sites and finds the reads that |
| 17 | +carry an edit, then asks whether ordinary 30x whole-genome sequencing of the *same* sample |
| 18 | +could have found those same edits on its own. |
| 19 | + |
| 20 | +Two arms, one join: **ECS is the truth arm**, **WGS is the arm under test**. |
| 21 | + |
| 22 | +--- |
| 23 | + |
| 24 | +## Setup (once) |
| 25 | + |
| 26 | +```bash |
| 27 | +# nextflow + apptainer on the RIS cluster |
| 28 | +S=/storage2/fs1/dspencer/Active/spencerlab/apps/modules/spack/software/linux-rhel9-x86_64_v4/gcc-13.3.0 |
| 29 | +export PATH="$S/nextflow-25.10.4-fz3mazhspz5uruz5ymyxhcwcauohjgfm/bin:$S/apptainer-1.4.5-ifu7outvqr62l343tqfxuxcguzftiuqe/bin:$PATH" |
| 30 | + |
| 31 | +# nf-test (once, anywhere on PATH) |
| 32 | +curl -fsSL https://code.askimed.com/install/nf-test | bash |
| 33 | +``` |
| 34 | + |
| 35 | +--- |
| 36 | + |
| 37 | +## Demo 1 — the whole test suite, real execution, ~1 minute |
| 38 | + |
| 39 | +```bash |
| 40 | +nf-test test --profile stub,apptainer |
| 41 | +``` |
| 42 | + |
| 43 | +``` |
| 44 | +OFFTARGET_METRICS |
| 45 | + Test 'computes the real metrics and emits every declared channel' PASSED |
| 46 | + Test '--betas reaches the script: a custom beta set changes the columns' PASSED |
| 47 | + Test '--negatives all_label0 is honoured as a distinct negative set' PASSED |
| 48 | +ECS_INDELS |
| 49 | + Test 'calls the edit and reports the right VAF — tagged BAM off by default' PASSED |
| 50 | + Test 'offtarget_tagged_bam=true emits a per-read tagged BAM and its index' PASSED |
| 51 | +
|
| 52 | +SUCCESS: Executed 5 tests |
| 53 | +``` |
| 54 | + |
| 55 | +These are not stubs. Each one starts the real container and runs the real Python. |
| 56 | + |
| 57 | +> **Why that distinction matters.** The repo also has `-stub-run` tests, and they are |
| 58 | +> useful — but a stub process body only `touch`es its output files. It never executes a |
| 59 | +> line of `bin/*.py`. A module that passes a flag the script does not accept stays green |
| 60 | +> through every stub test in the repo. That exact bug class once survived CI and then died |
| 61 | +> three hours into a cluster run. nf-test is the tier that catches it. |
| 62 | +
|
| 63 | +--- |
| 64 | + |
| 65 | +## The toy genome |
| 66 | + |
| 67 | +Everything below runs on `tests/fixtures/ecs/` — 5.5 KB, regenerable with |
| 68 | +`python tests/fixtures/make_ecs_fixture.py`: |
| 69 | + |
| 70 | +``` |
| 71 | +chr1, 3000 bp |
| 72 | + target A @ 1001 the EDITED site — 12 of 20 evaluable read pairs carry a 5 bp deletion |
| 73 | + target B @ 1201 a QUIET site — no edit |
| 74 | +
|
| 75 | +edited.cram 35 pairs / 70 records |
| 76 | +control.cram 20 clean pairs (the matched normal) |
| 77 | +``` |
| 78 | + |
| 79 | +The 35 pairs are awkward on purpose: |
| 80 | + |
| 81 | +| pairs | what they are | what should happen | |
| 82 | +|---|---|---| |
| 83 | +| 12 | carry a 5 bp deletion | counted as edited | |
| 84 | +| 8 | clean, span the target | counted as reference | |
| 85 | +| 2 | duplicate-flagged | **excluded** from the denominator | |
| 86 | +| 2 | MAPQ 3 | **excluded** (floor is 20) | |
| 87 | +| 2 | NM 6 | **excluded** (ceiling is 4) | |
| 88 | +| 9 | sit where the two ±150 bp windows **overlap** | visited by both targets, must be counted once | |
| 89 | + |
| 90 | +That last row is the interesting one: targets A and B are 200 bp apart, so their fetch |
| 91 | +windows overlap. Any per-read output has to survive being visited twice. |
| 92 | + |
| 93 | +--- |
| 94 | + |
| 95 | +## Demo 2 — find the edit |
| 96 | + |
| 97 | +```bash |
| 98 | +IMG=/storage2/fs1/dspencer/Active/spencerlab/abonney/apptainer_cache/ghcr.io-dhslab-docker-scge-offtarget-260710.img |
| 99 | +F=tests/fixtures/ecs |
| 100 | + |
| 101 | +apptainer exec -B /storage2 "$IMG" python bin/find_edited_reads.py \ |
| 102 | + --fasta $F/ref.fa \ |
| 103 | + --edited-bam $F/edited.cram --control-bam $F/control.cram \ |
| 104 | + --target-file $F/targets.vcf \ |
| 105 | + -o demo.tsv |
| 106 | +``` |
| 107 | + |
| 108 | +``` |
| 109 | +chrom start end total_reads indel_reads indel_fraction control_reads control_indel_reads |
| 110 | +chr1 1000 1001 20 12 0.6 20 0 |
| 111 | +chr1 1200 1201 23 0 0.0 20 0 |
| 112 | +``` |
| 113 | + |
| 114 | +Read that as: **at the edited site, 12 of 20 evaluable reads carry the deletion (VAF 0.60), |
| 115 | +and the matched control has none of them.** That last column is what makes the call somatic |
| 116 | +rather than germline. The quiet site stays quiet. |
| 117 | + |
| 118 | +Note `total_reads` is 20, not 26 — the duplicate, low-MAPQ and high-mismatch pairs were |
| 119 | +dropped before the denominator was formed. If a filter ever regresses, that number climbs |
| 120 | +and the VAF silently falls. The test asserts on it for exactly that reason. |
| 121 | + |
| 122 | +--- |
| 123 | + |
| 124 | +## Demo 3 — per-read tags for IGV *(new)* |
| 125 | + |
| 126 | +Add one flag and every read gets an `XC` tag naming how the caller classified it: |
| 127 | + |
| 128 | +```bash |
| 129 | +apptainer exec -B /storage2 "$IMG" python bin/find_edited_reads.py \ |
| 130 | + --fasta $F/ref.fa \ |
| 131 | + --edited-bam $F/edited.cram --control-bam $F/control.cram \ |
| 132 | + --target-file $F/targets.vcf \ |
| 133 | + --tagged-bam-out demo.tagged.bam -o demo.tsv |
| 134 | +``` |
| 135 | + |
| 136 | +``` |
| 137 | + 31 Unedited_WT |
| 138 | + 15 Skipped_NoSpan |
| 139 | + 12 Edited_Deletion_5bp |
| 140 | + 4 Skipped_Duplicate |
| 141 | + 4 Skipped_LowMapQ |
| 142 | + 4 Skipped_Mismatches |
| 143 | +────────────────────────── |
| 144 | + 70 records, 35 unique read names |
| 145 | +``` |
| 146 | + |
| 147 | +Load `demo.tagged.bam` in IGV and use **Color alignments by → tag → XC**. A reviewer can now |
| 148 | +see *why* a read was or wasn't counted, instead of taking the caller's word for it. |
| 149 | + |
| 150 | +Two things worth pointing at on a slide: |
| 151 | + |
| 152 | +- **70 records, 35 unique names.** Every read is written exactly once, even the 9 pairs in |
| 153 | + the window overlap that two different targets both visited. |
| 154 | +- **The `Skipped_` classes are visible.** Reads that were excluded are in the file and |
| 155 | + labelled, not silently missing. |
| 156 | + |
| 157 | +In the pipeline this is `--offtarget_tagged_bam`, and it is **off by default**: on a real |
| 158 | +AAVS1 sample (1149 targets, ~11000x) it produces an 800 MB BAM and pushes the caller's peak |
| 159 | +memory from 1.4 GB to 6.3 GB. Turn it on for review, not for a cohort sweep. |
| 160 | + |
| 161 | +--- |
| 162 | + |
| 163 | +## Demo 4 — clinical metrics *(new)* |
| 164 | + |
| 165 | +```bash |
| 166 | +apptainer exec -B /storage2 "$IMG" python bin/offtarget_metrics.py \ |
| 167 | + --training tests/fixtures/training_mini.tsv \ |
| 168 | + --betas 2,5 --out-json m.json --out-txt m.txt && cat m.txt |
| 169 | +``` |
| 170 | + |
| 171 | +``` |
| 172 | +-- RANKING (score as a continuous ranker; ECS-label denominator) -- |
| 173 | + n_pos / n_neg : 6 / 4 prevalence 0.6000 |
| 174 | + PR-AUC : 0.8552 (1.43x prevalence) |
| 175 | + excluded, no score : 1 (1 pos / 0 neg) — INSUFFICIENT COVERAGE, never filled with 0 |
| 176 | +
|
| 177 | +-- OPERATING POINT (verdict contains 'LIKELY EDIT'; ECS-label denominator) -- |
| 178 | + TP 5 FP 2 FN 1 TN 2 |
| 179 | + precision : 0.7143 recall : 0.8333 recall incl. unevaluable : 0.7143 |
| 180 | + F1 0.7692 F2 0.8065 F5 0.8280 |
| 181 | +``` |
| 182 | + |
| 183 | +Three things this output is built to stop you getting wrong: |
| 184 | + |
| 185 | +**1. F-beta is recall-weighted, and you can see it.** Here recall (0.833) beats precision |
| 186 | +(0.714), so F rises with beta: F1 < F2 < F5. For a screening assay a missed off-target edit |
| 187 | +costs far more than a followed-up false one, which is the whole argument for reporting F2/F5 |
| 188 | +rather than F1. F-beta is monotone in beta, so F1 is always an *endpoint* — if you ever see |
| 189 | +F1 in the middle, the weights got applied backwards. |
| 190 | + |
| 191 | +**2. The uncovered site is excluded, not counted as zero.** One credible positive had no WGS |
| 192 | +coverage at all. It is dropped from the ranking and reported separately, and `recall incl. |
| 193 | +unevaluable` (0.7143) is printed next to `recall` (0.8333) so the depth floor can never hide |
| 194 | +inside a good-looking number. |
| 195 | + |
| 196 | +**3. Precision is against the ECS label, never against manual review.** The report says so |
| 197 | +out loud, in the file. Manual review has 55 confirmed positives and 80,440 NaNs, where NaN |
| 198 | +means *unreviewed*, not *rejected* — it contains no confirmed negatives. Any precision |
| 199 | +computed from it would score every discovery the reviewers never reached as a false |
| 200 | +positive. Recall against manual review is reported separately, and as recall only, by |
| 201 | +`bin/validate_recall.py`. |
| 202 | + |
| 203 | +--- |
| 204 | + |
| 205 | +## Where the tests live and what each tier is for |
| 206 | + |
| 207 | +```bash |
| 208 | +bash run_offtarget_tests.sh # all five tiers, ~3 min |
| 209 | +``` |
| 210 | + |
| 211 | +| tier | what it runs | what it proves | catches | |
| 212 | +|---|---|---|---| |
| 213 | +| 0 | `ast.parse`, `nextflow -preview` | it is syntactically valid | typos | |
| 214 | +| 1 | `pytest tests/*.py` | **the science is right** | wrong math, wrong joins | |
| 215 | +| 2 | unpickle the shape model | the sklearn pin holds | version drift | |
| 216 | +| 3 | `nextflow -stub-run` | **the DAG wires up** | broken channels | |
| 217 | +| 4 | `nf-test` | **the two halves connect** | module ↔ script contract | |
| 218 | +| 5 | real AAVS1 cohort (SLURM) | it works on real data | everything else | |
| 219 | + |
| 220 | +Tiers 1 and 3 are the two halves that tier 4 joins: pytest runs the Python but never through |
| 221 | +Nextflow; the stub run goes through Nextflow but never runs the Python. |
| 222 | + |
| 223 | +Run one feature's tests only: |
| 224 | + |
| 225 | +```bash |
| 226 | +nf-test test --tag metrics --profile stub,apptainer |
| 227 | +nf-test test --tag tagged_bam --profile stub,apptainer |
| 228 | +``` |
| 229 | + |
| 230 | +--- |
| 231 | + |
| 232 | +## Honest limits |
| 233 | + |
| 234 | +- The fixtures are **synthetic**. They prove the plumbing and the arithmetic, not that the |
| 235 | + caller is right about real CRISPR biology — that is what the AAVS1 cohort run is for. |
| 236 | +- The metrics fixture is hand-built to produce a specific confusion matrix. Its PR-AUC is |
| 237 | + not a performance claim about the model; it is a fixed number chosen so the test can |
| 238 | + assert on it. |
| 239 | +- Two blind spots remain unmeasured on real data: the **sub-5% VAF floor**, and a **de-novo |
| 240 | + off-target positive control**. Neither is addressed by anything on this page. |
0 commit comments