From 6bc42e53b79e7bc94b8e55cad510a7e79d331ce2 Mon Sep 17 00:00:00 2001 From: Kemal Yaylali Date: Tue, 29 Sep 2026 10:30:12 +0100 Subject: [PATCH] CN: noise-corrected BAF split and complete-loss MF (CYTO-041); Pi tuning harness Mosaic fraction from the BAF split now subtracts the heterozygous BAF noise (median |BAF-0.5| of unsplit hets = 0.6745 sigma_b), and a mosaic loss lacking its expected heterozygotes gets MF 1 and is promoted to a constitutional DEL (CYTO-020). Adds scripts/v2/cn_tune.sh + cn_tune_report.py (GSE165845 rounds on the Pi). --- SIMURGARRAY_FUNCTIONAL_SPEC.md | 11 ++++++ scripts/v2/cn_tune.sh | 29 ++++++++++++++++ scripts/v2/cn_tune_report.py | 62 ++++++++++++++++++++++++++++++++++ src/simurg/cn/events.py | 24 +++++++++---- src/simurg/pipeline/cn.py | 2 +- tests/unit/test_cn.py | 20 ++++++++++- 6 files changed, 140 insertions(+), 8 deletions(-) create mode 100755 scripts/v2/cn_tune.sh create mode 100644 scripts/v2/cn_tune_report.py diff --git a/SIMURGARRAY_FUNCTIONAL_SPEC.md b/SIMURGARRAY_FUNCTIONAL_SPEC.md index 03f33fb..af25e14 100644 --- a/SIMURGARRAY_FUNCTIONAL_SPEC.md +++ b/SIMURGARRAY_FUNCTIONAL_SPEC.md @@ -2435,6 +2435,16 @@ about **algorithms or numerics**, the plan wins. Each such conflict found during Each yields valid, possibly empty, VCF and JSON output (EDGE-024). (src: E-35) +- **CYTO-041** `§10.3 · v2` **MF refinements [OURS].** Two refinements of CYTO-018, from real-data validation + (docs/validation/V2_cn.md, K-027): + - *noise-corrected split:* the median |BAF − 0.5| of unsplit heterozygotes is 0.6745·σ_b (σ_b = the sample's + heterozygous BAF SD), so `d = √max(0, d_obs² − (0.6745·σ_b)²)`. Without this, noise alone gave MF ≈ 0.2 in + diploid regions and waves were called as mosaic losses. + - *complete loss:* a mosaic loss with ≥ 20 finite-BAF probes, ≥ 5 expected heterozygotes (Σ 2p(1−p) over PFB p), + and fewer than 0.25 × expected probes with 0.1 < BAF < 0.9 has MF = 1. With the default + `--max-mosaic-fraction 1.0` it is therefore promoted (CYTO-020) to a constitutional DEL with `HIGHFRACTION`. + This catches constitutional deletions whose LRR is compressed towards the mosaic states. + --- ## 7. Targeted PGx copy number (v2) @@ -3710,6 +3720,7 @@ Generated from the item tags (`build.py`). Phases: P0–P10 as in plan §18; `v2 | CYTO-038 | Model building | §10 | v2 | A-G52 | | CYTO-039 | Batch scale | §10, §19 | v2 | D-9, E-35 | | CYTO-040 | Degenerate samples | §10 | v2 | E-35 | +| CYTO-041 | MF refinements | §10.3 | v2 | K-027 | | PCN-001 | Region catalogue | §11 | v2 | B-23, C-G25, D-10 | | PCN-002 | Genotype every region [EMU] | §11 | v2 | B-23 | | PCN-003 | Inputs | §11 | v2 | A-G37 | diff --git a/scripts/v2/cn_tune.sh b/scripts/v2/cn_tune.sh new file mode 100755 index 0000000..3340ade --- /dev/null +++ b/scripts/v2/cn_tune.sh @@ -0,0 +1,29 @@ +#!/usr/bin/env bash +# CN tuning round on the Pi (GSE165845): reference = replicate-3 arrays of all 120 individuals, test = replicates 1 +# and 2 of the first N individuals; truth = NYGC SV rows (data/v2cn/cloud1/truth.tsv). +# bash scripts/v2/cn_tune.sh TAG [extra `simurg cn call` options...] (env: N=20, MODEL=...) +set -euo pipefail +cd "$(dirname "$0")/../.." +TAG=${1:?tag}; shift +N=${N:-20} +G=data/v2cn/gse_gt; M=data/v2cn/gse_idats/samples_map.tsv; W=data/v2cn/tune +BPM=data/p8/gsa1/GSA-24v1-0_C2.bpm +mkdir -p $W/ref $W/test +if [[ ! -f $W/.dirs ]]; then + inds=$(tail -n +2 $M | cut -f3 | awk '!s[$0]++' | head -n "$N") + tail -n +2 $M | while IFS=$'\t' read -r sid gsm ind rep; do + [[ $rep == 3 ]] && ln -sfn "$(realpath $G/samples/$sid)" $W/ref/$sid + if [[ $rep != 3 ]] && grep -qx "$ind" <<< "$inds"; then ln -sfn "$(realpath $G/samples/$sid)" $W/test/$sid; fi + done + touch $W/.dirs +fi +MODEL=${MODEL:-$W/ref.simurg-cyto} +if [[ ! -f $MODEL ]]; then + MEM=${MEM:-5G} MIN_FREE_MB=1500 bash scripts/pi_run.sh .venv/bin/simurg cn build-model --bpm $BPM --fasta data/p6/GRCh38.primary.fa \ + --egt data/p8/gsa1/GSA-24v1-0_C1_ClusterFile.egt --reference-gtcs $W/ref --max-reference 120 --out $MODEL --quiet +fi +MEM=${MEM:-5G} MIN_FREE_MB=1500 bash scripts/pi_run.sh .venv/bin/simurg cn call --bpm $BPM --model $MODEL --gtcs $W/test \ + --summary $G/sample_summary.tsv --out $W/$TAG --threads 3 --force --quiet "$@" +.venv/bin/python scripts/v2/cn_validate_truth.py --calls $W/$TAG --map $M --truth data/v2cn/cloud1/truth.tsv --bpm $BPM \ + --min-qual 20 --out $W/$TAG/val > /dev/null +.venv/bin/python scripts/v2/cn_tune_report.py $W/$TAG $M diff --git a/scripts/v2/cn_tune_report.py b/scripts/v2/cn_tune_report.py new file mode 100644 index 0000000..0fa2d59 --- /dev/null +++ b/scripts/v2/cn_tune_report.py @@ -0,0 +1,62 @@ +#!/usr/bin/env python3 +"""One-line-per-call-type summary of a tuning round: calls, fraction in NYGC truth, fraction seen in the other +replicate; plus truth sensitivity by size bin (QUAL >= 20, PASS or SampleQuality).""" + +from __future__ import annotations + +import csv +import json +import sys +from collections import defaultdict +from pathlib import Path + + +def ro(a: tuple[int, int], b: tuple[int, int]) -> float: + ov = min(a[1], b[1]) - max(a[0], b[0]) + 1 + return 0.0 if ov <= 0 else min(ov / (a[1] - a[0] + 1), ov / (b[1] - b[0] + 1)) + + +def main() -> int: + run, mapf = Path(sys.argv[1]), Path(sys.argv[2]) + with mapf.open() as fh: + smap = {r["sample_id"]: r["individual"] for r in csv.DictReader(fh, delimiter="\t")} + truth: dict[str, list[tuple[str, int, int, str]]] = defaultdict(list) + with open("data/v2cn/cloud1/truth.tsv") as fh: + for line in fh: + c, s, e, t, i = line.rstrip("\n").split("\t") + truth[i].append((c.removeprefix("chr"), int(s), int(e), t)) + side = json.loads((run / "cn_call_manifest.json").read_text()) + calls = {s["sample_id"]: [e for e in s["events"] if e["quality"] >= 20 and set(e["filters"]) <= {"PASS", + "SampleQuality"}] for s in side["samples"] if s.get("status") == "OK"} # fmt: skip + cls = {"DEL": "DEL", "MDEL": "DEL", "DUP": "DUP", "MDUP": "DUP"} + by_ind: dict[str, list[str]] = defaultdict(list) + for sid in calls: + by_ind[smap[sid]].append(sid) + st: dict[str, list[int]] = defaultdict(lambda: [0, 0, 0, 0]) + for sid, evs in calls.items(): + ind = smap[sid] + others = [calls[x] for x in by_ind[ind] if x != sid] + for e in evs: + k = e["call"] + v = st[k] + v[0] += 1 + if k in cls: + v[1] += any(t == cls[k] and c == e["chrom"] and ro((s, en), (e["start"], e["end"])) >= 0.5 + for c, s, en, t in truth[ind]) # fmt: skip + for oc in others: + v[2] += any(o["chrom"] == e["chrom"] and o["call"] == k and ro((o["start"], o["end"]), + (e["start"], e["end"])) >= 0.5 for o in oc) # fmt: skip + v[3] += 1 + print(f"{run.name}: {len(calls)} samples") + for k, v in sorted(st.items()): + tr = f"{v[1] / v[0]:.2f}" if k in cls and v[0] else " - " + rp = f"{v[2] / v[3]:.2f}" if v[3] else " - " + print(f" {k:8s} calls {v[0]:5d} in-truth {tr} replicated {rp}") + val = json.loads((run / "val" / "summary.json").read_text()) + for t, bins in val["sensitivity"].items(): + print(f" sensitivity {t}: " + ", ".join(f"{b} {x['found']}/{x['truth']}" for b, x in sorted(bins.items()))) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/src/simurg/cn/events.py b/src/simurg/cn/events.py index 3d93454..af1738e 100644 --- a/src/simurg/cn/events.py +++ b/src/simurg/cn/events.py @@ -2,7 +2,7 @@ From a per-chromosome Viterbi path: runs of non-baseline states become events; statistics per CYTO-027; effective size per CYTO-010 (gaps > 150 × the chromosome's median spacing replaced by the median); QUAL from the mean posterior -of the called state class (CYTO-012); mosaic fraction from the BAF split or, failing that, LRR (CYTO-018/019); +of the called state class (CYTO-012); mosaic fraction from the BAF split or, failing that, LRR (CYTO-018/019/042); post-filters in the order of CYTO-015; per-type thresholds (CYTO-007/008) and promotion/low fraction (CYTO-020/021). """ @@ -20,6 +20,8 @@ F64 = npt.NDArray[np.float64] I64 = npt.NDArray[np.int64] GAP_FACTOR = 150 MIN_MF_PROBES = 20 +LOH_EXPECTED_HETS = 5.0 # CYTO-041: expected heterozygous probes needed to call a loss complete +LOH_HET_RATIO = 0.25 # CYTO-041: observed / expected heterozygotes below which a loss is complete TYPES = ("del", "dup", "loh", "gainloh", "mdel", "mdup", "mloh") DEFAULTS = {"del": (10, 25_000), "dup": (10, 50_000), "loh": (500, 3_000_000), "gainloh": (300, 500_000), "mdel": (10, 25_000), "mdup": (10, 50_000), "mloh": (500, 3_000_000)} # fmt: skip @@ -115,11 +117,21 @@ def extract(chrom: str, path: npt.NDArray[np.int32], states: list[State], pos: I return out -def _mf_from_baf(call: str, baf: F64) -> float | None: - b = baf[np.isfinite(baf) & (baf > 0.1) & (baf < 0.9)] +def _mf_from_baf(call: str, baf: F64, baf_sd: float = 0.0, pfb: F64 | None = None) -> float | None: + """CYTO-018 with the CYTO-041 refinements [OURS]: the split is corrected for BAF noise (the median |BAF − 0.5| of + unsplit heterozygotes is 0.6745·σ_b), and a loss whose expected heterozygotes are absent is complete (MF 1).""" + fin = np.isfinite(baf) + mid = fin & (baf > 0.1) & (baf < 0.9) + b = baf[mid] + if call == "MDEL" and pfb is not None and fin.sum() >= MIN_MF_PROBES: + p = np.clip(np.where(np.isfinite(pfb), pfb, 0.5), 0.0, 1.0)[fin] + expected = float((2 * p * (1 - p)).sum()) + if expected >= LOH_EXPECTED_HETS and len(b) < LOH_HET_RATIO * expected: + return 1.0 if len(b) < MIN_MF_PROBES: return None - d = float(np.median(np.abs(b - 0.5))) + d_obs = float(np.median(np.abs(b - 0.5))) + d = math.sqrt(max(0.0, d_obs**2 - (0.6745 * baf_sd) ** 2)) if call == "MDEL": f = 2 * d / (0.5 + d) elif call == "MDUP": @@ -139,7 +151,7 @@ def _mf_from_lrr(call: str, lrr: F64, gamma: float) -> float | None: def annotate_stats( ev: Event, pos: I64, lrr: F64, baf: F64, geno: npt.NDArray[np.uint8], post: F64, states: list[State], m: float, - gamma: float, haploid: bool, + gamma: float, haploid: bool, baf_sd: float = 0.0, pfb: F64 | None = None, ) -> None: # fmt: skip """CYTO-010/012/018/027 statistics in place (``lrr`` unsmoothed).""" sl = slice(ev.i0, ev.i1) @@ -158,7 +170,7 @@ def annotate_stats( p = float(post[sl][:, cls].sum(axis=1).mean()) if cls else 0.0 ev.quality = 60.0 if p >= 1.0 else float(min(60, max(0, round(-10 * math.log10(max(1e-300, 1 - p)))))) if ev.mosaic: - mf = None if haploid else _mf_from_baf(ev.call, baf[sl]) + mf = None if haploid else _mf_from_baf(ev.call, baf[sl], baf_sd, None if pfb is None else pfb[sl]) ev.mosaic_fraction = mf if mf is not None else _mf_from_lrr(ev.call, x, gamma) diff --git a/src/simurg/pipeline/cn.py b/src/simurg/pipeline/cn.py index 1aa9699..4d87904 100644 --- a/src/simurg/pipeline/cn.py +++ b/src/simurg/pipeline/cn.py @@ -392,7 +392,7 @@ def call_sample(si: SampleInput, ctx: dict[str, Any], opts: CnCallOptions) -> di post = hmm.forward_backward(e, sw, prior) evs = extract(chrom, path, states, p) for ev in evs: - annotate_stats(ev, p, x_raw, b, gt, post, states, m_sp, ctx["gamma"], haploid) + annotate_stats(ev, p, x_raw, b, gt, post, states, m_sp, ctx["gamma"], haploid, baf_sd, pfb[idx]) evs = post_filter(evs, [x_raw[ev.i0 : ev.i1] for ev in evs], th, counts) for s_, e_ in _runs(path): st = states[int(path[s_])] diff --git a/tests/unit/test_cn.py b/tests/unit/test_cn.py index bc1fb00..87a951d 100644 --- a/tests/unit/test_cn.py +++ b/tests/unit/test_cn.py @@ -12,7 +12,7 @@ from typer.testing import CliRunner from simurg.cli.main import app from simurg.cn import hmm -from simurg.cn.events import Event, Thresholds, apply_thresholds, effective_size, extract +from simurg.cn.events import Event, Thresholds, _mf_from_baf, apply_thresholds, effective_size, extract from simurg.cn.gcwave import gc_correct from simurg.cn.model import gc_fractions, population_bfreq from simurg.formats.gtc import Gtc, write_gtc @@ -228,3 +228,21 @@ def test_packaged_cn_data() -> None: # the built-in mask and cytobands ship wit m = read_mask(None, build) assert len(m) >= 24 and all(s <= e for iv in m.values() for s, e in iv) assert gzip.decompress(_builtin_file(f"cytoBand_{build}.txt.gz")).startswith(b"chr1\t0\t") + + +def test_mf_refinements() -> None: # CYTO-018 / CYTO-041 + rng = np.random.default_rng(3) + sd = 0.04 + # unsplit heterozygotes: noise alone must not look like a mosaic split + b = 0.5 + rng.normal(0, sd, 400) + mf = _mf_from_baf("MDEL", b, sd) + assert mf is not None and mf < 0.08 + assert (_mf_from_baf("MDEL", b, 0.0) or 0) > 0.05 # uncorrected, the old estimate was biased upwards + # a true 40 % mosaic loss is still recovered + d = 0.4 / (2 * (2 - 0.4)) + b = 0.5 + rng.choice([-d, d], 400) + rng.normal(0, sd, 400) + assert abs((_mf_from_baf("MDEL", b, sd) or 0) - 0.4) < 0.08 + # no heterozygotes where ~50 are expected: a complete (constitutional) loss + b = rng.choice([0.0, 1.0], 100) + rng.normal(0, 0.01, 100) + assert _mf_from_baf("MDEL", np.clip(b, 0, 1), sd, np.full(100, 0.3)) == 1.0 + assert _mf_from_baf("MDEL", np.clip(b, 0, 1), sd, np.full(100, 0.01)) is None # too few expected