diff --git a/THIRD_PARTY_NOTICES.md b/THIRD_PARTY_NOTICES.md index 86aea9c..89fe351 100644 --- a/THIRD_PARTY_NOTICES.md +++ b/THIRD_PARTY_NOTICES.md @@ -13,6 +13,7 @@ used only as external subprocess oracles and are not copied or linked. | [Illumina GTCtoVCF](https://github.com/Illumina/GTCtoVCF) — Copyright Illumina, Inc. | Apache-2.0 | Indel source-sequence handling ported to `simurg/harmonize/sites.py` (design-strand orientation, left shift, IUPAC flank matching, deletion/insertion context score); multi-probe consistency idea in `simurg/harmonize/combine.py`; test oracle (run under Python 3 with two local compatibility patches, never shipped) | | [Picard `picard.arrays`](https://github.com/broadinstitute/picard) | MIT | Second format reference | | [PyO3](https://github.com/PyO3/pyo3), [rust-numpy](https://github.com/PyO3/rust-numpy), [ndarray](https://github.com/rust-ndarray/ndarray) and their Rust dependencies | MIT / Apache-2.0 / BSD-2-Clause | Compiled into the optional `simurg_core` extension (`rust/simurg-core`); list with `cargo tree` | +| [UCSC Genome Browser](https://genome.ucsc.edu/) `cytoBand` tables (hg38, hg19) and `ncbiRefSeq` (RefSeq transcripts, NCBI) | Free for all uses (UCSC data terms; RefSeq is public domain) | `simurg/cn/data/cytoBand_GRCh3{7,8}.txt.gz` and the derived centromere/acrocentric mask BEDs are shipped; `scripts/build_annotation_bundle.sh` downloads cytoBand + ncbiRefSeq for gene lists | | [illuminaio](https://github.com/HenrikBengtsson/illuminaio) | GPL-2 | Documentation of IDAT field codes only; oracle via subprocess | ## BSD-2-Clause (BeadArrayFiles, verbatim from its README) diff --git a/scripts/build_annotation_bundle.sh b/scripts/build_annotation_bundle.sh new file mode 100755 index 0000000..3351760 --- /dev/null +++ b/scripts/build_annotation_bundle.sh @@ -0,0 +1,32 @@ +#!/usr/bin/env bash +# Build the cyto annotation bundle (spec CYTO-036, IN-054) from public UCSC tables: +# cytoBand (ISCN bands) and ncbiRefSeq (RefSeq transcripts with gene symbols, from NCBI RefSeq). +# Output: simurg-cyto-annot---.zip with cytoBand.txt.gz, genes.tsv.gz and manifest.json. +# bash scripts/build_annotation_bundle.sh GRCh38|GRCh37 [OUT_DIR] +set -euo pipefail +BUILD=${1:?usage: build_annotation_bundle.sh GRCh38|GRCh37 [OUT_DIR]} +OUT=${2:-data/annotation} +case $BUILD in GRCh38) UCSC=hg38;; GRCh37) UCSC=hg19;; *) echo "unknown build $BUILD"; exit 2;; esac +BASE=https://hgdownload.soe.ucsc.edu/goldenPath/$UCSC/database +W=$(mktemp -d); trap 'rm -rf "$W"' EXIT +mkdir -p "$OUT" +for t in cytoBand ncbiRefSeq; do curl -sSfL -o "$W/$t.txt.gz" "$BASE/$t.txt.gz"; done +curl -sSfL -o "$W/ncbiRefSeq.sql" "$BASE/ncbiRefSeq.sql" || true +# ncbiRefSeq columns: bin name chrom strand txStart txEnd cdsStart cdsEnd exonCount exonStarts exonEnds score name2 ... +zcat "$W/ncbiRefSeq.txt.gz" | awk -F'\t' 'BEGIN{OFS="\t"; print "chrom","start","end","gene","transcript"} + { c=$3; sub(/^chr/,"",c); if (c ~ /_/) next; print c, $5+1, $6, $13, $2 }' | gzip -9 > "$W/genes.tsv.gz" +DATE=$(date -u +%Y%m%d) +python3 - "$W" "$BUILD" "$UCSC" "$DATE" <<'PY' +import hashlib, json, sys, pathlib +w, build, ucsc, date = pathlib.Path(sys.argv[1]), *sys.argv[2:] +files = {n: hashlib.sha256((w / n).read_bytes()).hexdigest() for n in ("cytoBand.txt.gz", "genes.tsv.gz")} +m = {"schema": "simurg.cyto_annotation/1", "build": build, "created": date, + "sources": [{"name": "UCSC cytoBand", "version": ucsc, "release_date": None}, + {"name": "RefSeq (UCSC ncbiRefSeq)", "version": ucsc, "release_date": None}], + "files": files} +(w / "manifest.json").write_text(json.dumps(m, indent=1)) +PY +( cd "$W" && zip -q -X bundle.zip manifest.json cytoBand.txt.gz genes.tsv.gz ) +SHA=$(sha256sum "$W/bundle.zip" | cut -c1-12) +cp "$W/bundle.zip" "$OUT/simurg-cyto-annot-$BUILD-$DATE-$SHA.zip" +ls -la "$OUT/simurg-cyto-annot-$BUILD-$DATE-$SHA.zip" diff --git a/src/simurg/cli/main.py b/src/simurg/cli/main.py index 2813ee2..b9e618e 100644 --- a/src/simurg/cli/main.py +++ b/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") diff --git a/src/simurg/cn/annotate.py b/src/simurg/cn/annotate.py new file mode 100644 index 0000000..56ef77c --- /dev/null +++ b/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[] [](_)x`` per event, LOH as +``x2 hmz`` (gain-LOH ``x3 hmz``), mosaic as ``x~``; normal results ``arr[](1-22,X)x2`` (female), +``(1-22)x2,(XY)x1`` (male), ``(1-22)x2`` (sex unresolved). Traditional form ``del()()`` / +``dup()()`` 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 diff --git a/src/simurg/pipeline/cn.py b/src/simurg/pipeline/cn.py index 8d8214f..1cb5b18 100644 --- a/src/simurg/pipeline/cn.py +++ b/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__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] diff --git a/tests/unit/test_cn.py b/tests/unit/test_cn.py index 23067c7..ab2c172 100644 --- a/tests/unit/test_cn.py +++ b/tests/unit/test_cn.py @@ -171,3 +171,51 @@ def test_cn_build_and_call(tmp_path: Path) -> None: # CYTO-028 .. CYTO-032, CYT summ = (out / "cn_sample_summary.tsv").read_text().splitlines() assert summ[0].startswith("sample_id\tsample_name\tstatus") and len(summ) == 3 assert (out / "200000000001_R01C01.cn_segments.parquet").exists() + + +# ------------------------------------------------------------------------------------------------ annotate +def test_iscn_strings_and_genes() -> None: # CYTO-034 / CYTO-035 + from simurg.cn.annotate import Bands, Genes, iscn_sample + + bands = Bands.from_ucsc("chr1\t0\t2300000\tp36.33\tgneg\nchr1\t2300000\t5300000\tp36.32\tgpos25\n" + "chr9\t0\t2200000\tp24.3\tgneg\n") # fmt: skip + dele = {"chrom": "1", "start": 1_567_924, "end": 1_632_825, "call": "DEL", "call_type": "DEL", "copy_number": 1, + "mosaic": False} # fmt: skip + dup = {"chrom": "1", "start": 2_000_000, "end": 3_000_000, "call": "DUP", "call_type": "DUP", "copy_number": 3, + "mosaic": False} # fmt: skip + loh = {"chrom": "9", "start": 100, "end": 2_000_000, "call": "LOH", "call_type": "LOH", "copy_number": 2, + "mosaic": False} # fmt: skip + mos = {"chrom": "9", "start": 100, "end": 1_000, "call": "MDEL", "call_type": "DEL", "copy_number": 1, + "mosaic": True, "base_cn": 2} # fmt: skip + trad, micro = iscn_sample([loh, dup, dele], bands, "GRCh38", "F") + assert micro == ("arr[GRCh38] 1p36.33(1567924_1632825)x1,1p36.33p36.32(2000000_3000000)x3," + "9p24.3(100_2000000)x2 hmz") # fmt: skip + assert trad == "del(1)(p36.33p36.33),dup(1)(p36.33p36.32)" + assert iscn_sample([mos], bands, "GRCh38", "F")[1] == "arr[GRCh38] 9p24.3(100_1000)x2~1" + assert iscn_sample([], bands, "GRCh38", "M")[1] == "arr[GRCh38](1-22)x2,(XY)x1" + assert iscn_sample([], bands, "GRCh37", "F")[1] == "arr[GRCh37](1-22,X)x2" + assert iscn_sample([], bands, "GRCh38", "U")[1] == "arr[GRCh38](1-22)x2" + genes = Genes.from_rows([("1", 100, 5_000, "A"), ("1", 4_000, 4_500, "B"), ("1", 200, 300, "A"), + ("1", 9_000, 9_100, "C")]) # fmt: skip + assert genes.overlapping("1", 4_400, 8_999) == ["A", "B"] + + +def test_cn_annotate_end_to_end(tmp_path: Path) -> None: # CYTO-032 / CYTO-033 + bpm, fa, gtcs = _dataset(tmp_path) + model, out = tmp_path / "m.simurg-cyto", tmp_path / "cn" + run = CliRunner().invoke + assert run(app, ["cn", "annotate", "--vcfs", str(tmp_path)]).exit_code == 2 # E_CN_SIDECAR_MISSING + assert run(app, ["cn", "build-model", "--bpm", str(bpm), "--fasta", str(fa), "--out", str(model), + "--quiet"]).exit_code == 0 # fmt: skip + r = run(app, ["cn", "call", "--bpm", str(bpm), "--model", str(model), "--gtcs", str(gtcs), "--fasta", str(fa), + "--out", str(out), "--threads", "1", "--unknown-sex-policy", "female", "--quiet"]) # fmt: skip + assert r.exit_code == 0, r.output + r = run(app, ["cn", "annotate", "--vcfs", str(out), "--min-qual", "0", "--quiet"]) + assert r.exit_code == 0, r.output + doc = json.loads((out / "200000000001_R01C01.cyto.json").read_text()) + assert doc["schema"] == "simurg.cyto/1" and doc["annotation_type"] == "constitutional" + assert doc["iscn_microarray"].startswith("arr[") + assert any(e["call_type"] == "DEL" and e["chrom"] == "1" for e in doc["events"]) + assert {c["id"] for c in doc["chromosomes"]} == {"1", "2"} + normal = json.loads((out / "200000000001_R02C01.cyto.json").read_text()) + assert normal["iscn_microarray"] in ("arr[GRCh38](1-22,X)x2", "arr[GRCh37](1-22,X)x2") or normal["events"]