CN: noise-corrected BAF split and complete-loss MF (CYTO-041); Pi tuning harness
Mosaic fraction from the BAF split now subtracts the heterozygous BAF noise (median |BAF-0.5| of unsplit hets = 0.6745 sigma_b), and a mosaic loss lacking its expected heterozygotes gets MF 1 and is promoted to a constitutional DEL (CYTO-020). Adds scripts/v2/cn_tune.sh + cn_tune_report.py (GSE165845 rounds on the Pi).
This commit is contained in:
1 parent
98e16cf501
commit
6bc42e53b7
6 files changed
+140
-8
No files matched your search
+18
-6
@@ -2,7 +2,7 @@
|
||||
|
||||
From a per-chromosome Viterbi path: runs of non-baseline states become events; statistics per CYTO-027; effective
|
||||
size per CYTO-010 (gaps > 150 × the chromosome's median spacing replaced by the median); QUAL from the mean posterior
|
||||
of the called state class (CYTO-012); mosaic fraction from the BAF split or, failing that, LRR (CYTO-018/019);
|
||||
of the called state class (CYTO-012); mosaic fraction from the BAF split or, failing that, LRR (CYTO-018/019/042);
|
||||
post-filters in the order of CYTO-015; per-type thresholds (CYTO-007/008) and promotion/low fraction (CYTO-020/021).
|
||||
"""
|
||||
|
||||
@@ -20,6 +20,8 @@ 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
|
||||
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
|
||||
@@ -115,11 +117,21 @@ def extract(chrom: str, path: npt.NDArray[np.int32], states: list[State], pos: I
|
||||
return out
|
||||
|
||||
|
||||
def _mf_from_baf(call: str, baf: F64) -> float | None:
|
||||
b = baf[np.isfinite(baf) & (baf > 0.1) & (baf < 0.9)]
|
||||
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
|
||||
if len(b) < MIN_MF_PROBES:
|
||||
return None
|
||||
d = float(np.median(np.abs(b - 0.5)))
|
||||
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 call == "MDEL":
|
||||
f = 2 * d / (0.5 + d)
|
||||
elif call == "MDUP":
|
||||
@@ -139,7 +151,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,
|
||||
gamma: float, haploid: bool, baf_sd: float = 0.0, pfb: F64 | None = None,
|
||||
) -> None: # fmt: skip
|
||||
"""CYTO-010/012/018/027 statistics in place (``lrr`` unsmoothed)."""
|
||||
sl = slice(ev.i0, ev.i1)
|
||||
@@ -158,7 +170,7 @@ 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])
|
||||
mf = None if haploid else _mf_from_baf(ev.call, baf[sl], baf_sd, None if pfb is None else pfb[sl])
|
||||
ev.mosaic_fraction = mf if mf is not None else _mf_from_lrr(ev.call, x, gamma)
|
||||
|
||||
|
||||
|
||||
@@ -392,7 +392,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)
|
||||
annotate_stats(ev, p, x_raw, b, gt, post, states, m_sp, ctx["gamma"], haploid, baf_sd, 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