Feat/offtarget wgs - #18
Open
Bonney96 wants to merge 79 commits into
Open
Conversation
…eration to include new parameters for plots and improved HTML rendering in the report.
…a chromosome. Issue: The script calculated query windows by subtracting a window size from the genomic position (start = pos - window). For sites near the beginning of a chromosome (e.g., pos < 150), this resulted in a negative start coordinate (e.g., -65), causing pysam.fetch() to raise a ValueError: start out of range. Fix: Added max(0, ...) clamping to all fetch() calls to ensure start coordinates are never negative.
fixed bug that aggregated results for all sites in the same loop
…d_normal_counts - Move indel record list inside the per-interval loop so it resets each iteration; previously vcf_dicts accumulated across all intervals, causing duplicate/incorrect VCF output - Clamp fetch() start to max(0, ...) in add_normal_counts to prevent ValueError for variants near chromosome start - Rename vcf_results/vcf_dicts to all_indel_records/indel_vcf_records for clarity
Self-contained WGS off-target arm; leaves the default SCGE workflow untouched (selected via named entry OFFTARGET). Uses the pileup shape-model approach (legacy crispr_ml deprecated for off-target work). Flow: samplesheet(datatype ecs|wgs, guide) -> ECS_INDELS (find_edited_reads.py = error-corrected truth VAF) WGS_WORKLIST -> PON_OFFTARGET_FILTER (genome-wide homology-free worklist) HOTSPOT_TO_TABLE -> SCORE_HOTSPOTS -> BUILD_TRAINING_TABLE (WGS feats x ECS VAF) RECALL_VS_VAF (WGS recall vs ECS VAF + trustworthy-VAF floor) RECONCILE_OFFTARGET_REPORT (is_hotspot / ecs_confirmed annotation) New bin/: hotspot_to_table, join_training_table, recall_vs_vaf, reconcile_offtarget_report (+ core worklist_from_vcf/pon_filter/score/ features/pileup_snapshot). Model asset assets/models/wgs_shape_model.pkl. offtarget_* params in nextflow.config; publishDir via withName in conf/modules.config -> outdir/offtarget.
tests/: 5 pytest cases covering the OFFTARGET glue scripts end-to-end on synthetic fixtures (a real edit, a true negative, a below-floor low-VAF edit) - the ECS x WGS training-table join, recall-vs-VAF floor behaviour, worklist hotspot reconciliation, and the coordinate-mismatch guardrail. Run with `pytest tests/` under any interpreter with pandas. README: replace the untouched nf-core template (unfinished intro, TODO placeholders, stale samplesheet columns) with accurate pipeline summary, both entry points (SCGE + -entry OFFTARGET), correct samplesheet formats, params, containers, and testing/reference-data sections.
docs/offtarget_metro.mmd: Mermaid + %%metro source with three lines (ECS truth, WGS genome-wide discovery, WGS-only model). Rendered with nf-metro to light/dark SVGs and an interactive HTML (pan/zoom, per-line filter). Embedded in README and docs/OFFTARGET_WORKFLOW.md via a light/dark <picture>. Also fixed a stale branch name in the doc.
…ADME - Add whole-pipeline SCGE metro map (nf-metro) with animated line traffic, light/dark SVG + interactive HTML - Re-render off-target metro map with animation - README: swap static DAG for animated SCGE metro, drop stale reference-data section, replace DRAGEN container row with CRISPR_ML (TBD), link Spencer Lab
- README: restore the static scge_analysis_DAG.jpg for the default workflow (metro metaphor added little for a near-linear DAG) - Remove the now-unused SCGE metro sources/renders - Off-target metro: drop the redundant samplesheet->model edge that caused the green self-crossing, and widen station spacing so the model arm labels no longer collide
- features.check_sklearn_version(): raise on a backward load (runtime older than the model's pickle version, the silent mis-score case), warn on forward loads. Wired into score.py, worklist_from_vcf.py, find_edited_reads.py after joblib.load. - Point the 8 OFFTARGET modules at ghcr.io/dhslab/docker-scge-offtarget:260710 (dated tag for reproducible builds).
- README: replace the 'genome-wide homology-free discovery / recovers edits from WGS alone' claim with the measured scope — hotspot edit-confirmation + PoN worklist, WGS-only calls validated >=5% VAF (ROC-AUC 0.82, recall 0.92). - OFFTARGET_WORKFLOW.md: add 'what it is / is not', a measured-accuracy table (from read_cnn/pileup outputs), and the open gaps (sub-5% floor unmeasured; off-target discovery has no positive control; paired glue arm validated only end-to-end on the first RIS run).
…reference workflows/scge.nf validates params.refdir at parse time via a top-level statement, so it runs on EVERY entry — OFFTARGET and HOTSPOTS included, even though neither touches DRAGEN. The default points at /storage2, which exists on the cluster and nowhere else, so all three stub jobs died on GitHub Actions with 'DRAGEN hash table config not found' before a single process was scheduled. Predates this branch (added 2026-01-26); it is invisible locally precisely because the path resolves there. Nulling refdir in the stub profile short-circuits the check (Groovy && on a null left operand), leaving the production validation untouched. Verified against all five stub invocations the CI workflow runs, with refdir pointed at a nonexistent path.
Feat/offtarget prauc metrics attempt 2
…ork off-cluster Two fixes that both kept the metrics work from actually running. 1. workflows/offtarget.nf was left out of the original commit's `git add`, so OFFTARGET_METRICS was defined, parameterised and schema'd but never invoked — the process existed and nothing called it. The stub run appeared to exercise it only because the wiring was present uncommitted in the working tree. 2. Complete the stub reference-path fix. workflows/scge.nf validates every reference path at PARSE time through top-level `params.X ? Channel.fromPath(X, checkIfExists: true) : []` statements that run on EVERY entry, including OFFTARGET and HOTSPOTS, neither of which touches DRAGEN or VEP. All defaults point at /storage2, so a stub run dies before scheduling a process anywhere that path is absent. Nulling refdir alone just exposed the next path in the list; this covers the whole set. Two are deliberately not null: crisprme_index_dir is file()'d unconditionally when run_crisprme is on, and the OFFTARGET stub steps do not pass their own --fasta. These stub jobs have never passed on CI: the parse-time checks landed 2025-09-30 and the offtarget CI workflow 2026-07-21, and the failure reproduces at the commit that introduced the workflow. Verified across all six stub invocations the workflow runs, with /storage2 simulated absent, and re-checked with /storage2 present so cluster behaviour is unchanged.
…eview The ECS caller decides per read whether it supports an edit, aggregates to per-site counts and throws the per-read decision away, so a reviewer sees indel_fraction = 0.42 and has to re-derive by eye which reads made it up. Add --tagged-bam-out to bin/find_edited_reads.py: every read over a target window gets a string tag (default XC) naming its classification, so IGV can colour the pileup by the caller's own reasoning. Two constraints shaped the implementation: * The per-target read loop fetches with a +/-150 bp pad, so overlapping targets make it visit the same alignment record repeatedly. Writing in-loop would emit duplicate records, which IGV renders as inflated depth. Tags are therefore only accumulated in the loop, keyed on (query_name, flag, reference_start) -- query_name alone collides between mates and between primary/supplementary records -- and the BAM is written in a second pass, sorted and indexed. Conflicting calls for one record resolve by a documented precedence (edit > reference > unevaluable > padding). * Output is window-restricted and off by default. On one AAVS1 site14 ECS sample (1149 targets, ~11,000x) it is 800 MB / 21.8 M reads and lifts peak RSS from 1.4 GB to 6.3 GB, so ECS_INDELS gets 24 GB instead of 8 GB while the flag is set. The dominant classification is left implicit rather than stored, which keeps the tag map roughly 40% smaller. Verified on real data: with tagging off the TSV is byte-identical to the unmodified script; with it on, the full-panel BAM has 21,792,054 records, zero missing tags, zero duplicate keys, and is coordinate-sorted and indexable. Note the tag counts do not equal the TSV columns and are not meant to: tags are per alignment record, indel_reads is per fragment and post-filter. Documented. Wires params.offtarget_tagged_bam through ECS_INDELS and GET_INDELS, adds the schema entry, docs with the IGV recipe, and six tests on a new synthetic aligned-read fixture covering dedup, sort order, skip reasons and precedence.
offtarget: wire OFFTARGET_METRICS into the workflow; make stub runs w… attempt 3
offtarget: optional per-read XC tags on a target-window BAM for IGV r…
The plot lagged the CSV. Two things were wrong on the figure only:
- the y-axis read "WGS recall (score >= hi)", but detection is the reported
verdict, which includes high-evidence rescues scoring below that threshold.
- the per-bin annotation was n=<credible total>, which reads as "n/n detected"
and hides the depth floor: a bin of 4 credible edits with 2 uncovered plots
at recall 1.0 off 2 evaluable sites, not 4.
Label the evaluable count instead (with "+N uncov." when it differs), and stamp
the credible-set definition into the footer -- the figure travels further than
the CSV does, so it has to carry its own denominator.
No behaviour change; n_evaluable/n_unevaluable/n_excluded were already computed.
Agent tooling keeps per-repo settings and full worktree checkouts there -- each with its own .nextflow logs and work dirs, 663 MB in this checkout. Untracked and unignored, it was one `git add -A` away from being committed.
…tion, add audit columns Three gaps in the caller, all surfaced while reviewing the automated off-target queue. control_bnd_reads: BND calls had no matched-control counterpart, so a junction present in the unedited control looked identical to one that was not. Indels have carried control counts all along; breakends did not. BND position off-by-one: add_normal_counts was measuring breakend support one base off. The fix is deliberately BND-only -- a self-consistency run (--edited-bam == --control-bam) shows current handling is correct for indels, and applying the same +1 globally makes indels six times worse. On CTLA4_1: bnd_count 233 -> 210, control_bnd_reads 0 -> 197, with indel_count unchanged at 42. n_control_filtered / min_cut_distance: -x/--max-in-control silently removes events rather than reporting them, so a table gave no way to tell an absent site from a filtered one. These two columns make that visible without re-running. Header is now 19 columns.
…kends Lower the review gate from indel_reads>=10 & indel_fraction>=0.05 to >=2 & >=0.005. The old gate was the lab's historical Excel filter; 0.5% is the threshold we actually want to review at. On the 32-sample CAR-T WGS cohort this takes gated rows from 160 to 1,498 and the queue from 81 to 88, losing nothing. Note the reads>=2 half is inert on its own: 837 rows sit at exactly 2 indel reads and rule 3 needs 3 distinct lengths, so none can pass. See docs/THRESHOLD_CHANGE_IMPACT.md. pon_min_reads 3 -> 2, which must stay at or below review_min_reads. When the panel needs more evidence than the gate consuming it, there is a band where the filter admits a site the panel structurally cannot flag. That is how chr18:61,959,927 -- recurrent under 23 guides in 28 of 32 samples, 2 control reads -- reached the queue three times. Rule 6, systematic noise from the DRAGEN panel. A bare interval hit is not usable: it flags 33.1% of artifacts but also 1.2% of real on-target edits. Requiring both several panel donors and a D/I allele code makes it safe -- the D/I condition alone protects all 81 on-target sites at every donor threshold tested. Off unless --review_snv_noise is given, because it changes a finding: the off-target list goes 7 -> 6, dropping chr1:28,580,333 where 17 unrelated donors show deletion noise reaching 7.8% against our 9.4% call. REVIEW_FILTER_BND applies the same triage to breakends, replacing length diversity with breakpoint promiscuity since a junction has two ends. On the cohort: 2,404 junctions -> 36 gated -> 25 queued, being 23 multi-cut deletions across 7 samples spanning 413 bp to 125 kb and 2 deletions at a cut site. The promiscuity rule is currently non-binding -- the BND position fix dropped the worst hub from 25 partners to 3 -- and is kept as insurance for future cohorts. Only the WGS v3.1.0 SV panel is safe to wire into --review_sv_noise. IDPF_WGS v3.0.0 flags 25 of 25 real junctions and FF_Heme v3.1.0 flags 24 of 25; either would erase the entire breakend result. Also corrects schema defaults that had drifted from nextflow.config (review_min_reads, review_min_vaf, pon_min_reads), which nf-core lint compares.
A failed grid *submission* kills the whole Nextflow run, not just one task. The
2026-08-11 cohort run lost 160 pending tasks to a single
sbatch: error: Batch job submission failed: Connection reset by peer
after 104 tasks had already succeeded, when the RIS SLURM controller dropped
connections mid-run.
Nextflow's default executor.retry.reason matches ONLY "Socket timed out on
send/recv operation", so that error was never retried. Widen it to the
connection failures the controller actually emits, with a bounded backoff.
A 32-sample run makes hundreds of sbatch calls over many hours, so a controller
blip is a question of when, not if.
Rule 4 asked "is this site in the panel of normals?" -- an existence test. That
is the wrong shape for the question: a site can carry 1% background in the
controls and 40% signal in the treated sample, and a blacklist throws that edit
away. It also needs a cohort, which the single-sample user does not have.
Replace it with a test of expectation:
AQ = -10 log10 P(X >= k | n, background at this locus)
with the background an empirical-Bayes Beta posterior over the control's own
reads. Measured on the 32-sample CAR-T WGS cohort against the curated two-class
review label, --noise-model matched equals the panel of normals exactly:
matched control 91 rows 64 confirmed 9 rejected precision 0.877
panel of normals 88 rows 64 confirmed 9 rejected precision 0.877
while needing only the sample's OWN unedited control -- no cohort at all. AUC
0.977 / AP 0.944. Germline falls out for free (0 of 132 germline-like sites
reach the cut) without a dedicated rule.
The external DRAGEN panel cannot do this job and the measurement says why: only
3.6% of gated loci have an indel-capable record in it, so the floor rather than
the panel decides the outcome, and precision lands at 0.31.
Two hazards handled explicitly:
* The clean-control trap. A control with no alt reads at depth d shows the
background is below ~1/d, NOT that it is zero. Left alone the posterior mean
collapses to 3.5e-5 against a 170x control, ~100x past what 170 reads support,
inflating AQ at every locus the control never sampled deeply. apply_depth_floor
bounds the posterior mean at 1/control_depth. Harmless at the current
VAF>=0.005 gate (2 gated rows sit below a 170x control's 0.59% resolution) and
not harmless the moment that gate is lowered for a high-sensitivity run.
* Rule 4 silently having no source. No PoN, no noise model and one guide in the
invocation means cross-guide recurrence cannot fire, and the unfiltered queue
looks exactly like a precise one. That now warns loudly and can be made fatal
with --strict-fallback.
Default is off, so existing runs are unchanged: the 88 and 87 baselines
reproduce and the 87-row queue is byte-identical to the pre-change output.
Numbers, arms and reproduction commands in docs/NOISE_MODEL_EXPERIMENT.md.
…GB panel
Follow-on to the noise model. With the panel of normals gone, nothing left in
review_filter.py needs a cohort view, so the two-stage split had to justify
itself again. It does, on cost asymmetry:
REVIEW_FILTER, whole cohort, from TSVs 1.78 s
GET_INDELS, per sample, from CRAM 50 min - 1h 28m
Every threshold in this filter was calibrated by re-running it. Moving those
decisions into the caller turns each experiment from seconds into hours and
deletes the why_dropped audit trail -- including the 753 of 1,498 rows sitting
11-25 bp from the cut, which are the evidence for the 10 bp threshold. The
intermediate TSVs are 27 MB in total, so there is no I/O case either.
So computation moves to the caller and decisions stay here. add_features() now
consumes caller-supplied columns and derives them only as a fallback, which
keeps tables from older callers working -- verified byte-identical output on the
pre-annotation tables.
Also exposes the cut-distance trap rather than walking into it.
min_cut_distance (caller) and cut_dist_min (filter) are NOT the same quantity:
the first is min(|pos-PAM|, |pos+len(ref)-1-PAM|), the second is derived from
indel_info and measured from the anchor base only, so always larger. They
disagree on 162 of 1,498 gated rows -- rule 2 drops 838 vs 682 -- though both
retain all 64 confirmed edits. The caller's is the more correct metric, but the
10 bp threshold was calibrated against the other, so it is opt-in via
--cut-dist-source and defaults to the calibrated one.
bin/subset_noise_panel.py cuts the panel to the positions a guide can query:
97,851,753 records -> 4,891 (60 KB), byte-identical filter output. A target file
is per-guide, not per-sample, so this is what would make a per-sample panel
lookup affordable. It also retires the tabix question outright.
Three loose ends, closed together because they are the same decision: the pipeline had two cut-distance rules and two rule-4 implementations, and only one of each is needed. ONE DISTANCE RULE, OWNED BY THE CALLER -------------------------------------- find_edited_reads.py -d drops to 10 (was 25) and is now the pipeline's only cut-distance threshold. It filters on Distance = min(|pos-PAM|, |pos+len(ref)-1-PAM|), the true minimum, which is also what min_cut_distance reports -- and because -d is applied before min_cut_distance is computed, every surviving event is within 10 bp by construction. review_filter.py therefore stops deriving its own distance from indel_info. That was a different quantity (measured from the anchor base only, so always larger) wearing the same name, and the two disagreed on 162 of 1,498 gated rows. Rule 2 now reads the caller's value and is a consistency guard that drops nothing. Recalibrated before changing the default, against the curated review label: thresholds 6, 8, 10, 12 and 15 all hold 64/64 confirmed edits at precision 0.877; only 25 degrades it (0.831). 10 is inside that plateau and matches the cut biology. THE PANEL OF NORMALS IS OFF --------------------------- review_noise_model = 'matched' and review_auto_pon = false. The beta-binomial needs only the sample's own control, so the cohort-scale PoN buys nothing that a single-sample submission can use. PON_SCORE and BUILD_PON leave the DAG, which also removes a per-control scoring step from every run. Quality is unchanged on the same 32 tables -- precision 0.889, 64/64 confirmed, 70/70 ECS-confirmed edits including 12/12 below 5% VAF, 81/81 on-target -- for a queue of 96 rather than 87. The 9 extra rows are unlabelled candidates the curated review never adjudicated, i.e. the cost of the more permissive correct distance metric, not lost precision. Tables written before min_cut_distance existed still run: they fall back to the old derivation, warn that rule 2 will be stricter than the caller's, and reproduce the previous 87-row queue byte-identically. NOTE: changing -d alters the caller's task hash, so the next run re-executes all 32 GET_INDELS tasks rather than resuming them.
The beta-binomial against the sample's own matched control replaced it and matched it exactly (64/64 confirmed retained, 9 rejected, precision 0.877), so the cohort-scale PoN was dead weight: it needed a set of unedited samples that a single-sample submission does not have, and bought nothing the model does not. Deleted: bin/build_offtarget_pon.py modules/local/pon_score.nf (PON_SCORE, BUILD_PON) params review_auto_pon, offtarget_pon, pon_min_reads, pon_min_donors review_filter.py --pon and its rule-4 branch the BUILD_PON publishDir block This also removes a per-control scoring step from every run: PON_SCORE ran once per (control CRAM, target file) pair purely to build a file rule 4 no longer consults. NOT touched, despite the name: PON_OFFTARGET_FILTER / bin/pon_filter.py. That is a different thing entirely -- a pooled cross-donor PoN in the -entry OFFTARGET worklist arm, which pools all unedited normals to recover depth the single matched normal lacks. Nothing here replaces it. The rule-4 guardrail survives the removal and now reads "No --noise-model, and only 1 guide in this invocation", since the PoN is no longer one of the fixes it can suggest. Verified: filter output byte-identical to before the rip-out (96 rows on the 32-sample cohort), the default DAG builds with no PoN processes, the OFFTARGET arm still includes its own filter, 30 tests pass.
…removal The header still quoted the PoN-era result (88/87 queue) and the review_noise_model 'off' option still described keeping a boolean PoN that no longer exists. Replaced with the current measurement (96 queue, precision 0.889, 81/81 on-target, 70/70 ECS-confirmed incl. 12/12 below 5% VAF) and what 'off' actually does now, which is fall back to cross-guide recurrence.
README.md advertised an offtarget_pon.tsv output that is no longer produced and described a panel of normals the pipeline no longer builds. docs/OFFTARGET.md documented review_auto_pon as "on by default", listed removed params in its parameter table, showed --pon in the standalone example, and still quoted the old 10/0.05 entry gate. Replaced the "Panel of normals" section with "The noise model": what AQ is, why review_aq_min defaults to 5, what the depth floor protects against, that germline falls out for free, and what happens on a single guide with no model. docs/THRESHOLD_CHANGE_IMPACT.md is a historical measurement and stays, but two of its three recommendations were about tuning pon_min_reads and unioning the PoN with cross-guide recurrence -- neither of which exists now. Added a banner saying Finding 1 stands and recommendations 1 and 3 must not be acted on. Docs only: no bin/ changes, so the in-flight run's cache is untouched.
Rule 4 became a control-based statistical test, which makes rule 1 (matched control carries the indel at >=5% VAF) look like the crude version of the same idea. Measured, it is not: disabling rule 1 costs precision 0.889 -> 0.877 and admits another human-rejected site, while changing recall not at all. Of the 221 gated rows with control VAF in [5%,30%) -- all of which rule 1 drops -- rule 4 keeps 26, other rules catch 23, and 3 reach the queue. The clearest is CHD2-KO-DNA chr20:60,991,255: 26.0% in the edited sample against 22.6% in its own matched control, near-identical fractions, and the beta-binomial still passes it at AQ 5.33. That is structural. The binomial asks "is k surprising given p?", and at 77 reads a modest excess over a high background clears a low bar; rule 1 asks "does the control carry this at all?", which is the right question for germline. Raising aq_min cannot absorb it: catching all three needs aq_min > 17.13, and by 10 a confirmed edit and an on-target site are already lost, by 18 two on-target sites and an ECS-confirmed edit. Those germline sites rank ABOVE genuine edits in AQ order, so no cut separates them. Also documents a parsing trap in --keep-all output: why_dropped is empty for surviving rows and pandas reads that as NaN, so filtering with == "" silently returns nothing. It under-counted this analysis on the first pass.
… tests
bin/bnd_from_indels_to_vcf.py wrote ZERO records in every run ever performed.
Doubled escapes ('\\t', '\\n') meant the header never split, so the
'bnd_count' in hdr test was always False and every row hit continue. All 32
*.offtarget_svs.vcf were a single line of literal \n text. Fixed, and the
32-sample CAR-T cohort now emits 1022 breakend records (3-95 per sample).
Nothing was ever miscalled -- only the VCF serialisation step was broken, so
the caller does not need re-running.
This does NOT restore the report's SV section, and an earlier draft of the
docs claimed it would. COMPILE_REPORT_JSON does consume the VCF and
tables.bnd_vcf is now populated 32/32, but make_scge_report.qmd never
references bnd_vcf and neither does make_scge_excel.py. The report's SV panel
renders tables.on_target_sv_transgene, an unrelated transgene-junction table
that is empty in all 32 samples. Two independent defects; only the VCF one is
fixed here. docs/OFFTARGET.md says so explicitly.
pileup_snapshot.py was discarding the only evidence a junction has: it
filtered is_supplementary unconditionally, i.e. exactly the half of every
split read that proves a breakend. Supplementary retention is now opt-in via
keep_supplementary=, so indel snapshots are unchanged, and bnd_snapshots.py
renders one figure per JUNCTION rather than per row -- the caller reports each
junction from both ends at +/-4bp jitter under both orientations, so the
cohort's 25 queue rows are 8 junctions. A queue row is not an event, and a
row's read count is not the junction's support.
compile_report_json.nf had no stub: block, so -stub-run died there and had
never once reached the review arm. Added.
Tests now cover this path, which nothing under tests/ previously touched:
test_bnd_vcf.py (14 cases incl. the literal-escape regression),
test_bnd_snapshots.py (17 cases against the real 25-row queue), and
review_bnd_snapshots.nf.test (3 container cases against a 48 KB synthetic
split-read CRAM). nf-test 8/8; pytest 55 passed, 6 pre-existing edlib
ModuleNotFoundError failures in test_offtarget_glue.py.
…as filters Answers four questions put to the noise model: where the background comes from and how AQ is calculated, whether the DRAGEN SNV panel can filter indels, the same for the SV panels and breakends, and whether distance from the cut can be modelled rather than thresholded. docs/NOISE_MODEL.md derives the model: a Binomial sampling process with a Beta mixing distribution over the per-locus rate, giving a Beta-Binomial whose prior is fitted once per run by method of moments and updated per locus to alpha = a0 + control_alt, beta = b0 + control_ref. AQ is stated as what it is -- a survival function in Phred units, P(X >= k | H0) -- and explicitly NOT a likelihood ratio; the genuine likelihood ratio in this work is the spatial one, P(d|signal)/P(d|background), one section later. docs/NOISE_MODEL_VALIDATION.md tests whether that model describes the data. The two-process objection is correct and confirmed: 99.40% of 95,073 control observations sit at exactly zero, a second mode sits at 0.179, and there are zero observations in between, so a single Beta is fitted across a population that does not exist. It is nevertheless statistically inert here. Under a genuine null -- each donor's control scored against the other donors' at the same locus, 71,755 observations -- the shipped prior is calibrated at KS D ~ 0.004 against 0.140 for a plain Binomial, a 35x gap stable across seeds. The reason is that a < 1 makes the Beta monotone decreasing, so its shape already mimics zero inflation; the explicit mixture is unidentifiable (pi swings 0.536 -> 0.993 while held-out likelihood moves 0.000002 nats) and its epsilon collapses to the denormal floor, so the depth floor stays. No model change is recommended: every Beta-Binomial variant drops the same single row. Because the calibration test uses randomised p-values, its p-value is one draw. --calib-seeds sweeps it (default 20) and the docs quote D, not p: over 20 seeds the shipped model runs D 0.0028-0.0062 but p 0.009-0.624, with one seed rejecting at 0.05. randomised_tails() returns the seed-independent P(X > k) and P(X = k) so sweeping only repeats the cheap uniform draw. All models are scored on shared draws; an earlier unpaired comparison made MML look better calibrated than MOM (D 0.0025 vs 0.0047) and that gap was an artifact of the random state, not a real difference. The real gap is sequence context: sites in a >=9bp homopolymer carry control indel evidence 35.4% of the time against 0.53% in unique sequence. It costs nothing today -- no queue row sits in one -- which is an argument for fixing it before it costs something. docs/PANEL_AS_FILTER.md measures the panels. 152 of 99,238 indel calls (0.15%) touch an indel-capable noisy SNV locus and only 19 clear the shipped >=3-donor bar; the review queue is panel-clean. As a standalone filter the panel is a specificity instrument, flagging 1 of 69 confirmed edits but reaching sensitivity 0.50 at best, and it cannot stand in for the matched control: 64% of the germline calls rule 1 catches are invisible to it, because a population panel structurally cannot represent one donor's genotype. On the SV side the mechanism is simpler than previously documented -- IDPF v3.0.0's intervals cover 62.66% of hg38, so it flags 83% of all junctions and 25/25 real ones. That is blanket coverage, not population mismatch, and no min_donors threshold rescues it (all 25 hits carry >=3 donors). Only WGS_hg38_v3.1.0, at 1.09% coverage, discriminates. bin/panel_overlap.py restores a lost capability: the three-panel comparison published in CALLER_INTEGRATION_LOG.md had no surviving generator, so nothing could be re-measured when the caller changed underneath it. Its "11/11 artifacts" column is not reproducible on any run still on disk and should be re-derived rather than quoted. Corrects three stale in-repo claims, including review_filter_bnd.py's assertion that the caller caps cut distance at 25bp -- the -d default is 10 and get_indels.nf passes no override, which is why rule 2 drops zero rows. Analysis only. Nothing the pipeline computes changes.
The committed figure was rendered before randomised_tails/draw_pvalues split the seed-independent tail from the uniform draw, so each model consumed its own draws from a shared RNG stream rather than being scored on shared ones. Redrawn on paired draws, the three Beta-Binomial curves lie on top of each other along the diagonal, which is what the text and the seed sweep already said: the variants are indistinguishable by calibration as well as by held-out likelihood. The Binomial's departure is unchanged. Figure only; no numbers in the docs move. noise_two_populations.png is byte-identical, as it carries no RNG dependence.
…which The panel's truth-label evaluation joined against all 99,238 site-rows, while NOISE_MODEL.md's cut-distance work joins against the 479 gated rows. Both are defensible and both were correct, but neither document named its population, so the two sets of figures invited being read against each other. panel_overlap.py now reports both, and PANEL_AS_FILTER.md prints both tables with a note on which is operational: review_filter.py only ever runs this rule on gated rows, so that is the number describing what the panel does in the pipeline, while the all-rows figure answers the broader question of whether it could replace the noise model outright. The gap is not small. On gated rows the panel reaches sensitivity 0.732 at the loosest definition and 0.437 indel-capable, against 0.495 and 0.171 on all site-rows -- roughly 1.5x. That is the denominator changing, not the panel improving: gated rows are enriched for the recurrent-locus artifacts a population panel can see, while the full table is dominated by sites carrying no indel evidence at all. The base rate moves with it, from 75% negative to 53%, so both tables now print what a do-nothing filter would score. Every conclusion survives either denominator, which is the reason for showing both: exactly one confirmed edit is flagged in each population, both indel-aware definitions flag zero in each, and the rule that ships catches 0.033 / 0.042 of artifacts and claims no queue rows at all.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
PR checklist
nf-core lint).nextflow run . -profile test,docker --outdir <OUTDIR>).nextflow run . -profile debug,test,docker --outdir <OUTDIR>).docs/usage.mdis updated.docs/output.mdis updated.CHANGELOG.mdis updated.README.mdis updated (including new tool citations and authors/contributors).