From 70001a0d6c3d64d85bec85c4e315623df5d487a3 Mon Sep 17 00:00:00 2001 From: Kemal Yaylali Date: Tue, 29 Sep 2026 01:30:07 +0100 Subject: [PATCH] v2 CN: synthetic implant evaluation (CYTO-022); reference GTC discovery follows symlinks 216 implanted events in 22 synthetic GSA samples (4 replicates per kind x fraction x size; reference model from 30 other samples): constitutional DEL/DUP 16/16 at 5-20 Mb, mosaic loss/gain/CN-LOH at 15 % / 5 Mb 12/12. --- docs/validation/V2_cn_synthetic_gsa36.json | 97 ++++++++++++++ scripts/v2/cn_synthetic_eval.py | 147 +++++++++++++++++++++ src/simurg/pipeline/cn.py | 3 +- 3 files changed, 246 insertions(+), 1 deletion(-) create mode 100644 docs/validation/V2_cn_synthetic_gsa36.json create mode 100644 scripts/v2/cn_synthetic_eval.py diff --git a/docs/validation/V2_cn_synthetic_gsa36.json b/docs/validation/V2_cn_synthetic_gsa36.json new file mode 100644 index 0000000..b77533f --- /dev/null +++ b/docs/validation/V2_cn_synthetic_gsa36.json @@ -0,0 +1,97 @@ +{ + "events": 216, + "detected": 126, + "table": { + "MDEL f=0.10": { + "5 Mb": "1/4", + "20 Mb": "3/4", + "1 Mb": "1/4" + }, + "MLOH f=0.10": { + "5 Mb": "0/4", + "20 Mb": "0/4", + "1 Mb": "0/4" + }, + "MDUP f=0.20": { + "1 Mb": "4/4", + "20 Mb": "4/4", + "5 Mb": "3/4" + }, + "MLOH f=0.15": { + "1 Mb": "0/4", + "5 Mb": "4/4", + "20 Mb": "4/4" + }, + "MLOH f=0.30": { + "1 Mb": "0/4", + "20 Mb": "4/4", + "5 Mb": "2/4" + }, + "DEL": { + "1 Mb": "4/4", + "5 Mb": "4/4", + "20 Mb": "4/4" + }, + "MLOH f=0.20": { + "1 Mb": "0/4", + "5 Mb": "4/4", + "20 Mb": "4/4" + }, + "MLOH f=0.05": { + "1 Mb": "0/4", + "5 Mb": "0/4", + "20 Mb": "0/4" + }, + "MDEL f=0.20": { + "1 Mb": "4/4", + "20 Mb": "4/4", + "5 Mb": "3/4" + }, + "MDEL f=0.30": { + "20 Mb": "4/4", + "5 Mb": "4/4", + "1 Mb": "4/4" + }, + "MDUP f=0.10": { + "20 Mb": "2/4", + "1 Mb": "0/4", + "5 Mb": "2/4" + }, + "MDUP f=0.30": { + "5 Mb": "4/4", + "20 Mb": "4/4", + "1 Mb": "4/4" + }, + "DUP": { + "1 Mb": "2/4", + "20 Mb": "4/4", + "5 Mb": "4/4" + }, + "MDEL f=0.05": { + "5 Mb": "0/4", + "20 Mb": "0/4", + "1 Mb": "0/4" + }, + "LOH": { + "5 Mb": "3/4", + "20 Mb": "4/4", + "1 Mb": "0/4" + }, + "MDUP f=0.15": { + "1 Mb": "1/4", + "5 Mb": "4/4", + "20 Mb": "4/4" + }, + "MDEL f=0.15": { + "20 Mb": "4/4", + "5 Mb": "4/4", + "1 Mb": "3/4" + }, + "MDUP f=0.05": { + "5 Mb": "0/4", + "1 Mb": "0/4", + "20 Mb": "0/4" + } + }, + "calls_total": 564 +} \ No newline at end of file diff --git a/scripts/v2/cn_synthetic_eval.py b/scripts/v2/cn_synthetic_eval.py new file mode 100644 index 0000000..3f42fa7 --- /dev/null +++ b/scripts/v2/cn_synthetic_eval.py @@ -0,0 +1,147 @@ +#!/usr/bin/env python3 +"""Sensitivity of `simurg cn call` on implanted events (spec CYTO-022, plan §16.4). + +Real GTCs of normal samples are copied and events are written into their LRR/BAF: + constitutional DEL (CN1: LRR − 0.45, heterozygous BAF → 0/1), DUP (CN3: LRR + 0.3, het → 1/3 or 2/3), LOH (het → 0/1); + mosaic loss / gain / CN-LOH with cell fraction f (LRR shifted by f × the constitutional shift, heterozygous BAF + moved to 0.5 ± d(f), sign random per probe as in unphased data), keeping each probe's original noise. +Every synthetic sample gets several events on different chromosomes. `simurg cn call` then runs on them and an +implanted event counts as detected when a called event of the same class (loss / gain / LOH) overlaps ≥ 50 % of it. + + python scripts/v2/cn_synthetic_eval.py --bpm M.bpm --model MODEL --gtcs BASE_DIR --out OUT [--threads N] +""" + +from __future__ import annotations + +import argparse +import copy +import gzip +import json +import subprocess +import sys +from pathlib import Path + +import numpy as np + +from simurg.formats.bpm import read_bpm +from simurg.formats.discovery import discover_gtcs +from simurg.formats.gtc import read_gtc, write_gtc + +KINDS = [("DEL", 1.0), ("DUP", 1.0), ("LOH", 1.0)] + [ + (k, f) for k in ("MDEL", "MDUP", "MLOH") for f in (0.05, 0.10, 0.15, 0.20, 0.30) +] +SIZES = (1_000_000, 5_000_000, 20_000_000) +CLASS = {"DEL": "loss", "MDEL": "loss", "DUP": "gain", "MDUP": "gain", "LOH": "loh", "MLOH": "loh", "GAINLOH": "gain"} +CHROMS = [str(c) for c in range(1, 21)] # leave 21, 22, X, Y untouched + + +def implant(lrr: np.ndarray, baf: np.ndarray, idx: np.ndarray, kind: str, f: float, rng: np.random.Generator) -> None: + het = (baf[idx] > 0.15) & (baf[idx] < 0.85) + h = idx[het] + noise = baf[h] - 0.5 + sign = rng.choice([-1.0, 1.0], len(h)) + if kind in ("DEL", "LOH"): + if kind == "DEL": + lrr[idx] -= 0.45 + baf[h] = np.clip(np.where(sign > 0, 1.0, 0.0) + noise * 0.2, 0, 1) + elif kind == "DUP": + lrr[idx] += 0.30 + baf[h] = np.clip(np.where(sign > 0, 2 / 3, 1 / 3) + noise, 0, 1) + else: + d = {"MDEL": f / (2 * (2 - f)), "MDUP": f / (2 * (2 + f)), "MLOH": f / 2}[kind] + lrr[idx] += {"MDEL": -0.45 * f, "MDUP": 0.30 * f, "MLOH": 0.0}[kind] + baf[h] = np.clip(0.5 + sign * d + noise, 0, 1) + + +def main() -> int: + ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("--bpm", required=True) + ap.add_argument("--model", required=True) + ap.add_argument("--gtcs", required=True, type=Path, help="normal base samples (not in the model's reference)") + ap.add_argument("--out", required=True, type=Path) + ap.add_argument("--threads", type=int, default=0) + ap.add_argument("--seed", type=int, default=7) + ap.add_argument("--reps", type=int, default=4, help="replicates of every (kind, fraction, size) cell") + a = ap.parse_args() + rng = np.random.default_rng(a.seed) + man = read_bpm(a.bpm) + chroms = np.array(man.chroms) + pos = man.map_infos.astype(np.int64) + base = sorted(discover_gtcs(a.gtcs).values()) + syn = a.out / "gtcs" + syn.mkdir(parents=True, exist_ok=True) + events = [(k, f, s) for k, f in KINDS for s in SIZES for _ in range(a.reps)] + rng.shuffle(events) + truth = [] + per_sample = len(CHROMS) // 2 # one event per chromosome pair slot, spaced out + for n, start in enumerate(range(0, len(events), per_sample)): + chunk = events[start : start + per_sample] + src = base[n % len(base)] + g = copy.deepcopy(read_gtc(src)) + lrr = g.logr_ratios.astype(np.float64) + baf = g.b_allele_freqs.astype(np.float64) + chosen = rng.permutation(CHROMS)[: len(chunk)] + sid = f"syn{n:03d}" + for (kind, f, size), c in zip(chunk, chosen, strict=True): + sel = np.flatnonzero((chroms == c) & (pos > 0)) + lo, hi = int(pos[sel].min()), int(pos[sel].max()) + s0 = int(rng.integers(lo + 2_000_000, max(lo + 2_000_001, hi - size - 2_000_000))) + idx = sel[(pos[sel] >= s0) & (pos[sel] < s0 + size)] + implant(lrr, baf, idx, kind, f, rng) + truth.append({"sample": sid, "chrom": c, "start": s0, "end": s0 + size - 1, "kind": kind, "fraction": f, + "size": size, "probes": len(idx), "base": src.stem}) # fmt: skip + g.logr_ratios = lrr.astype(np.float32) + g.b_allele_freqs = baf.astype(np.float32) + (syn / sid).mkdir(exist_ok=True) + write_gtc(g, syn / sid / f"{sid}.gtc") + (a.out / "truth.json").write_text(json.dumps(truth, indent=1)) + simurg = str(Path(sys.executable).with_name("simurg")) + cmd = [simurg, "cn", "call", "--bpm", a.bpm, "--model", a.model, "--gtcs", str(syn), "--out", str(a.out / "call"), + "--threads", str(a.threads), "--unknown-sex-policy", "female", "--force", "--quiet"] # fmt: skip + subprocess.run(cmd, check=True) + calls: dict[str, list[tuple[str, int, int, str, str]]] = {} + for vcf in (a.out / "call").glob("*.cyto.cnv.vcf.gz"): + sid = vcf.name.split(".")[0] + with gzip.open(vcf, "rt") as fh: + lines = [ln for ln in fh if not ln.startswith("#")] + for line in lines: + c, p, vid, _, _alt, _q, filt, info = line.split("\t")[:8] + end = int(next(x for x in info.split(";") if x.startswith("END="))[4:]) + call = vid.split(":")[0] + cls = {"DEL": "loss", "MosaicDEL": "loss", "DUP": "gain", "MosaicDUP": "gain", "AOH": "loh", + "CNLOH": "loh", "GAINLOH": "gain"}[call] # fmt: skip + calls.setdefault(sid, []).append((c, int(p) + 1, end, cls, filt)) + for t in truth: + cover = 0 + seen = [] + for c, s, e, cls, filt in calls.get(t["sample"], []): + ov = max(0, min(e, t["end"]) - max(s, t["start"]) + 1) if c == t["chrom"] else 0 + if ov: + seen.append(f"{cls}:{filt}:{ov / (t['end'] - t['start'] + 1):.2f}") + if cls == CLASS[t["kind"]]: + cover += ov # fragments of one event add up + frac = cover / (t["end"] - t["start"] + 1) + t["detected"] = frac >= 0.5 + t["overlap"] = round(frac, 3) + t["overlapping_calls"] = seen + cells: dict[str, dict[str, list[bool]]] = {} + for t in truth: + key = t["kind"] if t["fraction"] == 1.0 else f"{t['kind']} f={t['fraction']:.2f}" + cells.setdefault(key, {}).setdefault(f"{t['size'] // 1_000_000} Mb", []).append(t["detected"]) + table = {k: {s: f"{sum(v)}/{len(v)}" for s, v in d.items()} for k, d in cells.items()} + n_calls = sum(len(v) for v in calls.values()) + implanted_hits = sum(t["detected"] for t in truth) + summary = {"events": len(truth), "detected": implanted_hits, "table": table, "calls_total": n_calls} + (a.out / "summary.json").write_text(json.dumps(summary, indent=1)) + (a.out / "truth_detected.json").write_text(json.dumps(truth, indent=1)) + for k, v in table.items(): + print(f"{k:14s} " + " ".join(f"{s}: {v.get(s, '-')}" for s in ("1 Mb", "5 Mb", "20 Mb"))) + print(f"detected {implanted_hits}/{len(truth)}; {n_calls} calls in {len(calls)} samples") + for t in truth: + if not t["detected"] and t["kind"] in ("DUP", "DEL", "LOH") and t["size"] >= 5_000_000: + print(" missed", t["kind"], t["size"], t["probes"], "probes:", t["overlapping_calls"]) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/src/simurg/pipeline/cn.py b/src/simurg/pipeline/cn.py index 2ea9a0f..8d8214f 100644 --- a/src/simurg/pipeline/cn.py +++ b/src/simurg/pipeline/cn.py @@ -88,9 +88,10 @@ def run_build_model( n_ref = 0 if reference_gtcs is not None: from simurg.cn.model import reference_stats + from simurg.formats.discovery import discover_gtcs from simurg.formats.gtc import read_gtc - paths = sorted(reference_gtcs.rglob("*.gtc")) + paths = sorted(discover_gtcs(reference_gtcs).values()) # E_NO_INPUT when there is none if len(paths) > max_reference: idx = np.random.default_rng(42).choice(len(paths), max_reference, replace=False) paths = [paths[i] for i in sorted(idx)]