From 803f4aaec31d7710b43794670a9822f7fe6eb2ee Mon Sep 17 00:00:00 2001 From: Kemal Yaylali Date: Tue, 29 Sep 2026 10:56:21 +0100 Subject: [PATCH] CN: MF split normalised by per-probe reference BAF noise; genotype-based complete loss (CYTO-041) The mosaic BAF split is now estimated over heterozygous calls against each probe's expected noise (the reference heterozygous spread scaled to the sample), inverting the median of |N(+-d, sigma)|. A mosaic loss lacking the heterozygous calls its PFB predicts gets MF 1 and is promoted (CYTO-020). --- SIMURGARRAY_FUNCTIONAL_SPEC.md | 16 ++++---- src/simurg/cn/events.py | 67 +++++++++++++++++++++++++--------- src/simurg/pipeline/cn.py | 23 +++++++++++- tests/unit/test_cn.py | 15 +++++--- 4 files changed, 90 insertions(+), 31 deletions(-) diff --git a/SIMURGARRAY_FUNCTIONAL_SPEC.md b/SIMURGARRAY_FUNCTIONAL_SPEC.md index af25e14..4e47c04 100644 --- a/SIMURGARRAY_FUNCTIONAL_SPEC.md +++ b/SIMURGARRAY_FUNCTIONAL_SPEC.md @@ -2437,13 +2437,15 @@ about **algorithms or numerics**, the plan wins. Each such conflict found during - **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. + - *noise-corrected split:* per-probe expected noise σ_i = the reference's heterozygous BAF spread of the probe × + the sample's scale (median |BAF − 0.5| / spread over autosomal heterozygous calls, ÷ 0.6745; the sample's + heterozygous BAF SD without a reference spread). Over the event's heterozygous calls, m = median(|BAF − 0.5| / + σ_i); the split in σ units δ solves Φ(m − δ) − Φ(−m − δ) = ½ (δ = 0 when m ≤ 0.6745) and d = median(σ_i)·δ. + Without this, noise alone gave MF ≈ 0.1–0.3 and genomic waves were called as mosaic losses. + - *complete loss:* a mosaic loss whose called genotypes hold fewer than a quarter of the heterozygotes the PFB + predicts (Σ 2p(1−p) ≥ 3) 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. --- diff --git a/src/simurg/cn/events.py b/src/simurg/cn/events.py index af1738e..3335939 100644 --- a/src/simurg/cn/events.py +++ b/src/simurg/cn/events.py @@ -20,8 +20,19 @@ 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 +LOSS_EXPECTED_HETS = 3.0 # CYTO-041: expected heterozygous calls needed to call a loss complete +LOSS_HET_RATIO = 0.25 # CYTO-041: observed / expected heterozygous calls below which a loss is complete + + +def complete_loss(geno: npt.NDArray[np.uint8], pfb: F64) -> bool: + """CYTO-041 [OURS]: a loss is complete when its called genotypes lack the heterozygotes the PFB predicts + (Σ 2p(1 − p) ≥ 3 expected, fewer than a quarter observed); a mosaic loss keeps its heterozygotes.""" + called = geno > 0 + p = np.clip(np.where(np.isfinite(pfb), pfb, 0.5), 0.0, 1.0)[called] + expected = float((2 * p * (1 - p)).sum()) + return expected >= LOSS_EXPECTED_HETS and float((geno[called] == 2).sum()) < LOSS_HET_RATIO * expected + + 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 @@ -117,21 +128,39 @@ def extract(chrom: str, path: npt.NDArray[np.int32], states: list[State], pos: I return out -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 +def _phi(x: float) -> float: + return 0.5 * (1.0 + math.erf(x / math.sqrt(2.0))) + + +def _split_from_median(m: float) -> float: + """δ ≥ 0 with median |N(±δ, 1)| = m (bisection; 0 when m ≤ 0.6745).""" + if _phi(m) - _phi(-m) <= 0.5: + return 0.0 + lo, hi = 0.0, m + 10.0 + for _ in range(60): + mid = 0.5 * (lo + hi) + if _phi(m - mid) - _phi(-m - mid) > 0.5: + lo = mid + else: + hi = mid + return 0.5 * (lo + hi) + + +def _mf_from_baf(call: str, baf: F64, sigma: F64 | float = 0.0) -> float | None: + """CYTO-018 with the CYTO-041 noise correction [OURS]. ``baf``: heterozygous-call BAF of the event; ``sigma``: + each probe's expected BAF noise (reference heterozygous spread × the sample's scale). With m = median(|BAF − 0.5| + / σ), the split δ (in σ units) solves Φ(m − δ) − Φ(−m − δ) = ½ (the median of |N(±δ, 1)|; δ = 0 when + m ≤ 0.6745, i.e. no more spread than noise), and d = median(σ)·δ.""" + sig = np.broadcast_to(np.asarray(sigma, dtype=np.float64), baf.shape) + ok = np.isfinite(baf) & (baf > 0.1) & (baf < 0.9) & np.isfinite(sig) + b, sg = baf[ok], sig[ok] if len(b) < MIN_MF_PROBES: return None - 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 float(np.median(sg)) > 0: + sg = np.maximum(sg, 1e-3) + d = float(np.median(sg)) * _split_from_median(float(np.median(np.abs(b - 0.5) / sg))) + else: + d = float(np.median(np.abs(b - 0.5))) if call == "MDEL": f = 2 * d / (0.5 + d) elif call == "MDUP": @@ -151,7 +180,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, baf_sd: float = 0.0, pfb: F64 | None = None, + gamma: float, haploid: bool, sigma: F64 | None = None, pfb: F64 | None = None, ) -> None: # fmt: skip """CYTO-010/012/018/027 statistics in place (``lrr`` unsmoothed).""" sl = slice(ev.i0, ev.i1) @@ -170,7 +199,11 @@ 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], baf_sd, None if pfb is None else pfb[sl]) + if not haploid and ev.call == "MDEL" and pfb is not None and complete_loss(g, pfb[sl]): + ev.mosaic_fraction = 1.0 + return + hb = np.where(g == 2, baf[sl], np.nan) if sigma is not None else baf[sl] + mf = None if haploid else _mf_from_baf(ev.call, hb, 0.0 if sigma is None else sigma[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 4d87904..c4a7ea4 100644 --- a/src/simurg/pipeline/cn.py +++ b/src/simurg/pipeline/cn.py @@ -288,6 +288,26 @@ def _resolve_sex(si: SampleInput, chroms: npt.NDArray[np.str_], lrr: F64, geno: return "U", "unresolved", ["sex_unresolved"] +def _baf_noise( + baf: npt.NDArray[np.float64], geno: npt.NDArray[np.uint8], auto: npt.NDArray[np.bool_], + ref_het_sd: npt.NDArray[np.float32] | None, baf_sd: float, +) -> npt.NDArray[np.float64]: # fmt: skip + """CYTO-041: expected heterozygous BAF noise per probe, the reference spread scaled to this sample + (median |BAF − 0.5| / spread over autosomal heterozygous calls = 0.6745 × scale); ``baf_sd`` without a reference + spread or where a probe has none.""" + out = np.full(len(baf), baf_sd, dtype=np.float64) + if ref_het_sd is None: + return out + h = np.asarray(ref_het_sd, dtype=np.float64) + hc = auto & (geno == 2) & np.isfinite(baf) & np.isfinite(h) & (h > 0) + if hc.sum() < 100: + return out + scale = float(np.median(np.abs(baf[hc] - 0.5) / np.maximum(h[hc], 1e-3))) / 0.6745 + ok = np.isfinite(h) & (h > 0) + out[ok] = np.maximum(h[ok], 1e-3) * scale + return out + + def call_sample(si: SampleInput, ctx: dict[str, Any], opts: CnCallOptions) -> dict[str, Any]: """One sample: preprocessing, per-region HMM, events, VCF + segments. Returns the summary row + events.""" from simurg.cn import hmm @@ -346,6 +366,7 @@ def call_sample(si: SampleInput, ctx: dict[str, Any], opts: CnCallOptions) -> di seg_sd = _robust_sd(lrr[auto], 0.05) b_het = baf[auto & (baf > 0.3) & (baf < 0.7)] baf_sd = _robust_sd(b_het, 0.02) + sigma = _baf_noise(baf, geno, auto, model.ref_het_sd, baf_sd) # CYTO-041 dip, hap = hmm.diploid_states(opts.mosaic_prior), hmm.haploid_states(opts.mosaic_prior) counts: dict[str, int] = {} events: list[Any] = [] @@ -392,7 +413,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, baf_sd, pfb[idx]) + annotate_stats(ev, p, x_raw, b, gt, post, states, m_sp, ctx["gamma"], haploid, sigma[idx], 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 87a951d..7e0494e 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, _mf_from_baf, apply_thresholds, effective_size, extract +from simurg.cn.events import Event, Thresholds, _mf_from_baf, apply_thresholds, complete_loss, 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 @@ -237,12 +237,15 @@ def test_mf_refinements() -> None: # CYTO-018 / CYTO-041 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 + assert (_mf_from_baf("MDEL", b) 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 + # no heterozygous calls where ~12 are expected: a complete (constitutional) loss + hom = np.full(30, 1, dtype=np.uint8) + assert complete_loss(hom, np.full(30, 0.3)) + assert not complete_loss(hom, np.full(30, 0.02)) # too few expected + mixed = hom.copy() + mixed[:8] = 2 + assert not complete_loss(mixed, np.full(30, 0.3))