"""``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 ref_het_sd: npt.NDArray[np.float32] | None = None # per-probe robust SD of 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.discovery import discover_gtcs from simurg.formats.gtc import read_gtc paths = sorted(discover_gtcs(reference_gtcs).values()) # E_NO_INPUT when there is none 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") from simurg.cn.model import reference_het_spread 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_het_sd": reference_het_spread(baf), "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"), ref_het_sd=col("ref_het_sd")) # 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 mosaic_prior: float = 0.004 # [OURS] tuning knobs (hmm.MOSAIC_PRIOR, hmm.SWITCH_D) switch_d: float = 5e8 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 _baf_noise( baf: npt.NDArray[np.float64], geno: npt.NDArray[np.uint8], auto: npt.NDArray[np.bool_], ref_het_sd: npt.NDArray[np.float32] | None, baf_sd: float, ) -> npt.NDArray[np.float64]: # fmt: skip """CYTO-041: expected heterozygous BAF noise per probe, the reference spread scaled to this sample (median |BAF − 0.5| / spread over autosomal heterozygous calls = 0.6745 × scale); ``baf_sd`` without a reference spread or where a probe has none.""" out = np.full(len(baf), baf_sd, dtype=np.float64) if ref_het_sd is None: return out h = np.asarray(ref_het_sd, dtype=np.float64) hc = auto & (geno == 2) & np.isfinite(baf) & np.isfinite(h) & (h > 0) if hc.sum() < 100: return out scale = float(np.median(np.abs(baf[hc] - 0.5) / np.maximum(h[hc], 1e-3))) / 0.6745 ok = np.isfinite(h) & (h > 0) out[ok] = np.maximum(h[ok], 1e-3) * scale return out 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])) # CYTO-014 gates on the noise of the LRR the caller actually uses (after the batch reference and GC correction): # the summary's logr_dev also measures cluster-file mismatch, which the reference removes (K-027) gate_sd = float(np.nanstd(lrr[auto])) if model.ref_lrr is not None else logr_dev noisy = gate_sd > 0.2 whole_chrom_only = gate_sd > 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) sigma = _baf_noise(baf, geno, auto, model.ref_het_sd, baf_sd) # CYTO-041 dip, hap = hmm.diploid_states(opts.mosaic_prior), hmm.haploid_states(opts.mosaic_prior) 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, switch_d=opts.switch_d) het_sd = model.ref_het_sd[idx].astype(np.float64) if model.ref_het_sd is not None else None e = hmm.emissions(states, x_seg, np.full(len(b), np.nan) if haploid else b, pfb[idx], prm, het_sd) sw = hmm.switch_probs(p, opts.switch_d) 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, sigma[idx], pfb[idx]) evs = post_filter(evs, [x_raw[ev.i0 : ev.i1] for ev in evs], th, counts) for s_, e_ in _runs(path): st = states[int(path[s_])] 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["lrr_sd_gate"] = gate_sd row["baf_sd"] = baf_sd # CYTO-033 sample and chromosome statistics (unsmoothed LRR), for cn annotate la = lrr[auto & np.isfinite(lrr)] row["logr_dev_robust"] = _robust_sd(la, 0.0) row["call_rate"] = float(g.call_rate) outside = np.ones(len(lrr), dtype=np.bool_) for ev in kept: sel = ctx["by_chrom"].get(ev.chrom, np.zeros(0, np.int64)) p_ = pos_all[sel] outside[sel[(p_ >= ev.start) & (p_ <= ev.end)]] = False bdev = {} for key, code in (("aa", 1), ("ab", 2), ("bb", 3)): v = baf[auto & outside & (geno == code) & np.isfinite(baf)] bdev[key] = float(np.std(v)) if len(v) > 1 else None row["baf_dev"] = bdev chrom_stats = [] for chrom in CALL_CHROMS: sel = ctx["by_chrom"].get(chrom) if sel is None or len(sel) == 0: continue x = lrr[sel][np.isfinite(lrr[sel])] gt = geno[sel] called = gt > 0 chrom_stats.append({ "id": chrom, "size": ctx["contig_len"].get(chrom, int(pos_all[sel].max())), "effective_size": effective_size(pos_all[sel], median_spacing(pos_all[sel])), "percent_het": float(100 * (gt == 2).sum() / called.sum()) if called.any() else None, "lrr_median": float(np.median(x)) if len(x) else None, "lrr_mean": float(np.mean(x)) if len(x) else None, "lrr_dev": float(np.std(x)) if len(x) > 1 else None, "min_lrr": float(x.min()) if len(x) else None, "max_lrr": float(x.max()) if len(x) else None, }) # fmt: skip row["chromosomes"] = chrom_stats 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, "base_cn": 1 if res["haploid"].get(id(ev)) else 2} 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, "contig_len": dict(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"], "reference": vcf_header["reference"], "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)) # ================================================================================================ cn annotate @dataclass class CnAnnotateOptions: vcfs: Path bundle: Path | None out: Path | None = None min_qual: int = 20 size_basis: str = "effective" iscn_min_size: int = 0 tighten: dict[str, int] = field(default_factory=dict) # min__probes / min__size_bp overrides force: bool = False def _load_bundle(path: Path | None, build: str) -> tuple[Any, Any, dict[str, Any] | None]: import gzip from simurg.cn.annotate import Bands, Genes if path is None: text = gzip.decompress(_builtin_file(f"cytoBand_{build}.txt.gz")).decode() return Bands.from_ucsc(text), None, None with zipfile.ZipFile(path) as z: man = json.loads(z.read("manifest.json")) bands = Bands.from_ucsc(gzip.decompress(z.read("cytoBand.txt.gz")).decode()) rows = [] for line in gzip.decompress(z.read("genes.tsv.gz")).decode().splitlines()[1:]: c, s_, e_, gene, tx = line.split("\t") rows.append((c, int(s_), int(e_), gene or tx)) man["file"] = path.name man["sha256"] = hashlib.sha256(path.read_bytes()).hexdigest() return bands, Genes.from_rows(rows), man def run_cn_annotate(opts: CnAnnotateOptions, rep: Reporter) -> int: """``simurg cn annotate`` (CLI-029): per-sample ``.cyto.json`` with ISCN, bands and genes.""" from simurg.cn.annotate import CHROM_ORDER, iscn_microarray_event, iscn_sample, iscn_traditional_event from simurg.core.atomic import atomic_write_text, check_outputs_absent side_p = opts.vcfs / "cn_call_manifest.json" if not side_p.exists(): raise SimurgError("E_CN_SIDECAR_MISSING", f"{side_p} not found; run simurg cn call first") side = json.loads(side_p.read_text()) out = opts.out or opts.vcfs samples = [s for s in side["samples"] if s.get("status") == "OK"] for s in samples: v = opts.vcfs / f"{s['sample_id']}.cyto.cnv.vcf.gz" if not (v.exists() and Path(str(v) + ".tbi").exists()): raise SimurgError("E_CN_VCF_UNINDEXED", f"{v} (+ .tbi) not found; cn call must run without --no-bgzip") check_outputs_absent([out / f"{s['sample_id']}.cyto.json" for s in samples], force=opts.force) th = dict(side["thresholds"]) for k, val in opts.tighten.items(): # CLI-029: can only tighten if val is not None and val > th.get(k, 0): th[k] = val build = side.get("genome_build", "GRCh38") bands, genes, bundle = _load_bundle(opts.bundle, build) if bundle is None: rep.warn("W_CN_NO_FASTA", "no --annotation-bundle: built-in cytobands only, no gene lists") key = {"DEL": "del", "DUP": "dup", "LOH": "loh", "GAINLOH": "gainloh", "MDEL": "mdel", "MDUP": "mdup", "MLOH": "mloh"} # fmt: skip out.mkdir(parents=True, exist_ok=True) for s in samples: events = [] for ev in s["events"]: k = key[ev["call"]] size = ev["effective_size"] if opts.size_basis == "effective" else ev["size"] if ev["probe_count"] < th[f"min_{k}_probes"] or size < th[f"min_{k}_size_bp"]: continue if ev["quality"] < opts.min_qual: continue e = dict(ev) e["start_cytoband"], e["end_cytoband"] = ( bands.band(ev["chrom"], ev["start"]), bands.band(ev["chrom"], ev["end"]), ) e["iscn_microarray"] = f"arr[{build}] " + iscn_microarray_event(ev, bands) e["iscn_traditional"] = iscn_traditional_event(ev, bands) gl = genes.overlapping(ev["chrom"], ev["start"], ev["end"]) if genes else [] e["genes"], e["gene_count"] = gl, len(gl) e["call_type"], e["loh_type"] = ev["call_type"], ev["loh_type"] events.append(e) reported = [e for e in events if e["filters"] == ["PASS"] and e["size"] >= opts.iscn_min_size] trad, micro = iscn_sample(reported, bands, build, s["sex_used"]) eff_genome = sum(c["effective_size"] for c in s.get("chromosomes", [])) loh = [e for e in events if e["call_type"] == "LOH"] def loh_stats(evs: list[dict[str, Any]], denom: int) -> dict[str, Any]: big = [e for e in evs if e["size"] > 1_000_000] tot = sum(e["effective_size"] for e in big) pct = round(100 * tot / denom, 3) if denom else None return {"num_loh_over_1mb": len(big), "num_loh_over_8mb": sum(e["size"] > 8_000_000 for e in evs), "total_size_loh_over_1mb": tot, "percent_loh": pct} # fmt: skip chroms = [] for c in s.get("chromosomes", []): ce = [e for e in events if e["chrom"] == c["id"]] base = 1 if (c["id"] == "Y" or (c["id"] == "X" and s["sex_used"] == "M")) else 2 cov = sum(e["effective_size"] for e in ce if e["call_type"] != "LOH") mean_cn = (base * max(c["effective_size"] - cov, 0) + sum(e["copy_number"] * e["effective_size"] for e in ce if e["call_type"] != "LOH")) / max(c["effective_size"], 1) # fmt: skip mfs = [e["mosaic_fraction"] for e in ce if e.get("mosaic") and e.get("mosaic_fraction") is not None] chroms.append({**c, "has_mosaicism": any(e.get("mosaic") for e in ce), **loh_stats([e for e in loh if e["chrom"] == c["id"]], c["effective_size"]), "copy_number_median": float(base if cov < c["effective_size"] / 2 else max((e["copy_number"] for e in ce), default=base)), "copy_number_mean": round(mean_cn, 4), "median_mosaic_fraction": float(np.median(mfs)) if mfs else None, "number_del": sum(e["call_type"] == "DEL" and not e["mosaic"] for e in ce), "number_dup": sum(e["call_type"] == "DUP" and not e["mosaic"] for e in ce), "number_loh": sum(e["call_type"] == "LOH" for e in ce), "number_gainloh": sum(e["call"] == "GAINLOH" for e in ce), "number_mosaic": sum(bool(e["mosaic"]) for e in ce)}) # fmt: skip flags = [f for f in str(s.get("qc_flags", "")).split(";") if f] doc = { "schema": "simurg.cyto/1", "software_version": __version__, "annotation_bundle": {k: bundle[k] for k in ("file", "sha256", "sources")} if bundle else None, "reference_genome": side.get("reference"), "genome_build": build, "annotation_type": "constitutional", "iscn_version": "ISCN 2020", "sample_id": s["sample_id"], "sample_name": s.get("sample_name"), "gc_correct": bool(s.get("gc_correct")), "gcwf": s.get("gcwf"), "parameters": {**th, "min_qual": opts.min_qual, "size_basis": opts.size_basis, "iscn_min_size": opts.iscn_min_size}, "overall_ploidy": round(float(s["overall_ploidy"]), 3), "call_rate": s.get("call_rate"), "logr_dev": s.get("logr_dev"), "logr_dev_robust": s.get("logr_dev_robust"), "baf_dev": s.get("baf_dev"), **loh_stats(loh, eff_genome), "copy_number_median": _weighted_median_cn(chroms, events, s["sex_used"]), "sex_estimate": s["sex_used"], "sex_source": s["sex_source"], "qc_flags": flags, "iscn_traditional": trad, "iscn_microarray": micro, "chromosomes": chroms, "events": sorted(events, key=lambda e: (CHROM_ORDER.get(e["chrom"], 99), e["start"])), } # fmt: skip atomic_write_text(out / f"{s['sample_id']}.cyto.json", json.dumps(doc, indent=1, default=_json_default)) rep.info(f"annotated {len(samples)} sample(s) -> {out}") return ExitCode.OK def _weighted_median_cn(chroms: list[dict[str, Any]], events: list[dict[str, Any]], sex: str) -> float: """CYTO-026: length-weighted median CN over all segments (baseline stretches and CNV events).""" parts: list[tuple[float, int]] = [] for c in chroms: base = 1 if (c["id"] == "Y" or (c["id"] == "X" and sex == "M")) else 2 ce = [e for e in events if e["chrom"] == c["id"] and e["call_type"] != "LOH"] cov = sum(e["effective_size"] for e in ce) parts.append((float(base), max(c["effective_size"] - cov, 0))) parts += [(float(e["copy_number"]), e["effective_size"]) for e in ce] if not parts: return 2.0 parts.sort() half = sum(w for _, w in parts) / 2 acc = 0 for cn, w in parts: acc += w if acc >= half: return cn return parts[-1][0]