v2 CN/LOH (WIP): cn build-model and cn call
Cyto model (.simurg-cyto: GC 1 kb/100 kb/1 Mb, PFB from the cluster file, optional sex-aware batch reference from normal GTCs; IN-051, CYTO-038); GC-wave correction (CYTO-003); smoothing (CYTO-004); HMM with constitutional, gain-LOH, mosaic and haploid states, forward-backward QUAL and Viterbi segments (CYTO-005/012); event statistics, effective size, post-filters, per-type thresholds, promotion/LowFraction, built-in centromere mask from UCSC cytoBand (CYTO-006..021); sex resolution (CYTO-024), ploidy; CNV VCF with bgzip+tabix (CYTO-028..031), segments parquet, cn_sample_summary.tsv and the cn_call_manifest.json sidecar.
This commit is contained in:
kemal committed
1 parent
5164124438
commit
edb59e4ee2
8 files changed
+1452
No files matched your search
@@ -0,0 +1,245 @@
|
||||
"""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);
|
||||
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
|
||||
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) -> float | None:
|
||||
b = baf[np.isfinite(baf) & (baf > 0.1) & (baf < 0.9)]
|
||||
if len(b) < MIN_MF_PROBES:
|
||||
return None
|
||||
d = float(np.median(np.abs(b - 0.5)))
|
||||
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,
|
||||
) -> 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])
|
||||
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"]
|
||||
@@ -0,0 +1,38 @@
|
||||
"""GC-wave correction of LRR — spec CYTO-003, plan §10.1. Pure functions, no I/O.
|
||||
|
||||
[LIT: Diskin et al. 2008 NAR 36:e126] LRR is regressed on the probe-centred GC fractions (1 kb, 100 kb, 1 Mb) by
|
||||
ordinary least squares over autosomal probes with |LRR| < 1 (CNV and outlier probes do not steer the fit); the GC
|
||||
component (fitted value minus its mean) is subtracted. [OURS] GCWF is reported as the standard deviation of that
|
||||
removed component, i.e. the size of the wave in LRR units.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from dataclasses import dataclass
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
F64 = npt.NDArray[np.float64]
|
||||
MIN_PROBES = 1000
|
||||
|
||||
|
||||
@dataclass
|
||||
class GcFit:
|
||||
applied: bool
|
||||
coef: F64 # intercept + one slope per window
|
||||
gcwf: float
|
||||
|
||||
|
||||
def gc_correct(lrr: F64, gc: npt.NDArray[np.float32], use: npt.NDArray[np.bool_]) -> tuple[F64, GcFit]:
|
||||
"""Return corrected LRR (NaN stays NaN) and the fit; unchanged LRR when too few probes qualify."""
|
||||
x = np.asarray(gc, dtype=np.float64)
|
||||
ok = use & np.isfinite(lrr) & (np.abs(lrr) < 1.0) & np.isfinite(x).all(axis=1)
|
||||
if ok.sum() < MIN_PROBES:
|
||||
return lrr.copy(), GcFit(False, np.zeros(x.shape[1] + 1), 0.0)
|
||||
xm = x[ok].mean(axis=0)
|
||||
design = np.column_stack([np.ones(int(ok.sum())), x[ok] - xm])
|
||||
coef, *_ = np.linalg.lstsq(design, lrr[ok], rcond=None)
|
||||
fitted = np.where(np.isfinite(x).all(axis=1), (np.nan_to_num(x) - xm) @ coef[1:], 0.0)
|
||||
out = lrr - fitted
|
||||
return out, GcFit(True, coef, float(np.std(fitted[ok])))
|
||||
@@ -0,0 +1,190 @@
|
||||
"""Genome-wide CN/LOH/mosaic hidden Markov model — spec CYTO-004/005/012/018, plan §10.2. Pure functions, no I/O.
|
||||
|
||||
[LIT: Wang et al. 2007 Genome Res 17:1665 (PennCNV: LRR/BAF emissions weighted by population B-allele frequency,
|
||||
distance-dependent transitions); Loh et al. 2018 Nature 559:350 (MoChA: mosaic BAF split); DOC: cnvPartition
|
||||
technical note (reference LRR state means CN0 −5, CN1 −0.45, CN2 0, CN3 0.3, CN4 0.75)] [OURS: state grid, mixture
|
||||
weights, switch model]
|
||||
|
||||
* Diploid regions: CN0, CN1, CN2, CN2-LOH, CN3, CN3-LOH, CN4 and mosaic loss / gain / CN-LOH at cell fractions
|
||||
0.2, 0.4, 0.6, 0.8. Haploid regions (male X outside PAR, male Y): CN0, CN1, CN2 and mosaic gain/loss relative to
|
||||
the region's median LRR.
|
||||
* LRR ~ (1−ε)·N(µ_state, σ) + ε·U(−5, 3); σ = 1.4826·MAD of the sample's autosomal LRR (floor 0.08); CN0 σ = 2.
|
||||
Mosaic means interpolate linearly between the constitutional means by the cell fraction.
|
||||
* BAF ~ (1−ε)·Σ_g w_g N(b_g, σ_b) + ε·U(0, 1); genotype weights from the probe's PFB p (binomial in the copy
|
||||
number; LOH: {0: 1−p, 1: p}); mosaic heterozygotes split into the two bands 0.5 ± d(f).
|
||||
* Transitions: switch probability 1 − exp(−d / D) for inter-probe distance d, D = 5·10⁸ bp; a switch lands on state
|
||||
j with probability π_j (the same target distribution from every state).
|
||||
* Viterbi path for segmentation; forward–backward posteriors for QUAL.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import math
|
||||
from dataclasses import dataclass, field
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
F64 = npt.NDArray[np.float64]
|
||||
FRACTIONS = (0.2, 0.4, 0.6, 0.8)
|
||||
LRR_MEAN = {0: -5.0, 1: -0.45, 2: 0.0, 3: 0.3, 4: 0.75}
|
||||
EPS = 0.01
|
||||
MOSAIC_PRIOR = 0.004 # total switch-target mass of the 12 diploid mosaic states
|
||||
SWITCH_D = 5e8
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class State:
|
||||
name: str
|
||||
cn: int # copy number (baseline-relative for haploid regions: 0, 1, 2)
|
||||
loh: bool
|
||||
mosaic: str # "" | "loss" | "gain" | "cnloh"
|
||||
fraction: float # cell fraction of mosaic states, 1.0 otherwise
|
||||
call: str # "" (baseline) | DEL | DUP | LOH | GAINLOH | MDEL | MDUP | MLOH
|
||||
lrr_mu: float
|
||||
baf: tuple[tuple[str, float], ...] # (genotype-weight key, BAF mean)
|
||||
prior: float
|
||||
|
||||
|
||||
def _d_split(kind: str, f: float) -> float:
|
||||
if kind == "loss":
|
||||
return f / (2 * (2 - f))
|
||||
if kind == "gain":
|
||||
return f / (2 * (2 + f))
|
||||
return f / 2
|
||||
|
||||
|
||||
def diploid_states() -> list[State]:
|
||||
s = [
|
||||
State("CN0", 0, False, "", 1.0, "DEL", LRR_MEAN[0], (), 0.01),
|
||||
State("CN1", 1, False, "", 1.0, "DEL", LRR_MEAN[1], (("h0", 0.0), ("h1", 1.0)), 0.06),
|
||||
State("CN2", 2, False, "", 1.0, "", 0.0, (("aa", 0.0), ("ab", 0.5), ("bb", 1.0)), 0.60),
|
||||
State("CN2LOH", 2, True, "", 1.0, "LOH", 0.0, (("h0", 0.0), ("h1", 1.0)), 0.06),
|
||||
State("CN3", 3, False, "", 1.0, "DUP", LRR_MEAN[3],
|
||||
(("t0", 0.0), ("t1", 1 / 3), ("t2", 2 / 3), ("t3", 1.0)), 0.06),
|
||||
State("CN3LOH", 3, True, "", 1.0, "GAINLOH", LRR_MEAN[3], (("h0", 0.0), ("h1", 1.0)), 0.02),
|
||||
State("CN4", 4, False, "", 1.0, "DUP", LRR_MEAN[4],
|
||||
(("q0", 0.0), ("q1", 0.25), ("q2", 0.5), ("q3", 0.75), ("q4", 1.0)), 0.03),
|
||||
] # fmt: skip
|
||||
for kind, call, mu_end in (("loss", "MDEL", LRR_MEAN[1]), ("gain", "MDUP", LRR_MEAN[3]), ("cnloh", "MLOH", 0.0)):
|
||||
for f in FRACTIONS:
|
||||
d = _d_split(kind, f)
|
||||
cn = {"loss": 1, "gain": 3, "cnloh": 2}[kind]
|
||||
bands = (("aa", 0.0), ("ab-", 0.5 - d), ("ab+", 0.5 + d), ("bb", 1.0))
|
||||
s.append(State(f"M{kind}{f:.1f}", cn, kind == "cnloh", kind, f, call, mu_end * f, bands, MOSAIC_PRIOR / 12))
|
||||
return s
|
||||
|
||||
|
||||
def haploid_states() -> list[State]:
|
||||
s = [
|
||||
State("H0", 0, False, "", 1.0, "DEL", LRR_MEAN[0], (), 0.02),
|
||||
State("H1", 1, False, "", 1.0, "", 0.0, (("h0", 0.0), ("h1", 1.0)), 0.80),
|
||||
State("H2", 2, False, "", 1.0, "DUP", 0.45, (("h0", 0.0), ("h1", 1.0)), 0.06),
|
||||
]
|
||||
for f in FRACTIONS:
|
||||
s.append(State(f"Hgain{f:.1f}", 2, False, "gain", f, "MDUP", 0.45 * f, (("h0", 0.0), ("h1", 1.0)),
|
||||
MOSAIC_PRIOR / 8)) # fmt: skip
|
||||
s.append(State(f"Hloss{f:.1f}", 0, False, "loss", f, "MDEL", -1.0 * f, (("h0", 0.0), ("h1", 1.0)),
|
||||
MOSAIC_PRIOR / 8)) # fmt: skip
|
||||
return s
|
||||
|
||||
|
||||
def _geno_weights(key: str, p: F64) -> F64:
|
||||
q = 1 - p
|
||||
table = {
|
||||
"aa": q * q, "ab": 2 * p * q, "bb": p * p, "ab-": p * q, "ab+": p * q,
|
||||
"h0": q, "h1": p,
|
||||
"t0": q**3, "t1": 3 * p * q * q, "t2": 3 * p * p * q, "t3": p**3,
|
||||
"q0": q**4, "q1": 4 * p * q**3, "q2": 6 * p * p * q * q, "q3": 4 * p**3 * q, "q4": p**4,
|
||||
} # fmt: skip
|
||||
return table[key]
|
||||
|
||||
|
||||
def _norm(x: F64, mu: float, sd: float) -> F64:
|
||||
return np.exp(-0.5 * ((x - mu) / sd) ** 2) / (sd * math.sqrt(2 * math.pi))
|
||||
|
||||
|
||||
@dataclass
|
||||
class Params:
|
||||
lrr_sd: float
|
||||
baf_sd: float
|
||||
lrr_offset: float = 0.0 # added to every state mean (haploid regions: the region's median LRR)
|
||||
switch_d: float = SWITCH_D
|
||||
eps: float = EPS
|
||||
extra: dict[str, float] = field(default_factory=dict)
|
||||
|
||||
|
||||
def emissions(states: list[State], lrr: F64, baf: F64, pfb: F64, par: Params) -> F64:
|
||||
"""Per-probe, per-state emission likelihood (n × S); missing LRR/BAF contribute 1."""
|
||||
n = len(lrr)
|
||||
out = np.ones((n, len(states)), dtype=np.float64)
|
||||
lok = np.isfinite(lrr)
|
||||
bok = np.isfinite(baf)
|
||||
x = np.where(lok, lrr, 0.0)
|
||||
b = np.clip(np.where(bok, baf, 0.5), 0.0, 1.0)
|
||||
p = np.clip(np.where(np.isfinite(pfb), pfb, 0.5), 0.01, 0.99)
|
||||
for j, st in enumerate(states):
|
||||
sd = 2.0 if st.cn == 0 and not st.mosaic else par.lrr_sd
|
||||
le = (1 - par.eps) * _norm(x, st.lrr_mu + par.lrr_offset, sd) + par.eps / 8.0
|
||||
if st.baf:
|
||||
be = np.zeros(n)
|
||||
for key, mean in st.baf:
|
||||
be += _geno_weights(key, p) * _norm(b, mean, par.baf_sd)
|
||||
be = (1 - par.eps) * be + par.eps
|
||||
else:
|
||||
be = np.ones(n)
|
||||
out[:, j] = np.where(lok, le, 1.0) * np.where(bok, be, 1.0)
|
||||
return out
|
||||
|
||||
|
||||
def switch_probs(pos: npt.NDArray[np.int64], switch_d: float = SWITCH_D) -> F64:
|
||||
"""Switch probability before each probe (the first probe gets 1: it draws from the prior)."""
|
||||
d = np.diff(np.asarray(pos, dtype=np.float64), prepend=np.nan)
|
||||
p = 1.0 - np.exp(-np.maximum(d, 0.0) / switch_d)
|
||||
p[0] = 1.0
|
||||
return np.clip(np.nan_to_num(p, nan=1.0), 1e-12, 1.0)
|
||||
|
||||
|
||||
def forward_backward(e: F64, sw: F64, prior: F64) -> F64:
|
||||
"""Posterior state probabilities (n × S) under the switch model."""
|
||||
n, s = e.shape
|
||||
alpha = np.empty((n, s))
|
||||
c = np.empty(n)
|
||||
a = prior * e[0]
|
||||
c[0] = a.sum() or 1e-300
|
||||
alpha[0] = a / c[0]
|
||||
for t in range(1, n):
|
||||
a = ((1 - sw[t]) * alpha[t - 1] + sw[t] * prior) * e[t]
|
||||
c[t] = a.sum() or 1e-300
|
||||
alpha[t] = a / c[t]
|
||||
beta = np.ones(s)
|
||||
post = np.empty((n, s))
|
||||
post[n - 1] = alpha[n - 1]
|
||||
for t in range(n - 1, 0, -1):
|
||||
eb = e[t] * beta
|
||||
beta = ((1 - sw[t]) * eb + sw[t] * float(prior @ eb)) / c[t]
|
||||
g = alpha[t - 1] * beta
|
||||
post[t - 1] = g / (g.sum() or 1e-300)
|
||||
return post
|
||||
|
||||
|
||||
def viterbi(e: F64, sw: F64, prior: F64) -> npt.NDArray[np.int32]:
|
||||
"""Most likely state path under the switch model (log space)."""
|
||||
n, s = e.shape
|
||||
le = np.log(np.maximum(e, 1e-300))
|
||||
lp = np.log(np.maximum(prior, 1e-300))
|
||||
back = np.empty((n, s), dtype=np.int32)
|
||||
ld = lp + le[0]
|
||||
idx = np.arange(s, dtype=np.int32)
|
||||
for t in range(1, n):
|
||||
stay = np.log(np.maximum(1 - sw[t] + sw[t] * prior, 1e-300)) + ld
|
||||
best = int(np.argmax(ld))
|
||||
jump = math.log(max(sw[t], 1e-300)) + lp + ld[best]
|
||||
take_jump = jump > stay
|
||||
back[t] = np.where(take_jump, best, idx)
|
||||
ld = np.where(take_jump, jump, stay) + le[t]
|
||||
path = np.empty(n, dtype=np.int32)
|
||||
path[-1] = int(np.argmax(ld))
|
||||
for t in range(n - 1, 0, -1):
|
||||
path[t - 1] = back[t, path[t]]
|
||||
return path
|
||||
@@ -0,0 +1,73 @@
|
||||
"""Cyto model covariates — spec CYTO-038, IN-051 (plan §10.1). Pure functions, no I/O.
|
||||
|
||||
GC fraction of windows centred on each probe [LIT: Diskin et al. 2008 NAR 36:e126 (GC-wave), windows 1 kb / 100 kb /
|
||||
1 Mb]: (G+C) / (A+C+G+T) over the window, N bases excluded; NaN when the window has no called base. Population
|
||||
B-allele frequency per probe (PFB) [LIT: Wang et al. 2007 Genome Res 17:1665 (PennCNV)] from cluster-file counts:
|
||||
(n_AB + 2·n_BB) / (2·N); 0.5 when unknown.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
GC_WINDOWS = (1_000, 100_000, 1_000_000)
|
||||
_C, _G, _A, _T = ord("C"), ord("G"), ord("A"), ord("T")
|
||||
|
||||
|
||||
def gc_fractions(
|
||||
seq: npt.NDArray[np.uint8], positions: npt.NDArray[np.int64], windows: tuple[int, ...] = GC_WINDOWS
|
||||
) -> npt.NDArray[np.float32]:
|
||||
"""GC fraction per (probe, window) for 1-based ``positions`` on one contig (upper-case ASCII ``seq``)."""
|
||||
n = len(seq)
|
||||
gc = np.zeros(n + 1, dtype=np.int32)
|
||||
called = np.zeros(n + 1, dtype=np.int32)
|
||||
step = 1 << 24
|
||||
run_gc = run_called = 0
|
||||
for s in range(0, n, step): # chunked cumulative sums keep temporaries small
|
||||
chunk = np.asarray(seq[s : s + step])
|
||||
g = (chunk == _C) | (chunk == _G)
|
||||
c = g | (chunk == _A) | (chunk == _T)
|
||||
cg = np.cumsum(g, dtype=np.int32) + run_gc
|
||||
cc = np.cumsum(c, dtype=np.int32) + run_called
|
||||
gc[s + 1 : s + 1 + len(chunk)] = cg
|
||||
called[s + 1 : s + 1 + len(chunk)] = cc
|
||||
run_gc, run_called = int(cg[-1]), int(cc[-1])
|
||||
out = np.full((len(positions), len(windows)), np.nan, dtype=np.float32)
|
||||
centre = np.asarray(positions, dtype=np.int64) - 1 # 0-based
|
||||
for j, w in enumerate(windows):
|
||||
lo = np.clip(centre - w // 2, 0, n)
|
||||
hi = np.clip(centre - w // 2 + w, 0, n)
|
||||
num = gc[hi] - gc[lo]
|
||||
den = called[hi] - called[lo]
|
||||
ok = den > 0
|
||||
out[ok, j] = num[ok] / den[ok]
|
||||
return out
|
||||
|
||||
|
||||
def population_bfreq(counts: npt.NDArray[np.int32] | None, n_loci: int) -> npt.NDArray[np.float32]:
|
||||
"""PFB from cluster counts (n, 3) = (AA, AB, BB); 0.5 where there are no counts."""
|
||||
if counts is None:
|
||||
return np.full(n_loci, 0.5, dtype=np.float32)
|
||||
c = np.asarray(counts, dtype=np.float64)
|
||||
tot = c.sum(axis=1)
|
||||
with np.errstate(invalid="ignore", divide="ignore"):
|
||||
p = (c[:, 1] + 2 * c[:, 2]) / (2 * tot)
|
||||
return np.where(tot > 0, p, 0.5).astype(np.float32)
|
||||
|
||||
|
||||
def reference_stats(
|
||||
lrr: npt.NDArray[np.float32], baf: npt.NDArray[np.float32]
|
||||
) -> tuple[npt.NDArray[np.float32], npt.NDArray[np.float32], npt.NDArray[np.int32]]:
|
||||
"""Batch reference (plan §10.1) from a samples × probes stack of normal samples: per-probe median LRR, per-probe
|
||||
median heterozygous BAF (samples with 0.2 < BAF < 0.8; 0.5 when none) and the number of such samples."""
|
||||
with np.errstate(all="ignore"):
|
||||
med_lrr = np.nanmedian(np.where(np.isfinite(lrr), lrr, np.nan), axis=0)
|
||||
het = np.where((baf > 0.2) & (baf < 0.8), baf, np.nan)
|
||||
n_het = np.isfinite(het).sum(axis=0).astype(np.int32)
|
||||
het_c = np.where(n_het > 0, np.nanmedian(het, axis=0), 0.5)
|
||||
return (
|
||||
np.nan_to_num(med_lrr, nan=0.0).astype(np.float32),
|
||||
np.nan_to_num(het_c, nan=0.5).astype(np.float32),
|
||||
n_het,
|
||||
)
|
||||
Reference in new issue
Block a user