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

@@ -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 |

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")

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"]

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])))

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

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,
+    )

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"),

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=<ID={c},length={n}>" for c, n in contigs),
+        '##ALT=<ID=DEL,Description="Copy number loss">',
+        '##ALT=<ID=DUP,Description="Copy number gain">',
+        '##ALT=<ID=LOH,Description="Absence/loss of heterozygosity (AOH/LOH/ROH)">',
+        '##FILTER=<ID=PASS,Description="All filters passed">',
+        f'##FILTER=<ID=Q{th.filter_qual},Description="Quality below {th.filter_qual}">',
+        '##FILTER=<ID=SampleQuality,Description="Sample flagged as potentially low quality: LRR SD above 0.2">',
+        f'##FILTER=<ID=LowFraction,Description="Mosaic event with estimated mosaic fraction below '
+        f'{th.min_mosaic_fraction:.2f}">',
+        '##INFO=<ID=END,Number=1,Type=Integer,Description="Last base of the event (1-based, inclusive)">',
+        '##INFO=<ID=SVTYPE,Number=1,Type=String,Description="DEL, DUP or LOH">',
+        '##INFO=<ID=SVLEN,Number=1,Type=Integer,Description="Event length in bp; negative for deletions">',
+        '##INFO=<ID=EFFLEN,Number=1,Type=Integer,Description="Effective length in bp (gaps above 150x median probe '
+        'spacing replaced by the median spacing)">',
+        '##INFO=<ID=PROBE,Number=1,Type=Integer,Description="Number of probes in the event">',
+        '##INFO=<ID=MOSAIC,Number=0,Type=Flag,Description="Mosaic event">',
+        '##INFO=<ID=HIGHFRACTION,Number=0,Type=Flag,Description="Mosaic event promoted to constitutional because its '
+        'mosaic fraction reached the maximum">',
+        '##INFO=<ID=LOHTYPE,Number=1,Type=String,Description="AOH (constitutional, copy-neutral or gain), CNLOH '
+        '(mosaic copy-neutral), GAINLOH (mosaic copy-gain)">',
+        '##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype (see SimurgArray docs, CYTO-031)">',
+        f'##FORMAT=<ID=CN,Number=1,Type=Integer,Description="Copy number; the maximum value {opts.cn_cap} means '
+        f'{opts.cn_cap} or more">',
+        '##FORMAT=<ID=NR,Number=1,Type=Float,Description="Aggregated normalized intensity (median 2^LRR; about CN/2)">',
+        '##FORMAT=<ID=LRD,Number=1,Type=Float,Description="Standard deviation of LRR within the event">',
+        '##FORMAT=<ID=MF,Number=1,Type=Float,Description="Estimated mosaic fraction; missing if not mosaic or not '
+        'estimable">',
+        f"#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\t{sample_name}",
+    ]
+    order = {c: i for i, (c, _) in enumerate(contigs)}
+    recs = []
+    for ev in sorted(events, key=lambda e: (order.get(e.chrom, 999), e.start, e.end)):
+        s_, e_ = ev.start, ev.end
+        pos = s_ - 1 if s_ > 1 else 1
+        ref = ref_base(ev.chrom, pos)
+        length = e_ - s_ + 1
+        svlen = -length if ev.svtype == "DEL" and opts.svlen_mode == "signed" else length
+        info = [f"END={e_}", f"SVTYPE={ev.svtype}", f"SVLEN={svlen}", f"EFFLEN={ev.effective_size}",
+                f"PROBE={ev.probes}"]  # fmt: skip
+        if ev.mosaic:
+            info.append("MOSAIC")
+        if ev.high_fraction:
+            info.append("HIGHFRACTION")
+        lt = {"LOH": "AOH", "GAINLOH": "AOH", "MLOH": "CNLOH"}.get(ev.call)
+        if ev.call == "GAINLOH" and ev.mosaic:
+            lt = "GAINLOH"
+        if lt:
+            info.append(f"LOHTYPE={lt}")
+        cn = min(ev.cn, opts.cn_cap)
+        mf = ev.mosaic_fraction if ev.mosaic or ev.high_fraction else None
+        line = (f"{ev.chrom}\t{pos}\t{ev.id_prefix}:{ev.chrom}:{s_}-{e_}\t{ref}\t<{ev.svtype}>\t{int(ev.quality)}\t"
+                f"{';'.join(ev.filters)}\t{';'.join(info)}\tGT:CN:NR:LRD:MF\t"
+                f"{_gt(ev, haploid.get(id(ev), False), opts.gt_mode)}:{cn}:{_fmt(ev.nr, 4)}:{_fmt(ev.lrd, 4)}:"
+                f"{_fmt(mf, 3)}")  # fmt: skip
+        recs.append((ev.chrom, pos - 1, e_, line))
+    return "\n".join(h) + "\n", recs
+
+
+def write_vcf(path: Path, header: str, recs: list[tuple[str, int, int, str]], bgzip: bool) -> None:
+    from simurg.core.atomic import atomic_open
+    from simurg.formats.bgzf import BgzfWriter, TabixBuilder, write_tbi
+
+    tbx = TabixBuilder()
+    with atomic_open(path, "wb") as raw:
+        if bgzip:
+            w = BgzfWriter(raw)  # type: ignore[arg-type]
+            w.write_record(header.encode())
+            for chrom, beg, end, line in recs:
+                vs, ve = w.write_record((line + "\n").encode())
+                tbx.add(chrom, beg, end, vs, ve)
+            w.close()
+        else:
+            raw.write((header + "".join(r[3] + "\n" for r in recs)).encode())  # type: ignore[arg-type]
+    if bgzip:
+        write_tbi(Path(str(path) + ".tbi"), tbx)
+
+
+# ------------------------------------------------------------------------------------------------ batch driver
+_W: dict[str, Any] = {}
+
+
+def _init_worker(path: str) -> None:
+    import pickle
+
+    with open(path, "rb") as fh:
+        _setup(pickle.load(fh))
+
+
+def _setup(state: dict[str, Any]) -> None:
+    _W.update(state)
+    fasta = state["opts"].fasta
+    if fasta is not None:
+        genome = Genome(fasta)
+        _W["genome"] = genome
+        _W["cmap"] = ContigMap([e.name for e in genome.entries], "auto")
+
+
+def _work(si: SampleInput) -> dict[str, Any]:
+    ctx, opts = _W["ctx"], _W["opts"]
+    try:
+        res = call_sample(si, ctx, opts)
+    except Exception as e:  # CYTO-040: one sample never stops the batch
+        code = e.code if isinstance(e, SimurgError) else "E_CN_FAILED"
+        return {"row": {"sample_id": si.sample_id, "sample_name": si.sample_name, "status": "FAILED",
+                        "qc_flags": f"{code}: {e}"}, "events": [], "error": code}  # fmt: skip
+    row = res["row"]
+    genome = _W.get("genome")
+
+    def ref_base(chrom: str, pos: int) -> str:
+        if genome is None:
+            return "N"
+        contig = _W["cmap"].get(chrom)
+        return genome.fetch(contig, pos - 1, pos) or "N" if contig else "N"
+
+    header = {**_W["vcf_header"], "sex": row["sex_used"], "sex_source": row["sex_source"],
+              "ploidy": str(row["overall_ploidy"]), "gc": str(row["gc_correct"])}  # fmt: skip
+    htxt, recs = vcf_text(si.sample_name, res["events"], res["haploid"], ctx["contigs"], header, opts, ref_base)
+    out = opts.out / f"{si.sample_id}.cyto.cnv.vcf{'.gz' if opts.bgzip else ''}"
+    write_vcf(out, htxt, recs, opts.bgzip)
+    seg = pa.Table.from_pylist(res["segments"]) if res["segments"] else pa.table({"chrom": pa.array([], pa.string())})
+    pq.write_table(seg, opts.out / f"{si.sample_id}.cn_segments.parquet", compression="zstd")
+    evs = [
+        {"id": f"{ev.id_prefix}:{ev.chrom}:{ev.start}-{ev.end}", "chrom": ev.chrom, "start": ev.start, "end": ev.end,
+         "call_type": ev.svtype, "call": ev.call, "loh_type": ev.loh_type, "mosaic": ev.mosaic,
+         "mosaic_fraction": ev.mosaic_fraction, "high_fraction": ev.high_fraction, "copy_number": ev.cn,
+         "quality": round(ev.quality, 1), "filters": ev.filters, "size": ev.size, "effective_size": ev.effective_size,
+         "probe_count": ev.probes, "percent_het": ev.percent_het, "lrr_median": ev.lrr_median, "lrr_dev": ev.lrd,
+         "baf_dev": ev.baf_dev, "nr": ev.nr, "flags": ev.flags}
+        for ev in res["events"]
+    ]  # fmt: skip
+    return {"row": row, "events": evs, "counts": res["counts"]}
+
+
+def _read_summary(path: Path | None) -> dict[str, dict[str, str]]:
+    if path is None or not path.exists():
+        return {}
+    lines = path.read_text(encoding="utf-8").splitlines()
+    head = lines[0].split("\t")
+    return {r[0]: dict(zip(head, r, strict=False)) for r in (ln.split("\t") for ln in lines[1:] if ln)}
+
+
+def run_cn_call(opts: CnCallOptions, rep: Reporter) -> int:
+    """``simurg cn call`` (CLI-028)."""
+    import multiprocessing as mp
+    import os
+    import pickle
+    import tempfile
+    from concurrent.futures import ProcessPoolExecutor
+
+    from simurg.core.atomic import atomic_write_text, check_outputs_absent
+    from simurg.core.provenance import RunManifest
+    from simurg.formats.discovery import discover_gtcs
+    from simurg.pipeline.qc import write_tsv
+
+    rm = RunManifest(command="cn call", argv=opts.argv)
+    rep.run_id = rm.run_id
+    if opts.gtcs is None:
+        raise SimurgError("E_NO_INPUT", "give --gtcs (folder of GTC files)")
+    gtcs = discover_gtcs(opts.gtcs)
+    summary_path = opts.summary or (opts.gtcs / "sample_summary.tsv")
+    summary = _read_summary(summary_path)
+    ext = ".cyto.cnv.vcf.gz" if opts.bgzip else ".cyto.cnv.vcf"
+    check_outputs_absent(
+        [opts.out / f"{k}{ext}" for k in gtcs] + [opts.out / "cn_sample_summary.tsv"], force=opts.force
+    )
+    th = opts.thresholds
+    for k in list(th.probes):
+        if th.probes[k] < 10:
+            rep.warn("W_CN_PROBE_FLOOR", f"--min-{k}-probes {th.probes[k]} raised to 10")
+            th.probes[k] = 10
+    man = read_bpm(opts.bpm)
+    model = read_model(opts.model, man)
+    build = str(model.header.get("genome_build") or "GRCh38")
+    chroms = np.array(["X" if c == "XY" else (c[3:] if c.lower().startswith("chr") else c) for c in man.chroms])
+    pos = man.map_infos.astype(np.int64)
+    by_chrom = {}
+    for c in CALL_CHROMS:
+        sel = np.flatnonzero((chroms == c) & (pos > 0))
+        by_chrom[c] = sel[np.argsort(pos[sel], kind="stable")]
+    auto = np.isin(chroms, CALL_CHROMS[:22]) & (pos > 0)
+    genome = None
+    cmap = None
+    lengths: dict[str, int] = {}
+    if opts.fasta is not None:
+        genome = Genome(opts.fasta)
+        cmap = ContigMap([e.name for e in genome.entries], "auto")
+        for c in CALL_CHROMS:
+            contig = cmap.get(c)
+            if contig:
+                lengths[c] = genome.contig_length(contig)
+    else:
+        rep.warn("W_CN_NO_FASTA", "no --fasta: VCF REF padding bases are written as N")
+    contigs = [(c, lengths.get(c, int(pos[by_chrom[c]].max()) if len(by_chrom[c]) else 0)) for c in CALL_CHROMS
+               if len(by_chrom[c])]  # fmt: skip
+    ctx = {"chroms": chroms, "pos": pos, "model": model, "auto": auto, "by_chrom": by_chrom, "build": build,
+           "mask": read_mask(opts.mask_bed, build), "gamma": 0.5, "contigs": contigs}  # fmt: skip
+    model_sha = hashlib.sha256(opts.model.read_bytes()).hexdigest()
+    man_sha = hashlib.sha256(opts.bpm.read_bytes()).hexdigest()
+    schema = model.header.get("schema_version")
+    vcf_header = {
+        "reference": str(opts.fasta.resolve()) if opts.fasta else "none",
+        "build": build,
+        "manifest": f"{man.manifest_name};sha256={man_sha}",
+        "model": f"{opts.model.name};sha256={model_sha};schema={schema}",
+    }
+    samples = []
+    for k, path in gtcs.items():
+        r = summary.get(k, {})
+        sex = r.get("sex_est", "U") or "U"
+        src = "summary" if sex in ("M", "F") else "none"
+        if sex not in ("M", "F") and r.get("sex_reported") in ("M", "F"):
+            sex, src = r["sex_reported"], "sheet"
+        try:
+            ld = float(r.get("logr_dev", "nan"))
+        except ValueError:
+            ld = math.nan
+        samples.append(SampleInput(k, r.get("sample_name", k) or k, str(path), sex, ld, src))
+    opts.out.mkdir(parents=True, exist_ok=True)
+    state = {"ctx": ctx, "opts": opts, "vcf_header": vcf_header}
+    workers = max(1, min(opts.threads or (os.cpu_count() or 1), len(samples)))
+    rep.info(f"CN calling of {len(samples)} sample(s) with {workers} worker(s)")
+    if workers == 1:
+        _setup(state)
+        try:
+            results = [_work(s) for s in samples]
+        finally:
+            _W.clear()
+    else:
+        fd, sp = tempfile.mkstemp(prefix=".simurg-workers-", suffix=".pkl", dir=opts.out)
+        try:
+            with os.fdopen(fd, "wb") as fh:
+                pickle.dump(state, fh, protocol=5)
+            with ProcessPoolExecutor(workers, mp_context=mp.get_context("spawn"), initializer=_init_worker,
+                                     initargs=(sp,)) as ex:  # fmt: skip
+                results = list(ex.map(_work, samples))
+        finally:
+            Path(sp).unlink(missing_ok=True)
+    cols = ["sample_id", "sample_name", "status", "logr_dev", "overall_ploidy", "gc_correct", "gcwf", "sex_used",
+            "sex_source", "n_del", "n_dup", "n_loh", "n_gainloh", "n_mosaic", "n_filtered", "qc_flags"]  # fmt: skip
+    rows = []
+    for s, res in zip(samples, results, strict=True):
+        row = {**res["row"], "sample_folder": str(Path(s.gtc).parent)}
+        if row["status"] != "OK":
+            rep.error(res.get("error", "E_CN_FAILED"), f"{s.sample_id}: {row['qc_flags']}", sample_id=s.sample_id)
+        elif row["sex_source"] == "resolved":
+            rep.info(f"Sample {s.sample_id} sex updated from U to {row['sex_used']}", "I_CN_SEX_RESOLVED")
+        rows.append(row)
+    write_tsv(opts.out / "cn_sample_summary.tsv", cols, rows)
+    th_dict = {f"min_{k}_probes": th.probes[k] for k in th.probes} | {
+        f"min_{k}_size_bp": th.size_bp[k] for k in th.size_bp
+    }
+    sidecar = {
+        "schema": "simurg.cn_call/1", "software": f"SimurgArray {__version__}", "run_id": rm.run_id,
+        "genome_build": build, "model": vcf_header["model"], "manifest": vcf_header["manifest"],
+        "summary": str(summary_path) if summary else None,
+        "thresholds": th_dict | {"max_mosaic_fraction": th.max_mosaic_fraction,
+                                 "min_mosaic_fraction": th.min_mosaic_fraction, "filter_qual": th.filter_qual,
+                                 "smoothing": opts.smoothing, "cn_cap": opts.cn_cap, "size_basis": "effective"},
+        "samples": [{**r["row"], "events": r["events"]} for r in results],
+    }  # fmt: skip
+    atomic_write_text(opts.out / "cn_call_manifest.json", json.dumps(sidecar, indent=1, default=_json_default))
+    ok = sum(r["row"]["status"] == "OK" for r in results)
+    rm.summary = {"samples": len(results), "ok": ok}
+    rm.status = "completed"
+    rm.write(opts.out / "run_manifest.json")
+    rep.info(f"done: {ok}/{len(results)} samples; {sum(len(r['events']) for r in results)} events")
+    if ok == len(results):
+        return ExitCode.OK
+    return ExitCode.PARTIAL if ok else ExitCode.ALL_FAILED
+
+
+def _json_default(v: Any) -> Any:
+    if isinstance(v, np.generic):
+        return v.item()
+    if isinstance(v, float) and not math.isfinite(v):
+        return None
+    raise TypeError(type(v))