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).
This commit is contained in:
1 parent
4315f37f1f
commit
aa3031a18a
4 files changed
+90
-31
No files matched your search
@@ -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_])]
|
||||
|
||||
Reference in new issue
Block a user