kemalandClaude Opus 5.5
662addf5eb
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). Co-Authored-By: Claude Opus 5.5 <[email protected]> Claude-Session: https://claude.ai/code/session_01814GKVzZVL3FuF3DeiZXhv
252 lines
13 KiB
Python
252 lines
13 KiB
Python
"""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] == "<DEL>"]
|
|
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] == "<DEL>")
|
|
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))
|