"""Genome-wide CN/LOH/mosaic calling (spec CYTO-001 … CYTO-040).""" from __future__ import annotations import gzip import json from pathlib import Path import numpy as np import pytest 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, 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 from simurg.pipeline.cn import _mask_events, moving_median from simurg.synth.files import SynthLocus, encode_bpm # ------------------------------------------------------------------------------------------------ covariates def test_gc_fractions() -> None: # CYTO-038 seq = np.frombuffer(b"GGGGAAAAAANNNNCCCC", dtype=np.uint8) out = gc_fractions(seq, np.array([1, 13, 18]), (4, 36)) assert out[0, 0] == pytest.approx(1.0) # window [0, 2) clipped: "GG" assert np.isnan(out[1, 0]) # window [10, 14) = "NNNN": no called base assert out[2, 1] == pytest.approx(8 / 14) # whole sequence: 8 G/C over 14 called bases def test_population_bfreq() -> None: p = population_bfreq(np.array([[10, 0, 0], [0, 0, 10], [5, 10, 5], [0, 0, 0]], dtype=np.int32), 4) assert p.tolist() == pytest.approx([0.0, 1.0, 0.5, 0.5]) def test_gc_correct_removes_wave() -> None: # CYTO-003 rng = np.random.default_rng(1) gc = rng.uniform(0.3, 0.6, (5000, 3)).astype(np.float32) noise = rng.normal(0, 0.05, 5000) lrr = 0.8 * (gc[:, 2] - 0.45) + noise out, fit = gc_correct(lrr, gc, np.ones(5000, dtype=np.bool_)) assert fit.applied and fit.gcwf > 0.02 assert abs(np.corrcoef(out, gc[:, 2])[0, 1]) < 0.05 same, fit2 = gc_correct(lrr[:10], gc[:10], np.ones(10, dtype=np.bool_)) assert not fit2.applied and np.array_equal(same, lrr[:10]) def test_moving_median_truncates_ends() -> None: # CYTO-004 x = np.array([1.0, 9.0, 2.0, 3.0, np.nan, 4.0]) m = moving_median(x, 1) assert m[0] == pytest.approx(5.0) # window truncated at the start: median(1, 9) assert m[2] == pytest.approx(3.0) and np.isnan(m[4]) assert np.array_equal(moving_median(x, 0), x, equal_nan=True) # ------------------------------------------------------------------------------------------------ HMM + events def _simulate(n: int = 4000, seed: int = 3) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: rng = np.random.default_rng(seed) pos = np.arange(1, n + 1, dtype=np.int64) * 5000 pfb = np.full(n, 0.3) geno = rng.choice([0, 1, 2], n, p=[0.49, 0.42, 0.09]) baf = np.choose(geno, [0.0, 0.5, 1.0]) + rng.normal(0, 0.03, n) lrr = rng.normal(0, 0.15, n) lrr[1000:1400] -= 0.45 # CN1 deletion: heterozygotes disappear baf[1000:1400] = np.where(baf[1000:1400] > 0.25, 1.0, 0.0) + rng.normal(0, 0.02, 400) baf[2500:3300] = np.where(baf[2500:3300] > 0.25, 1.0, 0.0) + rng.normal(0, 0.02, 800) # copy-neutral LOH return pos, lrr, np.clip(baf, 0, 1), pfb def test_hmm_finds_deletion_and_loh() -> None: pos, lrr, baf, pfb = _simulate() st = hmm.diploid_states() prior = np.array([s.prior for s in st]) prior /= prior.sum() e = hmm.emissions(st, lrr, baf, pfb, hmm.Params(0.15, 0.03)) sw = hmm.switch_probs(pos) path = hmm.viterbi(e, sw, prior) evs = extract("1", path, st, pos) calls = {(ev.call, ev.i0 // 100, ev.i1 // 100) for ev in evs if ev.probes >= 100} assert ("DEL", 10, 14) in calls assert ("LOH", 25, 33) in calls post = hmm.forward_backward(e, sw, prior) assert np.allclose(post.sum(axis=1), 1.0) assert post[1200, [0, 1]].sum() > 0.9 def test_effective_size_and_thresholds() -> None: # CYTO-007/008/010 pos = np.array([100, 200, 300, 400, 100_400, 100_500], dtype=np.int64) assert effective_size(pos, 100.0) == (100_500 - 100 + 1) - (100_000 - 100) th = Thresholds.defaults() ev = Event("1", 1, 30_000, 0, 12, "DEL", 1, "CN1", effective_size=30_000) small = Event("1", 1, 20_000, 0, 12, "DEL", 1, "CN1", effective_size=20_000) counts: dict[str, int] = {} assert apply_thresholds([ev, small], th, counts) == [ev] and counts["below_threshold"] == 1 def test_mask_trimming() -> None: # CYTO-017 mask = {"1": [(1_000, 1_999)]} a = Event("1", 1_500, 10_000, 0, 20, "DEL", 1, "CN1") b = Event("1", 1_100, 2_500, 0, 20, "DEL", 1, "CN1") kept = _mask_events([a, b], mask, counts := {}) assert kept == [a] and a.start == 2_000 and "centromere_trimmed" in a.flags and counts["masked"] == 1 # ------------------------------------------------------------------------------------------------ CLI end to end def _dataset(tmp: Path) -> tuple[Path, Path, Path]: rng = np.random.default_rng(11) n_per = 3000 loci = [] for c in ("1", "2"): for k in range(n_per): loci.append(SynthLocus(f"rs{c}_{k}", c, 5000 * (k + 1), "[A/G]", len(loci) + 1, 0, 0)) bpm = tmp / "m.bpm" bpm.write_bytes(encode_bpm(loci)) fa = tmp / "g.fa" body = {c: "".join(rng.choice(list("ACGT"), n_per * 5000 + 10_000)) for c in ("1", "2")} fa.write_text("".join(f">{c}\n" + "\n".join(s[i : i + 60] for i in range(0, len(s), 60)) + "\n" for c, s in body.items())) # fmt: skip gtcs = tmp / "gtcs" n = len(loci) for sid, deletion in (("200000000001_R01C01", True), ("200000000001_R02C01", False)): geno = rng.choice([1, 2, 3], n, p=[0.49, 0.42, 0.09]).astype(np.uint8) baf = np.choose(geno - 1, [0.0, 0.5, 1.0]) + rng.normal(0, 0.03, n) lrr = rng.normal(0, 0.12, n) if deletion: lrr[500:900] -= 0.45 baf[500:900] = np.where(baf[500:900] > 0.25, 1.0, 0.0) geno[500:900] = np.where(geno[500:900] == 2, 3, geno[500:900]) g = Gtc(num_snps=n, genotypes=geno, genotype_scores=np.full(n, 0.9, np.float32), b_allele_freqs=np.clip(baf, 0, 1).astype(np.float32), logr_ratios=lrr.astype(np.float32), raw_x=np.full(n, 1000, np.uint16), raw_y=np.full(n, 1000, np.uint16), base_calls=np.zeros((n, 2), np.uint8), gender=b"F") # fmt: skip (gtcs / sid).mkdir(parents=True) write_gtc(g, gtcs / sid / f"{sid}.gtc") return bpm, fa, gtcs def test_cn_build_and_call(tmp_path: Path) -> None: # CYTO-028 .. CYTO-032, CYTO-037 bpm, fa, gtcs = _dataset(tmp_path) model = tmp_path / "m.simurg-cyto" r = CliRunner().invoke(app, ["cn", "build-model", "--bpm", str(bpm), "--fasta", str(fa), "--out", str(model), "--quiet"]) # fmt: skip assert r.exit_code == 0, r.output out = tmp_path / "cn" r = CliRunner().invoke(app, ["cn", "call", "--bpm", str(bpm), "--model", str(model), "--gtcs", str(gtcs), "--fasta", str(fa), "--out", str(out), "--threads", "1", "--min-del-probes", "5", "--unknown-sex-policy", "female", "--quiet"]) # fmt: skip assert r.exit_code == 0, r.output side = json.loads((out / "cn_call_manifest.json").read_text()) assert side["thresholds"]["min_del_probes"] == 10 # CYTO-008 floor with gzip.open(out / "200000000001_R01C01.cyto.cnv.vcf.gz", "rt") as fh: text = fh.read() assert (out / "200000000001_R01C01.cyto.cnv.vcf.gz.tbi").exists() head = [ln for ln in text.splitlines() if ln.startswith("##")] assert head[0] == "##fileformat=VCFv4.3" and any(h.startswith("##simurg_thresholds=") for h in head) recs = [ln.split("\t") for ln in text.splitlines() if not ln.startswith("#")] dels = [r_ for r_ in recs if r_[4] == ""] assert dels, text[-2000:] for r_ in recs: # CYTO-029 / EDGE-025 pos, rid, info = int(r_[1]), r_[2], r_[7] end = int(next(x for x in info.split(";") if x.startswith("END="))[4:]) s, e = rid.split(":")[2].split("-") assert int(s) == pos + 1 and int(e) == end assert r_[3] in "ACGT" svlen = int(next(x for x in info.split(";") if x.startswith("SVLEN="))[6:]) assert (svlen < 0) == (r_[4] == "") assert r_[8] == "GT:CN:NR:LRD:MF" d = dels[0] assert d[0] == "1" and 2_300_000 < int(d[1]) < 2_700_000 and d[9].split(":")[:2] == ["0/1", "1"] summ = (out / "cn_sample_summary.tsv").read_text().splitlines() assert summ[0].startswith("sample_id\tsample_name\tstatus") and len(summ) == 3 assert (out / "200000000001_R01C01.cn_segments.parquet").exists() # ------------------------------------------------------------------------------------------------ annotate def test_iscn_strings_and_genes() -> None: # CYTO-034 / CYTO-035 from simurg.cn.annotate import Bands, Genes, iscn_sample bands = Bands.from_ucsc("chr1\t0\t2300000\tp36.33\tgneg\nchr1\t2300000\t5300000\tp36.32\tgpos25\n" "chr9\t0\t2200000\tp24.3\tgneg\n") # fmt: skip dele = {"chrom": "1", "start": 1_567_924, "end": 1_632_825, "call": "DEL", "call_type": "DEL", "copy_number": 1, "mosaic": False} # fmt: skip dup = {"chrom": "1", "start": 2_000_000, "end": 3_000_000, "call": "DUP", "call_type": "DUP", "copy_number": 3, "mosaic": False} # fmt: skip loh = {"chrom": "9", "start": 100, "end": 2_000_000, "call": "LOH", "call_type": "LOH", "copy_number": 2, "mosaic": False} # fmt: skip mos = {"chrom": "9", "start": 100, "end": 1_000, "call": "MDEL", "call_type": "DEL", "copy_number": 1, "mosaic": True, "base_cn": 2} # fmt: skip trad, micro = iscn_sample([loh, dup, dele], bands, "GRCh38", "F") assert micro == ("arr[GRCh38] 1p36.33(1567924_1632825)x1,1p36.33p36.32(2000000_3000000)x3," "9p24.3(100_2000000)x2 hmz") # fmt: skip assert trad == "del(1)(p36.33p36.33),dup(1)(p36.33p36.32)" assert iscn_sample([mos], bands, "GRCh38", "F")[1] == "arr[GRCh38] 9p24.3(100_1000)x2~1" assert iscn_sample([], bands, "GRCh38", "M")[1] == "arr[GRCh38](1-22)x2,(XY)x1" assert iscn_sample([], bands, "GRCh37", "F")[1] == "arr[GRCh37](1-22,X)x2" assert iscn_sample([], bands, "GRCh38", "U")[1] == "arr[GRCh38](1-22)x2" genes = Genes.from_rows([("1", 100, 5_000, "A"), ("1", 4_000, 4_500, "B"), ("1", 200, 300, "A"), ("1", 9_000, 9_100, "C")]) # fmt: skip assert genes.overlapping("1", 4_400, 8_999) == ["A", "B"] def test_cn_annotate_end_to_end(tmp_path: Path) -> None: # CYTO-032 / CYTO-033 bpm, fa, gtcs = _dataset(tmp_path) model, out = tmp_path / "m.simurg-cyto", tmp_path / "cn" run = CliRunner().invoke assert run(app, ["cn", "annotate", "--vcfs", str(tmp_path)]).exit_code == 2 # E_CN_SIDECAR_MISSING assert run(app, ["cn", "build-model", "--bpm", str(bpm), "--fasta", str(fa), "--out", str(model), "--quiet"]).exit_code == 0 # fmt: skip r = run(app, ["cn", "call", "--bpm", str(bpm), "--model", str(model), "--gtcs", str(gtcs), "--fasta", str(fa), "--out", str(out), "--threads", "1", "--unknown-sex-policy", "female", "--quiet"]) # fmt: skip assert r.exit_code == 0, r.output r = run(app, ["cn", "annotate", "--vcfs", str(out), "--min-qual", "0", "--quiet"]) assert r.exit_code == 0, r.output doc = json.loads((out / "200000000001_R01C01.cyto.json").read_text()) assert doc["schema"] == "simurg.cyto/1" and doc["annotation_type"] == "constitutional" assert doc["iscn_microarray"].startswith("arr[") assert any(e["call_type"] == "DEL" and e["chrom"] == "1" for e in doc["events"]) assert {c["id"] for c in doc["chromosomes"]} == {"1", "2"} normal = json.loads((out / "200000000001_R02C01.cyto.json").read_text()) assert normal["iscn_microarray"] in ("arr[GRCh38](1-22,X)x2", "arr[GRCh37](1-22,X)x2") or normal["events"] def test_packaged_cn_data() -> None: # the built-in mask and cytobands ship with the package (CYTO-017, CYTO-035) from simurg.pipeline.cn import _builtin_file, read_mask for build in ("GRCh37", "GRCh38"): 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) 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 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)) def test_probe_lrr_noise() -> None: # CYTO-042 from simurg.cn.model import reference_lrr_spread from simurg.pipeline.cn import _lrr_noise rng = np.random.default_rng(5) sd = np.array([0.1] * 50 + [0.3] * 50) lrr = (rng.normal(0, 1, (40, 100)) * sd + rng.normal(0, 0.2, (40, 1))).astype(np.float32) # + sample offsets spread = reference_lrr_spread(lrr, np.zeros(100, dtype=np.float32)) assert abs(float(np.median(spread[:50])) - 0.1) < 0.02 and abs(float(np.median(spread[50:])) - 0.3) < 0.05 auto = np.ones(100, dtype=np.bool_) s = _lrr_noise(0.2, auto, spread) assert np.median(s[50:]) / np.median(s[:50]) > 2 assert np.all(_lrr_noise(0.2, auto, None) == 0.2)