From d5ab45886c705a99f89f72d900cfadac47911fec Mon Sep 17 00:00:00 2001 From: Kemal Yaylali Date: Tue, 29 Sep 2026 11:57:58 +0100 Subject: [PATCH] cn_validate_truth: sensitivity by carrier frequency (rare < 5 % is the primary metric, K-027) --- scripts/v2/cn_tune_report.py | 2 ++ scripts/v2/cn_validate_truth.py | 20 +++++++++++++++++++- 2 files changed, 21 insertions(+), 1 deletion(-) diff --git a/scripts/v2/cn_tune_report.py b/scripts/v2/cn_tune_report.py index 0fa2d59..2f8b9fa 100644 --- a/scripts/v2/cn_tune_report.py +++ b/scripts/v2/cn_tune_report.py @@ -55,6 +55,8 @@ def main() -> int: 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()))) + for t, bins in val.get("sensitivity_by_carrier_frequency", {}).items(): + print(f" by frequency {t}: " + ", ".join(f"{b} {x['found']}/{x['truth']}" for b, x in sorted(bins.items()))) return 0 diff --git a/scripts/v2/cn_validate_truth.py b/scripts/v2/cn_validate_truth.py index d6429c3..be10bb1 100644 --- a/scripts/v2/cn_validate_truth.py +++ b/scripts/v2/cn_validate_truth.py @@ -6,7 +6,9 @@ structural-variant callset (1KGP_3202 GATK-SV/svtools freeze V3) and between rep ``truth.tsv``: ``chrom start end svtype sample`` for non-reference DEL/DUP genotypes of the 120 individuals (``bcftools query``; chr prefix removed). Matching: reciprocal overlap ≥ 50 % and the same class (DEL ↔ losses -incl. mosaic, DUP ↔ gains). Only truth events with ≥ 10 array probes are counted as detectable. +incl. mosaic, DUP ↔ gains). Only truth events with ≥ 10 array probes are counted as detectable. Sensitivity is also +split by the event's carrier frequency among the individuals in ``truth.tsv`` (identical chrom/start/end/type): +common events are normalised away by a median-LRR batch reference (K-027), so ``rare`` (< 5 %) is the primary metric. """ from __future__ import annotations @@ -22,6 +24,7 @@ import numpy as np from simurg.formats.bpm import read_bpm +FREQ = [(0.0, 0.05, "rare<5%"), (0.05, 0.2, "5-20%"), (0.2, 1.01, "common>=20%")] BINS = [(25_000, 50_000), (50_000, 100_000), (100_000, 500_000), (500_000, 1_000_000), (1_000_000, 10**9)] CLASS = {"DEL": "DEL", "MDEL": "DEL", "DUP": "DUP", "MDUP": "DUP"} @@ -70,15 +73,23 @@ def main() -> int: ign = set(a.ignore_filter) | {"PASS"} calls[s["sample_id"]] = [e for e in s["events"] if e["quality"] >= a.min_qual and set(e["filters"]) <= ign] truth: dict[str, list[tuple[str, int, int, str]]] = defaultdict(list) + carriers: dict[tuple[str, int, int, str], set[str]] = defaultdict(set) with a.truth.open() as fh: for line in fh: c, s_, e_, t, ind = line.rstrip("\n").split("\t") c = c[3:] if c.startswith("chr") else c + carriers[(c, int(s_), int(e_), t)].add(ind) if t in ("DEL", "DUP") and int(e_) - int(s_) + 1 >= 25_000: truth[ind].append((c, int(s_), int(e_), t)) + n_ind = len({i for v in carriers.values() for i in v}) or 1 + + def freq_class(key: tuple[str, int, int, str]) -> str: + f = len(carriers[key]) / n_ind + return next(name for lo, hi, name in FREQ if lo <= f < hi) # sensitivity per class and size bin sens: dict[str, dict[str, list[int]]] = defaultdict(lambda: defaultdict(lambda: [0, 0])) + sens_freq: dict[str, dict[str, list[int]]] = defaultdict(lambda: defaultdict(lambda: [0, 0])) prec = defaultdict(lambda: [0, 0]) per_sample = [] for sid, evs in calls.items(): @@ -93,6 +104,9 @@ def main() -> int: cell = sens[t][bin_of(e_ - s_ + 1)] cell[0] += ok cell[1] += 1 + fc = sens_freq[t][freq_class((c, s_, e_, t))] + fc[0] += ok + fc[1] += 1 hit_t += ok hit_c = 0 cnv_calls = [e for e in evs if e["call"] in CLASS] @@ -149,6 +163,10 @@ def main() -> int: t: {b: {"found": v[0], "truth": v[1], "rate": rate(v)} for b, v in sorted(d.items())} for t, d in sens.items() }, + "sensitivity_by_carrier_frequency": { + t: {b: {"found": v[0], "truth": v[1], "rate": rate(v)} for b, v in sorted(d.items())} + for t, d in sens_freq.items() + }, "precision": {t: {"in_truth": v[0], "calls": v[1], "rate": rate(v)} for t, v in prec.items()}, "replicate_reproducibility": { k: {"matched": v[0], "comparisons": v[1], "rate": rate(v)} for k, v in rep.items()