kemal
cc8c934167
- 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).
73 lines
3.8 KiB
Python
73 lines
3.8 KiB
Python
#!/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())
|