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
f571ad8a2e
commit
724e87dc66
8 files changed
+1452
No files matched your search
@@ -8,6 +8,10 @@ Generated from `src/simurg/core/errors.py` — do not edit by hand.
|
||||
| `E_BPM_CSV_MISMATCH` | error | BPM and CSV manifests describe different loci | IN-042 |
|
||||
| `E_BPM_FORMAT` | error | BPM manifest is malformed or unsupported | IN-040 |
|
||||
| `E_BUILD_MISMATCH` | error | Manifest genome build differs from the FASTA build | IN-063 |
|
||||
| `E_CN_FAILED` | error | CN calling failed for the sample | CYTO-040 |
|
||||
| `E_CN_SEX_UNRESOLVED` | error | Sample sex unresolved and --unknown-sex-policy error | CYTO-024 |
|
||||
| `E_CN_SIDECAR_MISSING` | error | cn_call_manifest.json not found in the cn call output | CYTO-032 |
|
||||
| `E_CN_VCF_UNINDEXED` | error | CNV VCFs must be bgzipped and indexed for annotate | CYTO-032 |
|
||||
| `E_CONFIG_KEY` | error | Unknown key in --config YAML | CLI-019 |
|
||||
| `E_CONFIG_PARSE` | error | --config file is not valid YAML | CLI-019 |
|
||||
| `E_CSV_MANIFEST_FORMAT` | error | CSV manifest is malformed | IN-041 |
|
||||
@@ -32,6 +36,8 @@ Generated from `src/simurg/core/errors.py` — do not edit by hand.
|
||||
| `E_METH_NO_CONTROLS` | error | Missing control probes in manifest | MET-002 |
|
||||
| `E_METH_NO_IDATS` | error | No IDATs found | MET-004 |
|
||||
| `E_METH_PREP_FAILED` | error | Prep failed (sample) | MET-010 |
|
||||
| `E_MODEL_MANIFEST_MISMATCH` | error | Model was built for a different manifest probe set | IN-051 |
|
||||
| `E_MODEL_VERSION` | error | Model schema is too old; rebuild it | IN-051 |
|
||||
| `E_NO_INPUT` | error | Neither an input folder nor a sample sheet was given | CLI-006 |
|
||||
| `E_NO_INPUT_FOLDER` | error | Sample sheet has no RootFolder and no folder was given | CLI-006 |
|
||||
| `E_OUTPUT_EXISTS` | error | Outputs already exist; use --force or --resume | CLI-010 |
|
||||
@@ -59,7 +65,11 @@ Generated from `src/simurg/core/errors.py` — do not edit by hand.
|
||||
| `E_SHEET_NO_KEY_COLUMNS` | error | Sample sheet lacks SentrixBarcode_A/SentrixPosition_A | IN-002 |
|
||||
| `E_SHEET_RESERVED` | error | Sample sheet column collides with an output column | IN-005 |
|
||||
| `E_SHEET_TOO_LONG` | error | Sample sheet name or value longer than 500 characters | IN-006 |
|
||||
| `I_CN_SEX_RESOLVED` | info | Sample sex updated from U by the caller | CYTO-024 |
|
||||
| `I_UNLISTED_IDATS` | info | IDAT pairs not listed in the sample sheet are ignored | IN-016 |
|
||||
| `W_CN_NO_FASTA` | warning | No FASTA: VCF REF padding bases written as N | CYTO-029 |
|
||||
| `W_CN_PROBE_FLOOR` | warning | Minimum probe count below 10 raised to 10 | CYTO-008 |
|
||||
| `W_CN_THRESHOLD_EXCEEDS_GENOME` | warning | Size threshold exceeds the effective genome; no events | CYTO-009 |
|
||||
| `W_CTRL_DEFS_DIFFER` | warning | BPM and CSV control definitions differ; the BPM wins | QC-001 |
|
||||
| `W_CTRL_MISSING` | warning | Controls needed by a metric are absent; metric not evaluable | QC-007 |
|
||||
| `W_CTRL_ZERO_DIVISOR` | warning | Control metric has a zero denominator (NaN) | QC-008 |
|
||||
|
||||
@@ -666,6 +666,107 @@ def qc_controls_cmd(
|
||||
run_guarded(rep, lambda: run_qc_controls(opts, rep))
|
||||
|
||||
|
||||
# --------------------------------------------------------------------------------------------------------
|
||||
cn_app = typer.Typer(help="Genome-wide copy number, LOH and mosaicism (CYTO-001).")
|
||||
app.add_typer(cn_app, name="cn")
|
||||
|
||||
|
||||
@cn_app.command("build-model")
|
||||
def cn_build_model_cmd(
|
||||
bpm: Annotated[Path, typer.Option("--bpm", help="BPM manifest (required).")],
|
||||
fasta: Annotated[Path, typer.Option("--fasta", help="Reference FASTA of the manifest's build (required).")],
|
||||
out: Annotated[Path, typer.Option("--out", help="Output model file (.simurg-cyto).")],
|
||||
egt: Annotated[Path | None, typer.Option("--egt", help="Cluster file: population B-allele frequencies.")] = None,
|
||||
reference_gtcs: Annotated[
|
||||
Path | None, typer.Option("--reference-gtcs", help="GTCs of normal samples (>= 24) for the batch reference.")
|
||||
] = None,
|
||||
max_reference: Annotated[int, typer.Option("--max-reference", min=1, help="Random subset size.")] = 100,
|
||||
log_json: LogJson = False,
|
||||
debug: Debug = False,
|
||||
quiet: Quiet = False,
|
||||
) -> None:
|
||||
"""Cyto model: GC covariates, PFB and an optional batch reference per probe (CYTO-038, IN-051)."""
|
||||
from simurg.pipeline.cn import run_build_model
|
||||
|
||||
rep = _reporter("cn build-model", out.parent, log_json, debug, quiet)
|
||||
run_guarded(rep, lambda: run_build_model(bpm, fasta, out, rep, egt, "auto", reference_gtcs, max_reference))
|
||||
|
||||
|
||||
def _kb(v: int) -> int:
|
||||
if not 0 <= v <= 250_000:
|
||||
raise typer.BadParameter("sizes are 0-250000 kb")
|
||||
return v * 1000
|
||||
|
||||
|
||||
@cn_app.command("call")
|
||||
def cn_call_cmd(
|
||||
bpm: Annotated[Path, typer.Option("--bpm", help="BPM manifest (required).")],
|
||||
model: Annotated[Path, typer.Option("--model", help="Cyto model from cn build-model (required).")],
|
||||
out: Annotated[Path, typer.Option("--out", help="Output folder (required).")],
|
||||
gtcs: Annotated[Path | None, typer.Option("--gtcs", help="Folder searched recursively for GTC files.")] = None,
|
||||
summary: Annotated[
|
||||
Path | None, typer.Option("--summary", help="sample_summary.tsv (default: next to GTCs).")
|
||||
] = None,
|
||||
fasta: Annotated[Path | None, typer.Option("--fasta", help="FASTA for REF padding bases.")] = None,
|
||||
min_del_probes: Annotated[int, typer.Option("--min-del-probes", min=0, max=250000)] = 10,
|
||||
min_del_size: Annotated[int, typer.Option("--min-del-size", help="kb")] = 25,
|
||||
min_dup_probes: Annotated[int, typer.Option("--min-dup-probes", min=0, max=250000)] = 10,
|
||||
min_dup_size: Annotated[int, typer.Option("--min-dup-size", help="kb")] = 50,
|
||||
min_loh_probes: Annotated[int, typer.Option("--min-loh-probes", min=0, max=250000)] = 500,
|
||||
min_loh_size: Annotated[int, typer.Option("--min-loh-size", help="kb")] = 3000,
|
||||
min_gainloh_probes: Annotated[int, typer.Option("--min-gainloh-probes", min=0, max=250000)] = 300,
|
||||
min_gainloh_size: Annotated[int, typer.Option("--min-gainloh-size", help="kb")] = 500,
|
||||
min_mdel_probes: Annotated[int, typer.Option("--min-mdel-probes", min=0, max=250000)] = 10,
|
||||
min_mdel_size: Annotated[int, typer.Option("--min-mdel-size", help="kb")] = 25,
|
||||
min_mdup_probes: Annotated[int, typer.Option("--min-mdup-probes", min=0, max=250000)] = 10,
|
||||
min_mdup_size: Annotated[int, typer.Option("--min-mdup-size", help="kb")] = 50,
|
||||
min_mloh_probes: Annotated[int, typer.Option("--min-mloh-probes", min=0, max=250000)] = 500,
|
||||
min_mloh_size: Annotated[int, typer.Option("--min-mloh-size", help="kb")] = 3000,
|
||||
smoothing: Annotated[int, typer.Option("--smoothing", min=0, max=1000)] = 5,
|
||||
max_mosaic_fraction: Annotated[float, typer.Option("--max-mosaic-fraction")] = 1.0,
|
||||
min_mosaic_fraction: Annotated[float, typer.Option("--min-mosaic-fraction")] = 0.10,
|
||||
filter_qual: Annotated[int, typer.Option("--filter-qual", min=0, max=60)] = 10,
|
||||
cn_cap: Annotated[int, typer.Option("--cn-cap", min=3, max=8)] = 4,
|
||||
gc_correct: Annotated[bool, typer.Option("--gc-correct/--no-gc-correct")] = True,
|
||||
unknown_sex_policy: Annotated[str, typer.Option("--unknown-sex-policy")] = "skip-sex-chromosomes",
|
||||
gt_mode: Annotated[str, typer.Option("--gt-mode", help="vcf|compat")] = "vcf",
|
||||
svlen_mode: Annotated[str, typer.Option("--svlen-mode", help="signed|positive")] = "signed",
|
||||
mask_bed: Annotated[Path | None, typer.Option("--mask-bed", help="Mask BED (default: built-in).")] = None,
|
||||
no_bgzip: Annotated[bool, typer.Option("--no-bgzip", help="Plain VCF (cn annotate then refuses it).")] = False,
|
||||
threads: ThreadsOpt = 0,
|
||||
force: ForceOpt = False,
|
||||
log_json: LogJson = False,
|
||||
debug: Debug = False,
|
||||
quiet: Quiet = False,
|
||||
) -> None:
|
||||
"""GTCs -> per-sample CNV/LOH/mosaic VCF, segments and summary (CLI-028)."""
|
||||
from simurg.cn.events import Thresholds
|
||||
from simurg.pipeline.cn import CnCallOptions, run_cn_call
|
||||
|
||||
if unknown_sex_policy not in ("skip-sex-chromosomes", "female", "male", "error"):
|
||||
raise typer.BadParameter("--unknown-sex-policy: skip-sex-chromosomes|female|male|error")
|
||||
if gt_mode not in ("vcf", "compat") or svlen_mode not in ("signed", "positive"):
|
||||
raise typer.BadParameter("--gt-mode vcf|compat, --svlen-mode signed|positive")
|
||||
if not 0 < max_mosaic_fraction <= 1 or not 0 <= min_mosaic_fraction < 1:
|
||||
raise typer.BadParameter("mosaic fractions: 0 < max <= 1, 0 <= min < 1")
|
||||
th = Thresholds(
|
||||
probes={"del": min_del_probes, "dup": min_dup_probes, "loh": min_loh_probes, "gainloh": min_gainloh_probes,
|
||||
"mdel": min_mdel_probes, "mdup": min_mdup_probes, "mloh": min_mloh_probes},
|
||||
size_bp={"del": _kb(min_del_size), "dup": _kb(min_dup_size), "loh": _kb(min_loh_size),
|
||||
"gainloh": _kb(min_gainloh_size), "mdel": _kb(min_mdel_size), "mdup": _kb(min_mdup_size),
|
||||
"mloh": _kb(min_mloh_size)},
|
||||
max_mosaic_fraction=max_mosaic_fraction, min_mosaic_fraction=min_mosaic_fraction, filter_qual=filter_qual,
|
||||
) # fmt: skip
|
||||
rep = _reporter("cn call", out, log_json, debug, quiet)
|
||||
opts = CnCallOptions(
|
||||
out=out, bpm=bpm, model=model, gtcs=gtcs, summary=summary, fasta=fasta, thresholds=th, smoothing=smoothing,
|
||||
cn_cap=cn_cap, gc_correct=gc_correct, unknown_sex_policy=unknown_sex_policy, gt_mode=gt_mode,
|
||||
svlen_mode=svlen_mode, mask_bed=mask_bed, bgzip=not no_bgzip, threads=threads, force=force,
|
||||
argv=sys.argv[1:],
|
||||
) # fmt: skip
|
||||
run_guarded(rep, lambda: run_cn_call(opts, rep))
|
||||
|
||||
|
||||
# --------------------------------------------------------------------------------------------------------
|
||||
methyl_app = typer.Typer(help="Methylation arrays: control QC, detection, beta values (MET-001).")
|
||||
app.add_typer(methyl_app, name="methyl")
|
||||
|
||||
@@ -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,
|
||||
)
|
||||
@@ -135,6 +135,22 @@ CATALOGUE: dict[str, MessageDef] = dict(
|
||||
_m("W_QC_PEDIGREE", "warning", "Parent1/Parent2/Replicate value does not name a sample in the set", "IN-018"),
|
||||
_m("W_QC_NO_BAF", "warning", "GTC has no B-allele frequencies; contamination not estimated", "QC-064"),
|
||||
_m("W_QC_TOO_FEW", "warning", "Too few samples or sites for a cohort QC statistic", "QC-026"),
|
||||
# --- cytogenetics CNV/LOH (CYTO-001 .. CYTO-040) -------------------------------------------------
|
||||
_m("E_MODEL_MANIFEST_MISMATCH", "error", "Model was built for a different manifest probe set", "IN-051"),
|
||||
_m("E_MODEL_VERSION", "error", "Model schema is too old; rebuild it", "IN-051"),
|
||||
_m("W_CN_PROBE_FLOOR", "warning", "Minimum probe count below 10 raised to 10", "CYTO-008"),
|
||||
_m(
|
||||
"W_CN_THRESHOLD_EXCEEDS_GENOME",
|
||||
"warning",
|
||||
"Size threshold exceeds the effective genome; no events",
|
||||
"CYTO-009",
|
||||
),
|
||||
_m("I_CN_SEX_RESOLVED", "info", "Sample sex updated from U by the caller", "CYTO-024"),
|
||||
_m("E_CN_SEX_UNRESOLVED", "error", "Sample sex unresolved and --unknown-sex-policy error", "CYTO-024"),
|
||||
_m("W_CN_NO_FASTA", "warning", "No FASTA: VCF REF padding bases written as N", "CYTO-029"),
|
||||
_m("E_CN_SIDECAR_MISSING", "error", "cn_call_manifest.json not found in the cn call output", "CYTO-032"),
|
||||
_m("E_CN_VCF_UNINDEXED", "error", "CNV VCFs must be bgzipped and indexed for annotate", "CYTO-032"),
|
||||
_m("E_CN_FAILED", "error", "CN calling failed for the sample", "CYTO-040"),
|
||||
# --- methylation control QC (MET-001 .. MET-018) -------------------------------------------------
|
||||
_m("E_METH_NO_ASSAY", "error", "Probe section not found", "MET-002"),
|
||||
_m("E_METH_NO_CONTROLS", "error", "Missing control probes in manifest", "MET-002"),
|
||||
|
||||
@@ -0,0 +1,779 @@
|
||||
"""``simurg cn build-model`` / ``simurg cn call`` — genome-wide CNV/LOH/mosaic calling (spec CYTO-001 … CYTO-040).
|
||||
|
||||
I/O and orchestration only; the numerics live in ``simurg.cn``.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import hashlib
|
||||
import io
|
||||
import json
|
||||
import math
|
||||
import zipfile
|
||||
from dataclasses import dataclass, field
|
||||
from pathlib import Path
|
||||
from typing import Any
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
import pyarrow as pa
|
||||
import pyarrow.parquet as pq
|
||||
|
||||
from simurg import __version__
|
||||
from simurg.core.errors import ExitCode, SimurgError
|
||||
from simurg.core.log import Reporter
|
||||
from simurg.formats.bpm import Manifest, read_bpm
|
||||
from simurg.formats.fasta import ContigMap, Genome
|
||||
|
||||
MODEL_SCHEMA = 1
|
||||
F64 = npt.NDArray[np.float64]
|
||||
|
||||
|
||||
def probe_set_hash(manifest: Manifest) -> str:
|
||||
h = hashlib.sha256()
|
||||
for name, chrom, pos in zip(manifest.names, manifest.chroms, manifest.map_infos.tolist(), strict=True):
|
||||
h.update(f"{name}\t{chrom}\t{pos}\n".encode())
|
||||
return h.hexdigest()
|
||||
|
||||
|
||||
@dataclass
|
||||
class CytoModel:
|
||||
header: dict[str, Any]
|
||||
names: list[str]
|
||||
gc: npt.NDArray[np.float32] # (n, 3): 1 kb, 100 kb, 1 Mb
|
||||
pfb: npt.NDArray[np.float32]
|
||||
ref_lrr: npt.NDArray[np.float32] | None = None # batch reference: per-probe median LRR
|
||||
ref_lrr_f: npt.NDArray[np.float32] | None = None # same, X/Y from female reference samples
|
||||
ref_lrr_m: npt.NDArray[np.float32] | None = None # same, X/Y from male reference samples
|
||||
ref_het: npt.NDArray[np.float32] | None = None # per-probe median heterozygous BAF
|
||||
|
||||
|
||||
def run_build_model(
|
||||
bpm: Path,
|
||||
fasta: Path,
|
||||
out: Path,
|
||||
rep: Reporter,
|
||||
egt: Path | None = None,
|
||||
contig_alias: str = "auto",
|
||||
reference_gtcs: Path | None = None,
|
||||
max_reference: int = 100,
|
||||
) -> int:
|
||||
"""CYTO-038: GC covariates (and PFB when a cluster file is given) per manifest probe → ``.simurg-cyto``."""
|
||||
from simurg.cn.model import GC_WINDOWS, gc_fractions, population_bfreq
|
||||
|
||||
man = read_bpm(bpm)
|
||||
genome = Genome(fasta)
|
||||
cmap = ContigMap([e.name for e in genome.entries], contig_alias)
|
||||
chroms = np.array(man.chroms)
|
||||
pos = man.map_infos.astype(np.int64)
|
||||
gc = np.full((man.num_loci, len(GC_WINDOWS)), np.nan, dtype=np.float32)
|
||||
for chrom in sorted(set(man.chroms)):
|
||||
contig = cmap.get(chrom)
|
||||
if contig is None:
|
||||
continue
|
||||
sel = np.flatnonzero((chroms == chrom) & (pos > 0))
|
||||
if len(sel):
|
||||
gc[sel] = gc_fractions(genome._contig_array(contig), pos[sel])
|
||||
rep.info(f"GC windows for chromosome {chrom}: {len(sel)} probes")
|
||||
pfb = population_bfreq(None, man.num_loci)
|
||||
egt_name = ""
|
||||
if egt is not None:
|
||||
from simurg.formats.egt import read_egt
|
||||
from simurg.genotype.clusters import align_clusters
|
||||
|
||||
clusters = align_clusters(read_egt(egt), man, allow_missing=True)
|
||||
pfb = population_bfreq(np.where(clusters.present[:, None], clusters.n, 0), man.num_loci)
|
||||
egt_name = Path(egt).name
|
||||
ref_cols: dict[str, Any] = {}
|
||||
n_ref = 0
|
||||
if reference_gtcs is not None:
|
||||
from simurg.cn.model import reference_stats
|
||||
from simurg.formats.gtc import read_gtc
|
||||
|
||||
paths = sorted(reference_gtcs.rglob("*.gtc"))
|
||||
if len(paths) > max_reference:
|
||||
idx = np.random.default_rng(42).choice(len(paths), max_reference, replace=False)
|
||||
paths = [paths[i] for i in sorted(idx)]
|
||||
if len(paths) < 24:
|
||||
rep.warn("W_QC_TOO_FEW", f"batch reference from only {len(paths)} GTCs (plan: >= 24)")
|
||||
lrr = np.empty((len(paths), man.num_loci), dtype=np.float32)
|
||||
baf = np.empty((len(paths), man.num_loci), dtype=np.float32)
|
||||
sexes = []
|
||||
for i, p in enumerate(paths):
|
||||
g = read_gtc(p)
|
||||
if len(g.logr_ratios) != man.num_loci:
|
||||
raise SimurgError("E_MODEL_MANIFEST_MISMATCH", f"{p} has {len(g.logr_ratios)} loci, the manifest "
|
||||
f"{man.num_loci}") # fmt: skip
|
||||
lrr[i], baf[i] = g.logr_ratios, g.b_allele_freqs
|
||||
sexes.append(g.gender.decode(errors="replace")[:1].upper() if g.gender else "U")
|
||||
ref_lrr, ref_het, n_het = reference_stats(lrr, baf)
|
||||
# X and Y from same-sex reference samples only (CYTO-005 baselines differ by sex)
|
||||
sx = np.isin(chroms, ["X", "XY", "Y"])
|
||||
ref_by_sex = {}
|
||||
for sex in ("F", "M"):
|
||||
rows = np.array([s_ == sex for s_ in sexes])
|
||||
col = ref_lrr.copy()
|
||||
if rows.sum() >= 3:
|
||||
col[sx] = reference_stats(lrr[rows][:, sx], baf[rows][:, sx])[0]
|
||||
ref_by_sex[sex] = col
|
||||
n_f, n_m = sexes.count("F"), sexes.count("M")
|
||||
rep.info(f"batch reference sexes: {n_f} F, {n_m} M, {len(sexes) - n_f - n_m} U")
|
||||
ref_cols = {"ref_lrr": ref_lrr, "ref_lrr_f": ref_by_sex["F"], "ref_lrr_m": ref_by_sex["M"],
|
||||
"ref_het_baf": ref_het, "ref_n_het": n_het} # fmt: skip
|
||||
n_ref = len(paths)
|
||||
rep.info(f"batch reference from {n_ref} GTC(s)")
|
||||
header = {
|
||||
"schema": "simurg.cyto_model",
|
||||
"schema_version": MODEL_SCHEMA,
|
||||
"software": f"SimurgArray {__version__}",
|
||||
"manifest": man.manifest_name,
|
||||
"probe_set_sha256": probe_set_hash(man),
|
||||
"n_probes": man.num_loci,
|
||||
"genome_build": genome.build,
|
||||
"fasta_sha256": genome.sha256,
|
||||
"gc_windows_bp": list(GC_WINDOWS),
|
||||
"pfb_source": f"cluster file {egt_name}" if egt_name else "none (0.5)",
|
||||
"reference_samples": n_ref,
|
||||
}
|
||||
table = pa.table({
|
||||
"name": man.names,
|
||||
"gc_1kb": gc[:, 0], "gc_100kb": gc[:, 1], "gc_1mb": gc[:, 2],
|
||||
"pfb": pfb,
|
||||
**ref_cols,
|
||||
}) # fmt: skip
|
||||
buf = io.BytesIO()
|
||||
pq.write_table(table, buf, compression="zstd")
|
||||
out.parent.mkdir(parents=True, exist_ok=True)
|
||||
tmp = out.with_suffix(out.suffix + ".part")
|
||||
with zipfile.ZipFile(tmp, "w", compression=zipfile.ZIP_STORED) as z:
|
||||
z.writestr("header.json", json.dumps(header, indent=1))
|
||||
z.writestr("probes.parquet", buf.getvalue())
|
||||
tmp.replace(out)
|
||||
rep.info(f"cyto model for {man.num_loci} probes written to {out}")
|
||||
return ExitCode.OK
|
||||
|
||||
|
||||
def read_model(path: Path, manifest: Manifest) -> CytoModel:
|
||||
"""IN-051: check schema and probe set."""
|
||||
with zipfile.ZipFile(path) as z:
|
||||
header = json.loads(z.read("header.json"))
|
||||
table = pq.read_table(io.BytesIO(z.read("probes.parquet")))
|
||||
if int(header.get("schema_version", 0)) < MODEL_SCHEMA:
|
||||
raise SimurgError("E_MODEL_VERSION", f"{path}: cyto model schema {header.get('schema_version')} is too old",
|
||||
hint="rebuild it with simurg cn build-model") # fmt: skip
|
||||
if header.get("probe_set_sha256") != probe_set_hash(manifest):
|
||||
raise SimurgError("E_MODEL_MANIFEST_MISMATCH", f"{path} was built for manifest {header.get('manifest')!r}, "
|
||||
f"not {manifest.manifest_name!r} (probe set differs)") # fmt: skip
|
||||
gc = np.column_stack([table.column(c).to_numpy(zero_copy_only=False) for c in ("gc_1kb", "gc_100kb", "gc_1mb")])
|
||||
|
||||
def col(c: str) -> npt.NDArray[np.float32] | None:
|
||||
return table.column(c).to_numpy(zero_copy_only=False).astype(np.float32) if c in table.column_names else None
|
||||
|
||||
pfb = col("pfb")
|
||||
assert pfb is not None
|
||||
return CytoModel(header, table.column("name").to_pylist(), gc.astype(np.float32), pfb, ref_lrr=col("ref_lrr"),
|
||||
ref_lrr_f=col("ref_lrr_f"), ref_lrr_m=col("ref_lrr_m"), ref_het=col("ref_het_baf")) # fmt: skip
|
||||
|
||||
|
||||
# ================================================================================================ cn call
|
||||
PAR = { # 1-based inclusive pseudo-autosomal regions on X
|
||||
"GRCh38": ((10_001, 2_781_479), (155_701_383, 156_030_895)),
|
||||
"GRCh37": ((60_001, 2_699_520), (154_931_044, 155_260_560)),
|
||||
}
|
||||
CALL_CHROMS = [str(i) for i in range(1, 23)] + ["X", "Y"]
|
||||
|
||||
|
||||
@dataclass
|
||||
class CnCallOptions:
|
||||
out: Path
|
||||
bpm: Path
|
||||
model: Path
|
||||
gtcs: Path | None = None
|
||||
sheet: Path | None = None
|
||||
summary: Path | None = None
|
||||
fasta: Path | None = None
|
||||
thresholds: Any = None # simurg.cn.events.Thresholds
|
||||
smoothing: int = 5
|
||||
cn_cap: int = 4
|
||||
gc_correct: bool = True
|
||||
unknown_sex_policy: str = "skip-sex-chromosomes"
|
||||
gt_mode: str = "vcf"
|
||||
svlen_mode: str = "signed"
|
||||
mask_bed: Path | None = None
|
||||
bgzip: bool = True
|
||||
threads: int = 0
|
||||
force: bool = False
|
||||
argv: list[str] = field(default_factory=list)
|
||||
|
||||
|
||||
def _builtin_file(name: str) -> bytes:
|
||||
from importlib import resources
|
||||
|
||||
return resources.files("simurg.cn").joinpath("data", name).read_bytes()
|
||||
|
||||
|
||||
def read_mask(path: Path | None, build: str) -> dict[str, list[tuple[int, int]]]:
|
||||
"""CYTO-017: BED (0-based half-open) → per-chromosome intervals as 1-based inclusive."""
|
||||
text = path.read_text() if path is not None else _builtin_file(f"mask_{build}.bed").decode()
|
||||
out: dict[str, list[tuple[int, int]]] = {}
|
||||
for line in text.splitlines():
|
||||
if not line.strip() or line.startswith(("#", "track", "browser")):
|
||||
continue
|
||||
c, s, e = line.split("\t")[:3]
|
||||
c = c[3:] if c.lower().startswith("chr") else c
|
||||
out.setdefault(c, []).append((int(s) + 1, int(e)))
|
||||
return out
|
||||
|
||||
|
||||
def moving_median(x: F64, k: int) -> F64:
|
||||
"""CYTO-004: median over 2k+1 consecutive values, truncated (never wrapped) at the ends; NaN ignored."""
|
||||
if k <= 0 or len(x) == 0:
|
||||
return x.copy()
|
||||
from numpy.lib.stride_tricks import sliding_window_view
|
||||
|
||||
pad = np.concatenate([np.full(k, np.nan), x, np.full(k, np.nan)])
|
||||
with np.errstate(all="ignore"):
|
||||
import warnings
|
||||
|
||||
with warnings.catch_warnings():
|
||||
warnings.simplefilter("ignore", RuntimeWarning)
|
||||
med = np.nanmedian(sliding_window_view(pad, 2 * k + 1), axis=1)
|
||||
return np.where(np.isfinite(x), med, np.nan)
|
||||
|
||||
|
||||
@dataclass
|
||||
class SampleInput:
|
||||
sample_id: str
|
||||
sample_name: str
|
||||
gtc: str
|
||||
sex: str # M | F | U
|
||||
logr_dev: float
|
||||
sex_source: str
|
||||
|
||||
|
||||
def _robust_sd(x: F64, floor: float) -> float:
|
||||
x = x[np.isfinite(x)]
|
||||
if not len(x):
|
||||
return floor
|
||||
return max(floor, float(1.4826 * np.median(np.abs(x - np.median(x)))))
|
||||
|
||||
|
||||
def _resolve_sex(si: SampleInput, chroms: npt.NDArray[np.str_], lrr: F64, geno: npt.NDArray[np.uint8],
|
||||
policy: str) -> tuple[str, str, list[str]]: # fmt: skip
|
||||
"""CYTO-024."""
|
||||
if si.sex in ("M", "F"):
|
||||
return si.sex, si.sex_source, []
|
||||
x = chroms == "X"
|
||||
called = geno[x] > 0
|
||||
het = float((geno[x] == 2).sum() / called.sum()) if called.any() else math.nan
|
||||
med = float(np.nanmedian(lrr[x])) if x.any() else math.nan
|
||||
if math.isfinite(het) and math.isfinite(med):
|
||||
if het < 0.03 and med < -0.15:
|
||||
return "M", "resolved", []
|
||||
if het > 0.08 and med > -0.1:
|
||||
return "F", "resolved", []
|
||||
if policy == "female":
|
||||
return "F", "assumed", ["sex_assumed"]
|
||||
if policy == "male":
|
||||
return "M", "assumed", ["sex_assumed"]
|
||||
if policy == "error":
|
||||
raise SimurgError("E_CN_SEX_UNRESOLVED", f"{si.sample_id}: sex could not be resolved")
|
||||
return "U", "unresolved", ["sex_unresolved"]
|
||||
|
||||
|
||||
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
|
||||
from simurg.cn.events import (
|
||||
annotate_stats,
|
||||
apply_thresholds,
|
||||
effective_size,
|
||||
extract,
|
||||
median_spacing,
|
||||
post_filter,
|
||||
set_filters,
|
||||
)
|
||||
from simurg.cn.gcwave import gc_correct
|
||||
from simurg.formats.gtc import read_gtc
|
||||
|
||||
man_chroms: npt.NDArray[np.str_] = ctx["chroms"]
|
||||
pos_all: npt.NDArray[np.int64] = ctx["pos"]
|
||||
model: CytoModel = ctx["model"]
|
||||
th = opts.thresholds
|
||||
g = read_gtc(si.gtc)
|
||||
lrr = g.logr_ratios.astype(np.float64)
|
||||
baf = g.b_allele_freqs.astype(np.float64)
|
||||
geno = g.genotypes
|
||||
lrr[~np.isfinite(lrr)] = np.nan
|
||||
baf[~np.isfinite(baf)] = np.nan
|
||||
sex, sex_source, flags = _resolve_sex(si, man_chroms, lrr - float(np.nanmedian(lrr[ctx["auto"]])), geno,
|
||||
opts.unknown_sex_policy) # fmt: skip
|
||||
ref = {"F": model.ref_lrr_f, "M": model.ref_lrr_m}.get(sex) if model.ref_lrr_f is not None else None
|
||||
if ref is None:
|
||||
ref = model.ref_lrr
|
||||
if ref is not None:
|
||||
lrr = lrr - ref
|
||||
if model.ref_het is not None:
|
||||
het = (baf > 0.1) & (baf < 0.9)
|
||||
baf = np.where(het, np.clip(baf - (model.ref_het - 0.5), 0.0, 1.0), baf)
|
||||
auto = ctx["auto"]
|
||||
gcfit = None
|
||||
if opts.gc_correct:
|
||||
lrr, gcfit = gc_correct(lrr, model.gc, auto)
|
||||
lrr = lrr - float(np.nanmedian(lrr[auto]))
|
||||
logr_dev = si.logr_dev if math.isfinite(si.logr_dev) else float(np.nanstd(lrr[auto]))
|
||||
noisy = logr_dev > 0.2
|
||||
whole_chrom_only = logr_dev > 0.3
|
||||
if noisy:
|
||||
flags.append("high_noise")
|
||||
if whole_chrom_only:
|
||||
flags.append("whole_chromosome_only")
|
||||
pfb = model.pfb.astype(np.float64)
|
||||
lrr_s = np.full_like(lrr, np.nan)
|
||||
for sel in ctx["by_chrom"].values():
|
||||
lrr_s[sel] = moving_median(lrr[sel], opts.smoothing)
|
||||
# emission SD from the unsmoothed LRR: smoothed values are autocorrelated, their own SD would overstate evidence
|
||||
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)
|
||||
dip, hap = hmm.diploid_states(), hmm.haploid_states()
|
||||
counts: dict[str, int] = {}
|
||||
events: list[Any] = []
|
||||
segments: list[dict[str, Any]] = []
|
||||
ploidy_num = 0.0
|
||||
ploidy_den = 0
|
||||
par = PAR.get(ctx["build"], PAR["GRCh38"])
|
||||
for chrom in CALL_CHROMS:
|
||||
if chrom in ("X", "Y") and sex == "U":
|
||||
continue
|
||||
if chrom == "Y" and sex == "F":
|
||||
continue
|
||||
sel = ctx["by_chrom"].get(chrom)
|
||||
if sel is None or len(sel) == 0:
|
||||
continue
|
||||
pos = pos_all[sel]
|
||||
in_par = np.zeros(len(sel), dtype=np.bool_)
|
||||
if chrom == "X":
|
||||
for a, b in par:
|
||||
in_par |= (pos >= a) & (pos <= b)
|
||||
regions = [(np.ones(len(sel), dtype=np.bool_), False)]
|
||||
if chrom == "X" and sex == "M":
|
||||
regions = [(in_par, False), (~in_par, True)]
|
||||
elif chrom == "Y":
|
||||
regions = [(np.ones(len(sel), dtype=np.bool_), True)]
|
||||
m_sp = median_spacing(pos)
|
||||
for rmask, haploid in regions:
|
||||
if rmask.sum() < 2:
|
||||
continue
|
||||
idx = sel[rmask]
|
||||
p, x_raw, b, gt = pos_all[idx], lrr[idx], baf[idx], geno[idx]
|
||||
x_seg = lrr_s[idx]
|
||||
states = hap if haploid else dip
|
||||
prior = np.array([s.prior for s in states])
|
||||
prior /= prior.sum()
|
||||
offset = float(np.nanmedian(x_raw)) if haploid else 0.0
|
||||
if not math.isfinite(offset):
|
||||
offset = 0.0
|
||||
prm = hmm.Params(seg_sd, baf_sd, lrr_offset=offset)
|
||||
e = hmm.emissions(states, x_seg, np.full(len(b), np.nan) if haploid else b, pfb[idx], prm)
|
||||
sw = hmm.switch_probs(p)
|
||||
path = hmm.viterbi(e, sw, prior)
|
||||
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)
|
||||
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_])]
|
||||
segments.append({"chrom": chrom, "start": int(p[s_]), "end": int(p[e_ - 1]), "probes": e_ - s_,
|
||||
"state": st.name, "call": st.call, "cn": st.cn, "haploid": haploid,
|
||||
"posterior": float(post[s_:e_, int(path[s_])].mean())}) # fmt: skip
|
||||
base_cn = 1 if haploid else 2
|
||||
eff = effective_size(p, m_sp)
|
||||
if whole_chrom_only:
|
||||
evs = [ev for ev in evs if ev.effective_size >= 0.9 * eff]
|
||||
ploidy_num += base_cn * eff + sum((ev.cn - base_cn) * ev.effective_size for ev in evs)
|
||||
ploidy_den += eff
|
||||
events.extend((ev, haploid) for ev in evs)
|
||||
kept = apply_thresholds([ev for ev, _ in events], th, counts)
|
||||
kept = _mask_events(kept, ctx["mask"], counts)
|
||||
set_filters(kept, th, noisy)
|
||||
hap_of = {id(ev): h for ev, h in events}
|
||||
row = {
|
||||
"sample_id": si.sample_id, "sample_name": si.sample_name, "status": "OK", "logr_dev": logr_dev,
|
||||
"overall_ploidy": ploidy_num / max(ploidy_den, 1), "gc_correct": bool(gcfit and gcfit.applied),
|
||||
"gcwf": gcfit.gcwf if gcfit else math.nan, "sex_used": sex, "sex_source": sex_source,
|
||||
"n_del": sum(e.call == "DEL" for e in kept), "n_dup": sum(e.call == "DUP" for e in kept),
|
||||
"n_loh": sum(e.call == "LOH" for e in kept), "n_gainloh": sum(e.call == "GAINLOH" for e in kept),
|
||||
"n_mosaic": sum(e.mosaic for e in kept), "n_filtered": sum(counts.values()),
|
||||
"qc_flags": ";".join(flags),
|
||||
} # fmt: skip
|
||||
per_chr: dict[str, int] = {}
|
||||
for ev in kept:
|
||||
if ev.call in ("DEL", "DUP", "LOH"):
|
||||
per_chr[ev.chrom] = per_chr.get(ev.chrom, 0) + 1
|
||||
excess = [c for c, n in per_chr.items() if n > 10]
|
||||
if excess:
|
||||
row["qc_flags"] = ";".join([*flags, "excess_events:" + ",".join(excess)])
|
||||
row["lrr_sd_segmentation"] = seg_sd
|
||||
row["baf_sd"] = baf_sd
|
||||
return {"row": row, "events": kept, "haploid": hap_of, "segments": segments, "counts": counts}
|
||||
|
||||
|
||||
def _runs(path: npt.NDArray[np.int32]) -> list[tuple[int, int]]:
|
||||
if len(path) == 0:
|
||||
return []
|
||||
br = np.flatnonzero(np.diff(path) != 0) + 1
|
||||
return list(zip(np.r_[0, br].tolist(), np.r_[br, len(path)].tolist(), strict=True))
|
||||
|
||||
|
||||
def _mask_events(events: list[Any], mask: dict[str, list[tuple[int, int]]], counts: dict[str, int]) -> list[Any]:
|
||||
"""CYTO-017: trim to unmasked sequence; drop events ≥ 50 % masked."""
|
||||
out = []
|
||||
for ev in events:
|
||||
covered = 0
|
||||
for a, b in mask.get(ev.chrom, []):
|
||||
lo, hi = max(a, ev.start), min(b, ev.end)
|
||||
if lo <= hi:
|
||||
covered += hi - lo + 1
|
||||
if a <= ev.start <= b:
|
||||
ev.start = b + 1
|
||||
ev.flags.append("centromere_trimmed")
|
||||
if a <= ev.end <= b:
|
||||
ev.end = a - 1
|
||||
ev.flags.append("centromere_trimmed")
|
||||
if covered >= 0.5 * (ev.end - ev.start + 1 + covered) or ev.end < ev.start:
|
||||
counts["masked"] = counts.get("masked", 0) + 1
|
||||
continue
|
||||
out.append(ev)
|
||||
return out
|
||||
|
||||
|
||||
# ------------------------------------------------------------------------------------------------ VCF (CYTO-028..031)
|
||||
def _gt(ev: Any, haploid: bool, mode: str) -> str:
|
||||
if mode == "compat":
|
||||
return "1/1"
|
||||
if haploid:
|
||||
return "1"
|
||||
if ev.mosaic:
|
||||
return "0/1"
|
||||
if ev.call in ("LOH", "GAINLOH"):
|
||||
return "1/1"
|
||||
if ev.call == "DEL":
|
||||
return "1/1" if ev.cn == 0 else "0/1"
|
||||
return "./1" if ev.cn >= 4 else "0/1"
|
||||
|
||||
|
||||
def _fmt(v: float | None, nd: int) -> str:
|
||||
return "." if v is None or not math.isfinite(v) else f"{v:.{nd}f}"
|
||||
|
||||
|
||||
def vcf_text(
|
||||
sample_name: str,
|
||||
events: list[Any],
|
||||
haploid: dict[int, bool],
|
||||
contigs: list[tuple[str, int]],
|
||||
header: dict[str, str],
|
||||
opts: CnCallOptions,
|
||||
ref_base: Any,
|
||||
) -> tuple[str, list[tuple[str, int, int, str]]]:
|
||||
"""Header and records; also (chrom, beg0, end, line) for indexing."""
|
||||
th = opts.thresholds
|
||||
extra = {
|
||||
"smoothing": opts.smoothing,
|
||||
"max_mosaic_fraction": th.max_mosaic_fraction,
|
||||
"min_mosaic_fraction": f"{th.min_mosaic_fraction:.2f}",
|
||||
"filter_qual": th.filter_qual,
|
||||
"size_basis": "effective",
|
||||
"cn_cap": opts.cn_cap,
|
||||
"gt_mode": opts.gt_mode,
|
||||
"svlen_mode": opts.svlen_mode,
|
||||
}
|
||||
h = [
|
||||
"##fileformat=VCFv4.3",
|
||||
f"##source=SimurgArray {__version__} cn",
|
||||
f"##reference=file://{header['reference']}",
|
||||
f"##simurg_genome_build={header['build']}",
|
||||
f"##simurg_manifest={header['manifest']}",
|
||||
f"##simurg_cyto_model={header['model']}",
|
||||
"##simurg_thresholds=" + th.as_header(extra),
|
||||
f"##simurg_sex_used={header['sex']};source={header['sex_source']}",
|
||||
f"##OverallPloidy={float(header['ploidy']):.3f}",
|
||||
f"##GCCorrect={'true' if header['gc'] == 'True' else 'false'}",
|
||||
*(f"##contig=<ID={c},length={n}>" for c, n in contigs),
|
||||
'##ALT=<ID=DEL,Description="Copy number loss">',
|
||||
'##ALT=<ID=DUP,Description="Copy number gain">',
|
||||
'##ALT=<ID=LOH,Description="Absence/loss of heterozygosity (AOH/LOH/ROH)">',
|
||||
'##FILTER=<ID=PASS,Description="All filters passed">',
|
||||
f'##FILTER=<ID=Q{th.filter_qual},Description="Quality below {th.filter_qual}">',
|
||||
'##FILTER=<ID=SampleQuality,Description="Sample flagged as potentially low quality: LRR SD above 0.2">',
|
||||
f'##FILTER=<ID=LowFraction,Description="Mosaic event with estimated mosaic fraction below '
|
||||
f'{th.min_mosaic_fraction:.2f}">',
|
||||
'##INFO=<ID=END,Number=1,Type=Integer,Description="Last base of the event (1-based, inclusive)">',
|
||||
'##INFO=<ID=SVTYPE,Number=1,Type=String,Description="DEL, DUP or LOH">',
|
||||
'##INFO=<ID=SVLEN,Number=1,Type=Integer,Description="Event length in bp; negative for deletions">',
|
||||
'##INFO=<ID=EFFLEN,Number=1,Type=Integer,Description="Effective length in bp (gaps above 150x median probe '
|
||||
'spacing replaced by the median spacing)">',
|
||||
'##INFO=<ID=PROBE,Number=1,Type=Integer,Description="Number of probes in the event">',
|
||||
'##INFO=<ID=MOSAIC,Number=0,Type=Flag,Description="Mosaic event">',
|
||||
'##INFO=<ID=HIGHFRACTION,Number=0,Type=Flag,Description="Mosaic event promoted to constitutional because its '
|
||||
'mosaic fraction reached the maximum">',
|
||||
'##INFO=<ID=LOHTYPE,Number=1,Type=String,Description="AOH (constitutional, copy-neutral or gain), CNLOH '
|
||||
'(mosaic copy-neutral), GAINLOH (mosaic copy-gain)">',
|
||||
'##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype (see SimurgArray docs, CYTO-031)">',
|
||||
f'##FORMAT=<ID=CN,Number=1,Type=Integer,Description="Copy number; the maximum value {opts.cn_cap} means '
|
||||
f'{opts.cn_cap} or more">',
|
||||
'##FORMAT=<ID=NR,Number=1,Type=Float,Description="Aggregated normalized intensity (median 2^LRR; about CN/2)">',
|
||||
'##FORMAT=<ID=LRD,Number=1,Type=Float,Description="Standard deviation of LRR within the event">',
|
||||
'##FORMAT=<ID=MF,Number=1,Type=Float,Description="Estimated mosaic fraction; missing if not mosaic or not '
|
||||
'estimable">',
|
||||
f"#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\t{sample_name}",
|
||||
]
|
||||
order = {c: i for i, (c, _) in enumerate(contigs)}
|
||||
recs = []
|
||||
for ev in sorted(events, key=lambda e: (order.get(e.chrom, 999), e.start, e.end)):
|
||||
s_, e_ = ev.start, ev.end
|
||||
pos = s_ - 1 if s_ > 1 else 1
|
||||
ref = ref_base(ev.chrom, pos)
|
||||
length = e_ - s_ + 1
|
||||
svlen = -length if ev.svtype == "DEL" and opts.svlen_mode == "signed" else length
|
||||
info = [f"END={e_}", f"SVTYPE={ev.svtype}", f"SVLEN={svlen}", f"EFFLEN={ev.effective_size}",
|
||||
f"PROBE={ev.probes}"] # fmt: skip
|
||||
if ev.mosaic:
|
||||
info.append("MOSAIC")
|
||||
if ev.high_fraction:
|
||||
info.append("HIGHFRACTION")
|
||||
lt = {"LOH": "AOH", "GAINLOH": "AOH", "MLOH": "CNLOH"}.get(ev.call)
|
||||
if ev.call == "GAINLOH" and ev.mosaic:
|
||||
lt = "GAINLOH"
|
||||
if lt:
|
||||
info.append(f"LOHTYPE={lt}")
|
||||
cn = min(ev.cn, opts.cn_cap)
|
||||
mf = ev.mosaic_fraction if ev.mosaic or ev.high_fraction else None
|
||||
line = (f"{ev.chrom}\t{pos}\t{ev.id_prefix}:{ev.chrom}:{s_}-{e_}\t{ref}\t<{ev.svtype}>\t{int(ev.quality)}\t"
|
||||
f"{';'.join(ev.filters)}\t{';'.join(info)}\tGT:CN:NR:LRD:MF\t"
|
||||
f"{_gt(ev, haploid.get(id(ev), False), opts.gt_mode)}:{cn}:{_fmt(ev.nr, 4)}:{_fmt(ev.lrd, 4)}:"
|
||||
f"{_fmt(mf, 3)}") # fmt: skip
|
||||
recs.append((ev.chrom, pos - 1, e_, line))
|
||||
return "\n".join(h) + "\n", recs
|
||||
|
||||
|
||||
def write_vcf(path: Path, header: str, recs: list[tuple[str, int, int, str]], bgzip: bool) -> None:
|
||||
from simurg.core.atomic import atomic_open
|
||||
from simurg.formats.bgzf import BgzfWriter, TabixBuilder, write_tbi
|
||||
|
||||
tbx = TabixBuilder()
|
||||
with atomic_open(path, "wb") as raw:
|
||||
if bgzip:
|
||||
w = BgzfWriter(raw) # type: ignore[arg-type]
|
||||
w.write_record(header.encode())
|
||||
for chrom, beg, end, line in recs:
|
||||
vs, ve = w.write_record((line + "\n").encode())
|
||||
tbx.add(chrom, beg, end, vs, ve)
|
||||
w.close()
|
||||
else:
|
||||
raw.write((header + "".join(r[3] + "\n" for r in recs)).encode()) # type: ignore[arg-type]
|
||||
if bgzip:
|
||||
write_tbi(Path(str(path) + ".tbi"), tbx)
|
||||
|
||||
|
||||
# ------------------------------------------------------------------------------------------------ batch driver
|
||||
_W: dict[str, Any] = {}
|
||||
|
||||
|
||||
def _init_worker(path: str) -> None:
|
||||
import pickle
|
||||
|
||||
with open(path, "rb") as fh:
|
||||
_setup(pickle.load(fh))
|
||||
|
||||
|
||||
def _setup(state: dict[str, Any]) -> None:
|
||||
_W.update(state)
|
||||
fasta = state["opts"].fasta
|
||||
if fasta is not None:
|
||||
genome = Genome(fasta)
|
||||
_W["genome"] = genome
|
||||
_W["cmap"] = ContigMap([e.name for e in genome.entries], "auto")
|
||||
|
||||
|
||||
def _work(si: SampleInput) -> dict[str, Any]:
|
||||
ctx, opts = _W["ctx"], _W["opts"]
|
||||
try:
|
||||
res = call_sample(si, ctx, opts)
|
||||
except Exception as e: # CYTO-040: one sample never stops the batch
|
||||
code = e.code if isinstance(e, SimurgError) else "E_CN_FAILED"
|
||||
return {"row": {"sample_id": si.sample_id, "sample_name": si.sample_name, "status": "FAILED",
|
||||
"qc_flags": f"{code}: {e}"}, "events": [], "error": code} # fmt: skip
|
||||
row = res["row"]
|
||||
genome = _W.get("genome")
|
||||
|
||||
def ref_base(chrom: str, pos: int) -> str:
|
||||
if genome is None:
|
||||
return "N"
|
||||
contig = _W["cmap"].get(chrom)
|
||||
return genome.fetch(contig, pos - 1, pos) or "N" if contig else "N"
|
||||
|
||||
header = {**_W["vcf_header"], "sex": row["sex_used"], "sex_source": row["sex_source"],
|
||||
"ploidy": str(row["overall_ploidy"]), "gc": str(row["gc_correct"])} # fmt: skip
|
||||
htxt, recs = vcf_text(si.sample_name, res["events"], res["haploid"], ctx["contigs"], header, opts, ref_base)
|
||||
out = opts.out / f"{si.sample_id}.cyto.cnv.vcf{'.gz' if opts.bgzip else ''}"
|
||||
write_vcf(out, htxt, recs, opts.bgzip)
|
||||
seg = pa.Table.from_pylist(res["segments"]) if res["segments"] else pa.table({"chrom": pa.array([], pa.string())})
|
||||
pq.write_table(seg, opts.out / f"{si.sample_id}.cn_segments.parquet", compression="zstd")
|
||||
evs = [
|
||||
{"id": f"{ev.id_prefix}:{ev.chrom}:{ev.start}-{ev.end}", "chrom": ev.chrom, "start": ev.start, "end": ev.end,
|
||||
"call_type": ev.svtype, "call": ev.call, "loh_type": ev.loh_type, "mosaic": ev.mosaic,
|
||||
"mosaic_fraction": ev.mosaic_fraction, "high_fraction": ev.high_fraction, "copy_number": ev.cn,
|
||||
"quality": round(ev.quality, 1), "filters": ev.filters, "size": ev.size, "effective_size": ev.effective_size,
|
||||
"probe_count": ev.probes, "percent_het": ev.percent_het, "lrr_median": ev.lrr_median, "lrr_dev": ev.lrd,
|
||||
"baf_dev": ev.baf_dev, "nr": ev.nr, "flags": ev.flags}
|
||||
for ev in res["events"]
|
||||
] # fmt: skip
|
||||
return {"row": row, "events": evs, "counts": res["counts"]}
|
||||
|
||||
|
||||
def _read_summary(path: Path | None) -> dict[str, dict[str, str]]:
|
||||
if path is None or not path.exists():
|
||||
return {}
|
||||
lines = path.read_text(encoding="utf-8").splitlines()
|
||||
head = lines[0].split("\t")
|
||||
return {r[0]: dict(zip(head, r, strict=False)) for r in (ln.split("\t") for ln in lines[1:] if ln)}
|
||||
|
||||
|
||||
def run_cn_call(opts: CnCallOptions, rep: Reporter) -> int:
|
||||
"""``simurg cn call`` (CLI-028)."""
|
||||
import multiprocessing as mp
|
||||
import os
|
||||
import pickle
|
||||
import tempfile
|
||||
from concurrent.futures import ProcessPoolExecutor
|
||||
|
||||
from simurg.core.atomic import atomic_write_text, check_outputs_absent
|
||||
from simurg.core.provenance import RunManifest
|
||||
from simurg.formats.discovery import discover_gtcs
|
||||
from simurg.pipeline.qc import write_tsv
|
||||
|
||||
rm = RunManifest(command="cn call", argv=opts.argv)
|
||||
rep.run_id = rm.run_id
|
||||
if opts.gtcs is None:
|
||||
raise SimurgError("E_NO_INPUT", "give --gtcs (folder of GTC files)")
|
||||
gtcs = discover_gtcs(opts.gtcs)
|
||||
summary_path = opts.summary or (opts.gtcs / "sample_summary.tsv")
|
||||
summary = _read_summary(summary_path)
|
||||
ext = ".cyto.cnv.vcf.gz" if opts.bgzip else ".cyto.cnv.vcf"
|
||||
check_outputs_absent(
|
||||
[opts.out / f"{k}{ext}" for k in gtcs] + [opts.out / "cn_sample_summary.tsv"], force=opts.force
|
||||
)
|
||||
th = opts.thresholds
|
||||
for k in list(th.probes):
|
||||
if th.probes[k] < 10:
|
||||
rep.warn("W_CN_PROBE_FLOOR", f"--min-{k}-probes {th.probes[k]} raised to 10")
|
||||
th.probes[k] = 10
|
||||
man = read_bpm(opts.bpm)
|
||||
model = read_model(opts.model, man)
|
||||
build = str(model.header.get("genome_build") or "GRCh38")
|
||||
chroms = np.array(["X" if c == "XY" else (c[3:] if c.lower().startswith("chr") else c) for c in man.chroms])
|
||||
pos = man.map_infos.astype(np.int64)
|
||||
by_chrom = {}
|
||||
for c in CALL_CHROMS:
|
||||
sel = np.flatnonzero((chroms == c) & (pos > 0))
|
||||
by_chrom[c] = sel[np.argsort(pos[sel], kind="stable")]
|
||||
auto = np.isin(chroms, CALL_CHROMS[:22]) & (pos > 0)
|
||||
genome = None
|
||||
cmap = None
|
||||
lengths: dict[str, int] = {}
|
||||
if opts.fasta is not None:
|
||||
genome = Genome(opts.fasta)
|
||||
cmap = ContigMap([e.name for e in genome.entries], "auto")
|
||||
for c in CALL_CHROMS:
|
||||
contig = cmap.get(c)
|
||||
if contig:
|
||||
lengths[c] = genome.contig_length(contig)
|
||||
else:
|
||||
rep.warn("W_CN_NO_FASTA", "no --fasta: VCF REF padding bases are written as N")
|
||||
contigs = [(c, lengths.get(c, int(pos[by_chrom[c]].max()) if len(by_chrom[c]) else 0)) for c in CALL_CHROMS
|
||||
if len(by_chrom[c])] # fmt: skip
|
||||
ctx = {"chroms": chroms, "pos": pos, "model": model, "auto": auto, "by_chrom": by_chrom, "build": build,
|
||||
"mask": read_mask(opts.mask_bed, build), "gamma": 0.5, "contigs": contigs} # fmt: skip
|
||||
model_sha = hashlib.sha256(opts.model.read_bytes()).hexdigest()
|
||||
man_sha = hashlib.sha256(opts.bpm.read_bytes()).hexdigest()
|
||||
schema = model.header.get("schema_version")
|
||||
vcf_header = {
|
||||
"reference": str(opts.fasta.resolve()) if opts.fasta else "none",
|
||||
"build": build,
|
||||
"manifest": f"{man.manifest_name};sha256={man_sha}",
|
||||
"model": f"{opts.model.name};sha256={model_sha};schema={schema}",
|
||||
}
|
||||
samples = []
|
||||
for k, path in gtcs.items():
|
||||
r = summary.get(k, {})
|
||||
sex = r.get("sex_est", "U") or "U"
|
||||
src = "summary" if sex in ("M", "F") else "none"
|
||||
if sex not in ("M", "F") and r.get("sex_reported") in ("M", "F"):
|
||||
sex, src = r["sex_reported"], "sheet"
|
||||
try:
|
||||
ld = float(r.get("logr_dev", "nan"))
|
||||
except ValueError:
|
||||
ld = math.nan
|
||||
samples.append(SampleInput(k, r.get("sample_name", k) or k, str(path), sex, ld, src))
|
||||
opts.out.mkdir(parents=True, exist_ok=True)
|
||||
state = {"ctx": ctx, "opts": opts, "vcf_header": vcf_header}
|
||||
workers = max(1, min(opts.threads or (os.cpu_count() or 1), len(samples)))
|
||||
rep.info(f"CN calling of {len(samples)} sample(s) with {workers} worker(s)")
|
||||
if workers == 1:
|
||||
_setup(state)
|
||||
try:
|
||||
results = [_work(s) for s in samples]
|
||||
finally:
|
||||
_W.clear()
|
||||
else:
|
||||
fd, sp = tempfile.mkstemp(prefix=".simurg-workers-", suffix=".pkl", dir=opts.out)
|
||||
try:
|
||||
with os.fdopen(fd, "wb") as fh:
|
||||
pickle.dump(state, fh, protocol=5)
|
||||
with ProcessPoolExecutor(workers, mp_context=mp.get_context("spawn"), initializer=_init_worker,
|
||||
initargs=(sp,)) as ex: # fmt: skip
|
||||
results = list(ex.map(_work, samples))
|
||||
finally:
|
||||
Path(sp).unlink(missing_ok=True)
|
||||
cols = ["sample_id", "sample_name", "status", "logr_dev", "overall_ploidy", "gc_correct", "gcwf", "sex_used",
|
||||
"sex_source", "n_del", "n_dup", "n_loh", "n_gainloh", "n_mosaic", "n_filtered", "qc_flags"] # fmt: skip
|
||||
rows = []
|
||||
for s, res in zip(samples, results, strict=True):
|
||||
row = {**res["row"], "sample_folder": str(Path(s.gtc).parent)}
|
||||
if row["status"] != "OK":
|
||||
rep.error(res.get("error", "E_CN_FAILED"), f"{s.sample_id}: {row['qc_flags']}", sample_id=s.sample_id)
|
||||
elif row["sex_source"] == "resolved":
|
||||
rep.info(f"Sample {s.sample_id} sex updated from U to {row['sex_used']}", "I_CN_SEX_RESOLVED")
|
||||
rows.append(row)
|
||||
write_tsv(opts.out / "cn_sample_summary.tsv", cols, rows)
|
||||
th_dict = {f"min_{k}_probes": th.probes[k] for k in th.probes} | {
|
||||
f"min_{k}_size_bp": th.size_bp[k] for k in th.size_bp
|
||||
}
|
||||
sidecar = {
|
||||
"schema": "simurg.cn_call/1", "software": f"SimurgArray {__version__}", "run_id": rm.run_id,
|
||||
"genome_build": build, "model": vcf_header["model"], "manifest": vcf_header["manifest"],
|
||||
"summary": str(summary_path) if summary else None,
|
||||
"thresholds": th_dict | {"max_mosaic_fraction": th.max_mosaic_fraction,
|
||||
"min_mosaic_fraction": th.min_mosaic_fraction, "filter_qual": th.filter_qual,
|
||||
"smoothing": opts.smoothing, "cn_cap": opts.cn_cap, "size_basis": "effective"},
|
||||
"samples": [{**r["row"], "events": r["events"]} for r in results],
|
||||
} # fmt: skip
|
||||
atomic_write_text(opts.out / "cn_call_manifest.json", json.dumps(sidecar, indent=1, default=_json_default))
|
||||
ok = sum(r["row"]["status"] == "OK" for r in results)
|
||||
rm.summary = {"samples": len(results), "ok": ok}
|
||||
rm.status = "completed"
|
||||
rm.write(opts.out / "run_manifest.json")
|
||||
rep.info(f"done: {ok}/{len(results)} samples; {sum(len(r['events']) for r in results)} events")
|
||||
if ok == len(results):
|
||||
return ExitCode.OK
|
||||
return ExitCode.PARTIAL if ok else ExitCode.ALL_FAILED
|
||||
|
||||
|
||||
def _json_default(v: Any) -> Any:
|
||||
if isinstance(v, np.generic):
|
||||
return v.item()
|
||||
if isinstance(v, float) and not math.isfinite(v):
|
||||
return None
|
||||
raise TypeError(type(v))
|
||||
Reference in new issue
Block a user