kemal
004739da8e
- scripts/pgx/wgs_star_truth.py: star alleles from the 1000G 30x phased panel with the DB definitions (wildcard-aware haplotype assignment; checked against GeT-RM on the demo CEPH individuals). - validate_star_getrm.py: --all-genes, --individuals, --coverage (projection onto what the array can see); replicate_concordance.py. - 120 individuals x 3 GSA v1 arrays: array-resolvable 1637/1710 (95.7 %), top 1223/1800 (chip lacks many defining variants), GeT-RM subset 28/35, replicates 99.2 % without the failed array (V2_pgx_star.md, K-029).
200 lines
10 KiB
Python
200 lines
10 KiB
Python
#!/usr/bin/env python3
|
|
"""Compare `simurg pgx star-allele call` with GeT-RM consensus diplotypes (PGX-032).
|
|
|
|
python scripts/pgx/validate_star_getrm.py --calls STAR_DIR/pgx_calls.json --truth getrm_truth.tsv \
|
|
--identity scripts/pgx/demo36_identity.tsv --out summary.json
|
|
|
|
GeT-RM notation is normalised per allele: "*2(*35)" / "*2 (*35)" accept *2 or *35; a trailing sub-allele letter
|
|
("*1F") also matches the core allele ("*1") when the call has no letter; "(*68)+*4" tandems and "xN" duplications
|
|
reduce to the core allele (structural forms are reported separately); "(*28 + *60)" keeps the star alleles the DB
|
|
defines. "A/B or C/D" accepts either. UGT1A1 and SLCO1B1 are reported but left out of the total: GeT-RM uses pre-2023
|
|
names (UGT1A1 *60, SLCO1B1 *1B/*17/*21) that the CPIC definitions no longer carry. A sample's gene counts as **top**
|
|
concordant when the top solution matches, and as **equivalent**
|
|
concordant when the truth is among the array-equivalent members of any feasible solution. Replicate arrays of one
|
|
individual are counted once (the first).
|
|
"""
|
|
|
|
from __future__ import annotations
|
|
|
|
import argparse
|
|
import csv
|
|
import json
|
|
import re
|
|
from collections import defaultdict
|
|
from pathlib import Path
|
|
|
|
|
|
def truth_options(text: str) -> list[set[str]] | None:
|
|
"""Per haplotype, the set of acceptable core alleles; None when the text is not a diplotype."""
|
|
t = text.replace(" ", "")
|
|
if "/" not in t:
|
|
return None
|
|
out = []
|
|
for side in t.split("/", 1):
|
|
opts = set(re.findall(r"\*\d+[A-Z]?", side))
|
|
if not opts:
|
|
return None
|
|
cores = {re.sub(r"x\d+|xN", "", o) for o in opts}
|
|
cores |= {re.sub(r"[A-Z]$", "", o) for o in cores}
|
|
if "+" in side: # tandem / hybrid: the non-parenthesised last unit is the core allele
|
|
last = re.findall(r"\*\d+[A-Z]?", side.split("+")[-1])
|
|
cores |= set(last)
|
|
out.append(cores)
|
|
return out
|
|
|
|
|
|
def _gstm1(t: str) -> int: # GeT-RM GSTM1: DEL/DEL = no copy; *A/*B assays do not see a heterozygous deletion
|
|
return 0 if t.replace(" ", "").upper() == "DEL/DEL" else 1
|
|
|
|
|
|
def _gstt1(t: str) -> int: # "0" null, "(*A/0)" one null allele, "(*A/*A)" two copies
|
|
t = t.replace(" ", "")
|
|
return 0 if t in ("0", "0/0") else (1 if "/0" in t or "0/" in t else 2)
|
|
|
|
|
|
def _ugt2b17(t: str) -> int: # *2 = whole-gene deletion
|
|
return 2 - t.count("*2")
|
|
|
|
|
|
CN_TRUTH = {"GSTM1": _gstm1, "GSTT1": _gstt1, "UGT2B17": _ugt2b17} # capped: GSTM1 at 1, the others at 2
|
|
# GeT-RM uses older names or definitions than the current CPIC/PharmVar ones: UGT1A1 *60 and SLCO1B1 *1B/*17/*21
|
|
# (pre-2023), CYP1A2 lettered alleles (*1F …), CYP2C8 *19 (PharmVar now defines rs2071426 alone as *19; GeT-RM *1)
|
|
LEGACY_NOMENCLATURE = {"UGT1A1", "SLCO1B1", "CYP1A2", "CYP2C8"}
|
|
|
|
|
|
def truth_alternatives(text: str) -> list[list[set[str]]]:
|
|
"""GeT-RM "A/B or C/D" → one option list per alternative diplotype."""
|
|
alts = [truth_options(x) for x in re.split(r"\bor\b", text)]
|
|
return [x for x in alts if x is not None]
|
|
|
|
|
|
def matches(call: tuple[str, str], truth: list[list[set[str]]]) -> bool:
|
|
a, b = (re.sub(r"x\d+$", "", x.split("+")[-1]) for x in call) # tandem "*68+*4" → core "*4"
|
|
return any((a in t[0] and b in t[1]) or (a in t[1] and b in t[0]) for t in truth)
|
|
|
|
|
|
def main() -> int:
|
|
ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
|
|
ap.add_argument("--calls", required=True, type=Path)
|
|
ap.add_argument("--truth", required=True, type=Path)
|
|
ap.add_argument("--identity", required=True, type=Path)
|
|
ap.add_argument("--out", required=True, type=Path)
|
|
ap.add_argument("--all-genes", action="store_true",
|
|
help="truth with the same definitions (wgs_star_truth.py): no legacy-name exclusions") # fmt: skip
|
|
ap.add_argument("--individuals", type=Path, help="only these individuals (one per line)")
|
|
ap.add_argument("--coverage", type=Path, help="pgx_coverage.tsv of the array: also report array-resolvable "
|
|
"concordance (both diplotypes projected onto the alleles the array can tell apart)") # fmt: skip
|
|
a = ap.parse_args()
|
|
if a.all_genes:
|
|
LEGACY_NOMENCLATURE.clear()
|
|
only = set(a.individuals.read_text().split()) if a.individuals else None
|
|
proj: dict[tuple[str, str], str] = {}
|
|
refs: dict[str, str] = {}
|
|
if a.coverage:
|
|
with a.coverage.open() as fh:
|
|
for r in csv.DictReader((ln for ln in fh if not ln.startswith("#")), delimiter="\t"):
|
|
if r["allele_kind"] == "reference":
|
|
refs[r["gene"]] = r["allele"]
|
|
st = r["status"]
|
|
proj[(r["gene"], r["allele"])] = st.split(":", 1)[1] if st.startswith("equivalent_to:") else (
|
|
r["allele"] if st == "callable" else "") # fmt: skip
|
|
|
|
def project(gene: str, allele: str) -> str:
|
|
core = re.sub(r"x\d+$|xN$", "", allele.split("+")[-1])
|
|
v = proj.get((gene, core), core)
|
|
if v:
|
|
return v
|
|
ref = refs.get(gene, core) # not covered: indistinguishable from the reference's class
|
|
return proj.get((gene, ref)) or ref
|
|
|
|
with a.identity.open() as fh:
|
|
ident = {r["sample_id"]: r["individual"] for r in csv.DictReader((ln for ln in fh if not ln.startswith("#")),
|
|
delimiter="\t")} # fmt: skip
|
|
truth: dict[tuple[str, str], str] = {}
|
|
with a.truth.open() as fh:
|
|
for r in csv.DictReader(fh, delimiter="\t"):
|
|
if r.get("flag") == "incomplete" or (only is not None and r["individual"] not in only):
|
|
continue
|
|
truth[(r["individual"], r["gene"])] = r["diplotype"]
|
|
calls = json.loads(a.calls.read_text())
|
|
per_gene: dict[str, dict[str, int]] = defaultdict(lambda: defaultdict(int))
|
|
detail = []
|
|
seen = set()
|
|
for s in calls["samples"]:
|
|
ind = ident.get(s["sample_id"])
|
|
if ind is None or ind in seen:
|
|
continue
|
|
seen.add(ind)
|
|
for g in s["genes"]:
|
|
t = truth.get((ind, g["gene"]))
|
|
if t is None:
|
|
continue
|
|
st = per_gene[g["gene"]]
|
|
if g["gene"] in CN_TRUTH: # whole-gene deletion genes: compare copy-number classes (PGX-006)
|
|
want_cn = CN_TRUTH[g["gene"]](t)
|
|
cn = g.get("copy_number")
|
|
got_cn = None if cn is None else (min(cn, 2) if g["gene"] != "GSTM1" else min(cn, 1))
|
|
st["n"] += 1
|
|
ok = want_cn is not None and want_cn == got_cn
|
|
st["top"] += ok
|
|
st["equivalent"] += ok
|
|
if not ok:
|
|
detail.append({"individual": ind, "gene": g["gene"], "truth": t, "call": g["genotype"],
|
|
"within_equivalents": False}) # fmt: skip
|
|
continue
|
|
if g["gene"] == "VKORC1": # c.-1639G>A (minus strand) = rs9923231 C>T: count A = count T
|
|
want = t.count("A")
|
|
got = g["genotype"].split("/")[-1].count("T") if g["genotype"] != "Reference" else 0
|
|
st["n"] += 1
|
|
st["top"] += want == got
|
|
st["equivalent"] += want == got
|
|
if want != got:
|
|
detail.append({"individual": ind, "gene": "VKORC1", "truth": t, "call": g["genotype"],
|
|
"within_equivalents": False}) # fmt: skip
|
|
continue
|
|
opts = truth_alternatives(t)
|
|
if not opts:
|
|
st["not_comparable"] += 1
|
|
continue
|
|
st["n"] += 1
|
|
top = g["candidates"][0]["genotype"].split("/") if g["candidates"] else ["?", "?"]
|
|
ok_top = len(top) == 2 and matches((top[0], top[1]), opts)
|
|
feas = [c for c in g["candidates"] if c["genotype"] in g["feasible_solutions"]]
|
|
ok_eq = ok_top or any(matches(tuple(e.split("/")[:2]), opts) for c in feas for e in c["equivalents"]
|
|
if len(e.split("/")) >= 2) # fmt: skip
|
|
st["top"] += ok_top
|
|
st["equivalent"] += ok_eq
|
|
if a.coverage and len(top) == 2 and len(opts) == 1 and all(len(x) == 1 for x in opts[0]):
|
|
tt = sorted(project(g["gene"], next(iter(x))) for x in opts[0])
|
|
cc = sorted(project(g["gene"], x) for x in top)
|
|
st["array_resolvable_n"] += 1
|
|
st["array_resolvable"] += tt == cc
|
|
if not ok_top:
|
|
detail.append({"individual": ind, "gene": g["gene"], "truth": t, "call": g["genotype"],
|
|
"within_equivalents": ok_eq}) # fmt: skip
|
|
comparable = {k: v for k, v in per_gene.items() if k not in LEGACY_NOMENCLATURE}
|
|
tot_ar = [
|
|
sum(v.get("array_resolvable", 0) for v in comparable.values()),
|
|
sum(v.get("array_resolvable_n", 0) for v in comparable.values()),
|
|
]
|
|
summary = {"individuals": len(seen), "genes": {k: dict(v) for k, v in sorted(per_gene.items())},
|
|
"array_resolvable_total": {"concordant": tot_ar[0], "n": tot_ar[1]},
|
|
"legacy_nomenclature_genes": sorted(LEGACY_NOMENCLATURE),
|
|
"total": {k: sum(v.get(k, 0) for v in comparable.values()) for k in ("n", "top", "equivalent")},
|
|
"discordant": detail} # fmt: skip
|
|
a.out.write_text(json.dumps(summary, indent=1))
|
|
for g, v in sorted(per_gene.items()):
|
|
print(f"{g:8s} top {v['top']}/{v['n']} within equivalents {v['equivalent']}/{v['n']}"
|
|
+ (f" (not comparable {v['not_comparable']})" if v.get("not_comparable") else "")
|
|
+ (f" array-resolvable {v.get('array_resolvable', 0)}/{v.get('array_resolvable_n', 0)}"
|
|
if a.coverage else "")
|
|
+ (" [legacy nomenclature]" if g in LEGACY_NOMENCLATURE else "")) # fmt: skip
|
|
print("total (excluding legacy-nomenclature genes)", summary["total"])
|
|
if a.coverage:
|
|
print("array-resolvable total", summary["array_resolvable_total"])
|
|
return 0
|
|
|
|
|
|
if __name__ == "__main__":
|
|
raise SystemExit(main())
|