From 724e87dc66b2199beee6398c057b349c13503d0a Mon Sep 17 00:00:00 2001 From: Kemal Yaylali Date: Tue, 29 Sep 2026 01:13:51 +0100 Subject: [PATCH] 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. --- docs/messages.md | 10 + src/simurg/cli/main.py | 101 +++++ src/simurg/cn/events.py | 245 ++++++++++++ src/simurg/cn/gcwave.py | 38 ++ src/simurg/cn/hmm.py | 190 ++++++++++ src/simurg/cn/model.py | 73 ++++ src/simurg/core/errors.py | 16 + src/simurg/pipeline/cn.py | 779 ++++++++++++++++++++++++++++++++++++++ 8 files changed, 1452 insertions(+) create mode 100644 src/simurg/cn/events.py create mode 100644 src/simurg/cn/gcwave.py create mode 100644 src/simurg/cn/hmm.py create mode 100644 src/simurg/cn/model.py create mode 100644 src/simurg/pipeline/cn.py diff --git a/docs/messages.md b/docs/messages.md index b1d65c4..804104f 100644 --- a/docs/messages.md +++ b/docs/messages.md @@ -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 | diff --git a/src/simurg/cli/main.py b/src/simurg/cli/main.py index 097a840..2813ee2 100644 --- a/src/simurg/cli/main.py +++ b/src/simurg/cli/main.py @@ -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") diff --git a/src/simurg/cn/events.py b/src/simurg/cn/events.py new file mode 100644 index 0000000..3d93454 --- /dev/null +++ b/src/simurg/cn/events.py @@ -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"] diff --git a/src/simurg/cn/gcwave.py b/src/simurg/cn/gcwave.py new file mode 100644 index 0000000..1009c6d --- /dev/null +++ b/src/simurg/cn/gcwave.py @@ -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]))) diff --git a/src/simurg/cn/hmm.py b/src/simurg/cn/hmm.py new file mode 100644 index 0000000..e039744 --- /dev/null +++ b/src/simurg/cn/hmm.py @@ -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 diff --git a/src/simurg/cn/model.py b/src/simurg/cn/model.py new file mode 100644 index 0000000..15f7d39 --- /dev/null +++ b/src/simurg/cn/model.py @@ -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, + ) diff --git a/src/simurg/core/errors.py b/src/simurg/core/errors.py index 43b204a..404c278 100644 --- a/src/simurg/core/errors.py +++ b/src/simurg/core/errors.py @@ -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"), diff --git a/src/simurg/pipeline/cn.py b/src/simurg/pipeline/cn.py new file mode 100644 index 0000000..2ea9a0f --- /dev/null +++ b/src/simurg/pipeline/cn.py @@ -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=" for c, n in contigs), + '##ALT=', + '##ALT=', + '##ALT=', + '##FILTER=', + f'##FILTER=', + '##FILTER=', + f'##FILTER=', + '##INFO=', + '##INFO=', + '##INFO=', + '##INFO=', + '##INFO=', + '##INFO=', + '##INFO=', + '##INFO=', + '##FORMAT=', + f'##FORMAT=', + '##FORMAT=', + '##FORMAT=', + '##FORMAT=', + 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))