kemal
4315f37f1f
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).
258 lines
10 KiB
Python
258 lines
10 KiB
Python
"""CN/LOH event extraction, statistics and filters — spec CYTO-006 … CYTO-021, CYTO-027. Pure functions, no I/O.
|
||
|
||
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/042);
|
||
post-filters in the order of CYTO-015; per-type thresholds (CYTO-007/008) and promotion/low fraction (CYTO-020/021).
|
||
"""
|
||
|
||
from __future__ import annotations
|
||
|
||
import math
|
||
from dataclasses import dataclass, field
|
||
|
||
import numpy as np
|
||
import numpy.typing as npt
|
||
|
||
from simurg.cn.hmm import State
|
||
|
||
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
|
||
CALL_KEY = {"DEL": "del", "DUP": "dup", "LOH": "loh", "GAINLOH": "gainloh", "MDEL": "mdel", "MDUP": "mdup",
|
||
"MLOH": "mloh"} # fmt: skip
|
||
|
||
|
||
@dataclass
|
||
class Event:
|
||
chrom: str
|
||
start: int # first affected base S (1-based) = first probe position
|
||
end: int # last affected base E = last probe position
|
||
i0: int # probe index range in the chromosome arrays [i0, i1)
|
||
i1: int
|
||
call: str # DEL | DUP | LOH | GAINLOH | MDEL | MDUP | MLOH
|
||
cn: int
|
||
state: str
|
||
quality: float = 0.0
|
||
effective_size: int = 0
|
||
nr: float = math.nan
|
||
lrd: float = math.nan
|
||
baf_dev: float = math.nan
|
||
percent_het: float = math.nan
|
||
lrr_median: float = math.nan
|
||
mosaic_fraction: float | None = None
|
||
high_fraction: bool = False
|
||
filters: list[str] = field(default_factory=list)
|
||
flags: list[str] = field(default_factory=list)
|
||
|
||
@property
|
||
def mosaic(self) -> bool:
|
||
return self.call.startswith("M")
|
||
|
||
@property
|
||
def probes(self) -> int:
|
||
return self.i1 - self.i0
|
||
|
||
@property
|
||
def size(self) -> int:
|
||
return self.end - self.start + 1
|
||
|
||
@property
|
||
def svtype(self) -> str:
|
||
return {"MDEL": "DEL", "MDUP": "DUP", "MLOH": "LOH", "GAINLOH": "LOH"}.get(self.call, self.call)
|
||
|
||
@property
|
||
def loh_type(self) -> str | None:
|
||
if self.call == "LOH":
|
||
return "AOH"
|
||
if self.call == "GAINLOH":
|
||
return "AOH"
|
||
if self.call == "MLOH":
|
||
return "CNLOH"
|
||
return None
|
||
|
||
@property
|
||
def id_prefix(self) -> str:
|
||
return {"DEL": "DEL", "DUP": "DUP", "LOH": "AOH", "GAINLOH": "AOH", "MDEL": "MosaicDEL", "MDUP": "MosaicDUP",
|
||
"MLOH": "CNLOH"}[self.call] # fmt: skip
|
||
|
||
|
||
def median_spacing(pos: I64) -> float:
|
||
d = np.diff(np.sort(pos))
|
||
d = d[d > 0]
|
||
return float(np.median(d)) if len(d) else 1.0
|
||
|
||
|
||
def effective_size(pos: I64, m: float) -> int:
|
||
"""CYTO-010 for the probes of one event (sorted positions)."""
|
||
if len(pos) < 2:
|
||
return 1
|
||
s = np.diff(pos).astype(np.float64)
|
||
gap = s > GAP_FACTOR * m
|
||
return int((pos[-1] - pos[0] + 1) - float((s[gap] - m).sum()))
|
||
|
||
|
||
def extract(chrom: str, path: npt.NDArray[np.int32], states: list[State], pos: I64) -> list[Event]:
|
||
"""Runs of non-baseline states, merged across runs of the same call type (mosaic fraction steps)."""
|
||
out: list[Event] = []
|
||
if len(path) == 0:
|
||
return out
|
||
br = np.flatnonzero(np.diff(path) != 0) + 1
|
||
starts, ends = np.r_[0, br], np.r_[br, len(path)]
|
||
for s, e in zip(starts.tolist(), ends.tolist(), strict=True):
|
||
st = states[int(path[s])]
|
||
if not st.call:
|
||
continue
|
||
if out and out[-1].i1 == s and out[-1].call == st.call and out[-1].cn == st.cn:
|
||
out[-1].i1 = e
|
||
out[-1].end = int(pos[e - 1])
|
||
continue
|
||
out.append(Event(chrom, int(pos[s]), int(pos[e - 1]), s, e, st.call, st.cn, st.name))
|
||
return out
|
||
|
||
|
||
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_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":
|
||
f = 2 * d / (0.5 - d) if d < 0.5 else 1.0
|
||
else:
|
||
f = 2 * d
|
||
return float(min(1.0, max(0.0, f)))
|
||
|
||
|
||
def _mf_from_lrr(call: str, lrr: F64, gamma: float) -> float | None:
|
||
x = lrr[np.isfinite(lrr)]
|
||
if len(x) < MIN_MF_PROBES or gamma <= 0:
|
||
return None
|
||
f = abs(2 ** (float(np.median(x)) / gamma) - 1)
|
||
return float(min(1.0, max(0.0, f)))
|
||
|
||
|
||
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, 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)
|
||
x = lrr[sl]
|
||
fx = x[np.isfinite(x)]
|
||
ev.nr = float(np.median(2.0**fx)) if len(fx) else math.nan
|
||
ev.lrd = float(np.std(fx)) if len(fx) > 1 else math.nan
|
||
ev.lrr_median = float(np.median(fx)) if len(fx) else math.nan
|
||
g = geno[sl]
|
||
called = g > 0
|
||
ev.percent_het = float(100.0 * (g == 2).sum() / called.sum()) if called.any() else math.nan
|
||
b = baf[sl][(g == 2) & np.isfinite(baf[sl])]
|
||
ev.baf_dev = float(np.std(b)) if len(b) > 1 else math.nan
|
||
ev.effective_size = effective_size(pos[sl], m)
|
||
cls = [j for j, st in enumerate(states) if st.call == ev.call]
|
||
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], 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)
|
||
|
||
|
||
@dataclass
|
||
class Thresholds:
|
||
probes: dict[str, int]
|
||
size_bp: dict[str, int]
|
||
max_mosaic_fraction: float = 1.0
|
||
min_mosaic_fraction: float = 0.10
|
||
filter_qual: int = 10
|
||
lrr_direction_min_frac: float = 0.6
|
||
|
||
@classmethod
|
||
def defaults(cls) -> Thresholds:
|
||
return cls({k: v[0] for k, v in DEFAULTS.items()}, {k: v[1] for k, v in DEFAULTS.items()})
|
||
|
||
def as_header(self, extra: dict[str, object]) -> str:
|
||
parts = []
|
||
for k in ("del", "dup", "loh", "gainloh", "mdel", "mdup", "mloh"):
|
||
parts += [f"min_{k}_probes={self.probes[k]}", f"min_{k}_size_bp={self.size_bp[k]}"]
|
||
parts += [f"{k}={v}" for k, v in extra.items()]
|
||
return ";".join(parts)
|
||
|
||
|
||
def post_filter(events: list[Event], lrr_by_event: list[F64], th: Thresholds, counts: dict[str, int]) -> list[Event]:
|
||
"""CYTO-015 steps 1–3 and CYTO-020 promotion; per-type thresholds are applied by :func:`apply_thresholds`."""
|
||
kept: list[Event] = []
|
||
for ev, x in zip(events, lrr_by_event, strict=True):
|
||
fx = x[np.isfinite(x)]
|
||
if ev.call in ("DEL", "MDEL") and len(fx) and (fx < 0).mean() < th.lrr_direction_min_frac:
|
||
counts["lrr_direction"] = counts.get("lrr_direction", 0) + 1
|
||
continue
|
||
if ev.call in ("DUP", "MDUP", "GAINLOH") and len(fx) and (fx > 0).mean() < th.lrr_direction_min_frac:
|
||
counts["lrr_direction"] = counts.get("lrr_direction", 0) + 1
|
||
continue
|
||
kept.append(ev)
|
||
# short LOH between copy-neutral neighbours (no other event adjacent on either side)
|
||
out: list[Event] = []
|
||
for i, ev in enumerate(kept):
|
||
if ev.call in ("LOH", "MLOH") and ev.size < 500_000:
|
||
left = kept[i - 1] if i > 0 and kept[i - 1].chrom == ev.chrom else None
|
||
right = kept[i + 1] if i + 1 < len(kept) and kept[i + 1].chrom == ev.chrom else None
|
||
if not (left and left.i1 == ev.i0) and not (right and right.i0 == ev.i1):
|
||
counts["short_loh"] = counts.get("short_loh", 0) + 1
|
||
continue
|
||
out.append(ev)
|
||
merged: list[Event] = []
|
||
for ev in out:
|
||
prev = merged[-1] if merged else None
|
||
if (prev and prev.chrom == ev.chrom and prev.call == ev.call and prev.cn == ev.cn
|
||
and ev.i0 - prev.i1 <= 3 and (ev.start - prev.end) <= 0.1 * (ev.end - prev.start + 1)): # fmt: skip
|
||
prev.i1, prev.end = ev.i1, ev.end
|
||
counts["merged"] = counts.get("merged", 0) + 1
|
||
continue
|
||
merged.append(ev)
|
||
for ev in merged: # CYTO-020 promotion
|
||
if ev.mosaic and ev.mosaic_fraction is not None and ev.mosaic_fraction >= th.max_mosaic_fraction:
|
||
ev.call = {"MDEL": "DEL", "MDUP": "DUP", "MLOH": "LOH"}[ev.call]
|
||
ev.high_fraction = True
|
||
return merged
|
||
|
||
|
||
def apply_thresholds(events: list[Event], th: Thresholds, counts: dict[str, int]) -> list[Event]:
|
||
out = []
|
||
for ev in events:
|
||
k = CALL_KEY[ev.call]
|
||
if ev.probes < th.probes[k] or ev.effective_size < th.size_bp[k]:
|
||
counts["below_threshold"] = counts.get("below_threshold", 0) + 1
|
||
continue
|
||
out.append(ev)
|
||
return out
|
||
|
||
|
||
def set_filters(events: list[Event], th: Thresholds, sample_noisy: bool) -> None:
|
||
"""CYTO-013 FILTER values."""
|
||
for ev in events:
|
||
f = []
|
||
if ev.quality < th.filter_qual:
|
||
f.append(f"Q{th.filter_qual}")
|
||
if sample_noisy:
|
||
f.append("SampleQuality")
|
||
if ev.mosaic and ev.mosaic_fraction is not None and ev.mosaic_fraction < th.min_mosaic_fraction:
|
||
f.append("LowFraction")
|
||
ev.filters = f or ["PASS"]
|