v2 CN: cn annotate (JSON report, ISCN 2020, cytobands, genes) and the annotation bundle builder

cn call now records per-chromosome and sample statistics (percent het, LRR stats, BAF
deviations, robust LogR dev, call rate) in the sidecar; cn annotate writes
<sample>.cyto.json (CYTO-033) with ISCN microarray/traditional strings (CYTO-034), bands
and RefSeq gene lists (CYTO-035), LOH summaries, weighted median CN (CYTO-026); annotation
bundle from UCSC cytoBand + ncbiRefSeq (CYTO-036, IN-054). Tests for ISCN strings, gene
overlap, sidecar-missing error and an end-to-end call+annotate.

src/simurg/cli/main.py

@@ -767,6 +767,32 @@ def cn_call_cmd(
     run_guarded(rep, lambda: run_cn_call(opts, rep))
 
 
+@cn_app.command("annotate")
+def cn_annotate_cmd(
+    vcfs: Annotated[Path, typer.Option("--vcfs", help="Output folder of cn call (required).")],
+    bundle: Annotated[
+        Path | None, typer.Option("--annotation-bundle", help="Bundle from scripts/build_annotation_bundle.sh.")
+    ] = None,
+    out: Annotated[Path | None, typer.Option("--out", help="Output folder (default: --vcfs).")] = None,
+    min_qual: Annotated[int, typer.Option("--min-qual", min=0, max=60)] = 20,
+    size_basis: Annotated[str, typer.Option("--size-basis", help="effective|raw")] = "effective",
+    iscn_min_size: Annotated[int, typer.Option("--iscn-min-size", min=0, help="bp")] = 0,
+    force: ForceOpt = False,
+    log_json: LogJson = False,
+    debug: Debug = False,
+    quiet: Quiet = False,
+) -> None:
+    """cn call output -> per-sample JSON report with ISCN 2020 strings, cytobands and genes (CLI-029)."""
+    from simurg.pipeline.cn import CnAnnotateOptions, run_cn_annotate
+
+    if size_basis not in ("effective", "raw"):
+        raise typer.BadParameter("--size-basis effective|raw")
+    rep = _reporter("cn annotate", out or vcfs, log_json, debug, quiet)
+    opts = CnAnnotateOptions(vcfs=vcfs, bundle=bundle, out=out, min_qual=min_qual, size_basis=size_basis,
+                             iscn_min_size=iscn_min_size, force=force)  # fmt: skip
+    run_guarded(rep, lambda: run_cn_annotate(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/annotate.py

@@ -0,0 +1,130 @@
+"""CN event annotation: cytobands, genes, ISCN 2020 strings — spec CYTO-033 … CYTO-035. Pure functions, no I/O.
+
+[DOC: ISCN 2020 (An International System for Human Cytogenomic Nomenclature), microarray nomenclature section;
+UCSC cytoBand tables] Microarray form ``arr[<build>] <chr><band_start>[<band_end>](<S>_<E>)x<CN>`` per event, LOH as
+``x2 hmz`` (gain-LOH ``x3 hmz``), mosaic as ``x<base>~<CN>``; normal results ``arr[<build>](1-22,X)x2`` (female),
+``(1-22)x2,(XY)x1`` (male), ``(1-22)x2`` (sex unresolved). Traditional form ``del(<chr>)(<bands>)`` /
+``dup(<chr>)(<bands>)`` for CNVs only.
+"""
+
+from __future__ import annotations
+
+import bisect
+from collections.abc import Sequence
+from dataclasses import dataclass
+from typing import Any
+
+CHROM_ORDER = {str(i): i for i in range(1, 23)} | {"X": 23, "Y": 24}
+
+
+@dataclass
+class Bands:
+    """1-based inclusive band intervals per chromosome, sorted."""
+
+    starts: dict[str, list[int]]
+    ends: dict[str, list[int]]
+    names: dict[str, list[str]]
+
+    @classmethod
+    def from_ucsc(cls, text: str) -> Bands:
+        rows: dict[str, list[tuple[int, int, str]]] = {}
+        for line in text.splitlines():
+            f = line.split("\t")
+            if len(f) < 4:
+                continue
+            c = f[0][3:] if f[0].lower().startswith("chr") else f[0]
+            if "_" in c:
+                continue
+            rows.setdefault(c, []).append((int(f[1]) + 1, int(f[2]), f[3]))
+        s, e, n = {}, {}, {}
+        for c, v in rows.items():
+            v.sort()
+            s[c], e[c], n[c] = [x[0] for x in v], [x[1] for x in v], [x[2] for x in v]
+        return cls(s, e, n)
+
+    def band(self, chrom: str, pos: int) -> str:
+        st = self.starts.get(chrom)
+        if not st:
+            return ""
+        i = max(0, bisect.bisect_right(st, pos) - 1)
+        return self.names[chrom][i]
+
+
+@dataclass
+class Genes:
+    """Transcripts: per chromosome sorted by start, with the gene label (HGNC symbol or Ensembl/RefSeq ID)."""
+
+    starts: dict[str, list[int]]
+    ends: dict[str, list[int]]
+    labels: dict[str, list[str]]
+    max_len: dict[str, int]
+
+    @classmethod
+    def from_rows(cls, rows: Sequence[tuple[str, int, int, str]]) -> Genes:
+        by: dict[str, list[tuple[int, int, str]]] = {}
+        for c, s, e, lab in rows:
+            by.setdefault(c, []).append((s, e, lab))
+        st, en, lb, ml = {}, {}, {}, {}
+        for c, v in by.items():
+            v.sort()
+            st[c], en[c], lb[c] = [x[0] for x in v], [x[1] for x in v], [x[2] for x in v]
+            ml[c] = max(x[1] - x[0] for x in v)
+        return cls(st, en, lb, ml)
+
+    def overlapping(self, chrom: str, start: int, end: int) -> list[str]:
+        """CYTO-035: every transcript overlapping by ≥ 1 bp; labels deduplicated in transcript-start order."""
+        st = self.starts.get(chrom)
+        if not st:
+            return []
+        lo = bisect.bisect_left(st, start - self.max_len[chrom])
+        hi = bisect.bisect_right(st, end)
+        out: list[str] = []
+        seen: set[str] = set()
+        for i in range(lo, hi):
+            if self.ends[chrom][i] >= start and self.labels[chrom][i] not in seen:
+                seen.add(self.labels[chrom][i])
+                out.append(self.labels[chrom][i])
+        return out
+
+
+def _band_span(bands: Bands, ev: dict[str, Any]) -> tuple[str, str]:
+    return bands.band(ev["chrom"], ev["start"]), bands.band(ev["chrom"], ev["end"])
+
+
+def _chrom_label(c: str) -> str:
+    return c
+
+
+def iscn_microarray_event(ev: dict[str, Any], bands: Bands) -> str:
+    b1, b2 = _band_span(bands, ev)
+    span = b1 if b1 == b2 else f"{b1}{b2}"
+    cn = int(ev["copy_number"])
+    if ev.get("mosaic"):
+        base = int(ev.get("base_cn", 2))
+        copies = f"x{base}~{cn if ev['call'] != 'MLOH' else base}"
+        if ev["call"] == "MLOH":
+            copies += " hmz"
+    elif ev["call_type"] == "LOH":
+        copies = f"x{cn} hmz"
+    else:
+        copies = f"x{cn}"
+    return f"{_chrom_label(ev['chrom'])}{span}({ev['start']}_{ev['end']}){copies}"
+
+
+def iscn_traditional_event(ev: dict[str, Any], bands: Bands) -> str | None:
+    if ev["call_type"] not in ("DEL", "DUP"):
+        return None
+    b1, b2 = _band_span(bands, ev)
+    kind = "del" if ev["call_type"] == "DEL" else "dup"
+    return f"{kind}({ev['chrom']})({b1}{b2})"
+
+
+def iscn_sample(events: list[dict[str, Any]], bands: Bands, build: str, sex: str) -> tuple[str, str]:
+    """(traditional, microarray) strings for the reported events (already filtered)."""
+    evs = sorted(events, key=lambda e: (CHROM_ORDER.get(e["chrom"], 99), e["start"]))
+    if not evs:
+        normal = {"F": "(1-22,X)x2", "M": "(1-22)x2,(XY)x1"}.get(sex, "(1-22)x2")
+        return "", f"arr[{build}]{normal}"
+    micro = f"arr[{build}] " + ",".join(iscn_microarray_event(e, bands) for e in evs)
+    trad = ",".join(t for t in (iscn_traditional_event(e, bands) for e in evs) if t)
+    return trad, micro

src/simurg/pipeline/cn.py

@@ -418,6 +418,37 @@ def call_sample(si: SampleInput, ctx: dict[str, Any], opts: CnCallOptions) -> di
         row["qc_flags"] = ";".join([*flags, "excess_events:" + ",".join(excess)])
     row["lrr_sd_segmentation"] = seg_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}
 
 
@@ -628,7 +659,7 @@ def _work(si: SampleInput) -> dict[str, Any]:
          "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}
+         "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"]}
@@ -696,7 +727,8 @@ def run_cn_call(opts: CnCallOptions, rep: Reporter) -> int:
     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
+           "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")
@@ -755,6 +787,7 @@ def run_cn_call(opts: CnCallOptions, rep: Reporter) -> int:
     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,
@@ -778,3 +811,159 @@ def _json_default(v: Any) -> Any:
     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_<type>_probes / min_<type>_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 ``<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]