PGx CN: pgx copy-number train|call (CLI-031/032, PCN-001..017, IN-050), region BED v1, demo identities (K-028)

- simurg/pgx/cn.py: per-region CN mixture, mu_c = delta + gamma*log2((c+1)/3), free lower CN-0 component,
  Hardy-Weinberg-constrained weights over haplotype copy numbers, EM from a grid of label starts.
- pipeline/pgx_cn.py: model file (.simurg-pgxcn), batch rules (PCN-004), per-sample CN VCF + tbi, batch and
  sample summaries.
- The GSA v4 PGx demo arrays are 30 CEPH/1000G individuals (scripts/pgx/demo36_identity.tsv); against NYGC 30x
  SV RD_CN 224/240 region calls agree (GSTM1, UGT2B17 30/30), in-sample (docs/validation/V2_pgx_cn.md).

scripts/pgx/validate_cn_truth.py

@@ -0,0 +1,72 @@
+#!/usr/bin/env python3
+"""Validate `simurg pgx copy-number call` against the NYGC 30x SV callset (1KGP_3202 freeze V3) — PCN-017.
+
+    python scripts/pgx/validate_cn_truth.py --calls CALL_DIR --identity scripts/pgx/demo36_identity.tsv \
+        --sv 1KGP_3202.sv_freeze_V3.vcf.gz --bcftools BCFTOOLS --out summary.json
+
+Truth per region = RD_CN of one SV record per region (below, chosen as the record spanning the whole region; for
+GSTM1 the read-depth DUP/CNV record, whose RD_CN gives 0/1/2 in Hardy–Weinberg proportions — the DEL record there
+reports CN 1 for 27/30 individuals). CYP2D6 sub-regions share one gene-level record, so exon-9 conversions (*36) and
+2D7-2D6 hybrids (*13, *68) legitimately differ from it. GSTT1 lies on an alt contig with no SV calls: only its CN
+distribution is reported. Replicate arrays of one individual are counted once (the first array).
+"""
+
+from __future__ import annotations
+
+import argparse
+import collections
+import csv
+import json
+import subprocess
+from pathlib import Path
+
+TRUTH_RECORD = {"GSTM1": "HGSV_9691", "UGT2B17": "HGSV_63272", "SULT1A1": "HGSV_207805",
+                "CYP2A6.intron.7": "HGSV_232737", "CYP2A6.exon.1": "HGSV_232737", "CYP2D6.exon.9": "HGSV_247317",
+                "CYP2D6.intron.2": "HGSV_247317", "CYP2D6.p5": "HGSV_247317"}  # fmt: skip
+
+
+def main() -> int:
+    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
+    ap.add_argument("--calls", required=True, type=Path)
+    ap.add_argument("--identity", required=True, type=Path)
+    ap.add_argument("--sv", required=True)
+    ap.add_argument("--bcftools", default="bcftools")
+    ap.add_argument("--out", required=True, type=Path)
+    a = ap.parse_args()
+    ident = {}
+    with a.identity.open() as fh:
+        for r in csv.DictReader((ln for ln in fh if not ln.startswith("#")), delimiter="\t"):
+            ident[r["sample_id"]] = r["individual"]
+    calls: dict[str, dict[str, str]] = {}
+    with (a.calls / "pgx_cn_sample_summary.tsv").open() as fh:
+        for r in csv.DictReader(fh, delimiter="\t"):
+            ind = ident.get(r["sample_id"])
+            if ind and ind not in calls:
+                calls[ind] = {k[3:]: v for k, v in r.items() if k.startswith("cn_")}
+    ids = sorted(set(TRUTH_RECORD.values()))
+    q = subprocess.run([a.bcftools, "query", "-H", "-s", ",".join(sorted(calls)), "--force-samples", "-i",
+                        "ID=@/dev/stdin", "-f", "%ID[\t%RD_CN]\n", a.sv], input="\n".join(ids), capture_output=True,
+                       text=True, check=True).stdout.splitlines()  # fmt: skip
+    samples = [h.split("]")[1].split(":")[0] for h in q[0].split("\t")[1:]]
+    truth = {ln.split("\t")[0]: dict(zip(samples, ln.split("\t")[1:], strict=True)) for ln in q[1:]}
+    res: dict[str, dict[str, object]] = {}
+    tot = [0, 0]
+    for region, rid in TRUTH_RECORD.items():
+        pairs = [(truth[rid][ind], c[region]) for ind, c in calls.items() if ind in truth.get(rid, {})]
+        pairs = [(t, c) for t, c in pairs if t not in (".", "")]
+        ok = sum(t == c for t, c in pairs)
+        tot[0] += ok
+        tot[1] += len(pairs)
+        res[region] = {"record": rid, "concordant": ok, "n": len(pairs),
+                       "rate": round(ok / len(pairs), 4) if pairs else None,
+                       "truth_to_call_mismatches": dict(collections.Counter(f"{t}->{c}" for t, c in pairs if t != c))}  # fmt: skip
+    summary = {"individuals": len(calls), "regions": res, "overall": {"concordant": tot[0], "n": tot[1],
+               "rate": round(tot[0] / tot[1], 4) if tot[1] else None},
+               "gstt1_cn_distribution": dict(sorted(collections.Counter(c.get("GSTT1", ".") for c in calls.values()).items()))}  # fmt: skip
+    a.out.write_text(json.dumps(summary, indent=1))
+    print(json.dumps(summary, indent=1))
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())