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).

Co-Authored-By: Claude Opus 5.5 <[email protected]>
Claude-Session: https://claude.ai/code/session_01814GKVzZVL3FuF3DeiZXhv

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

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())