#!/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 LEGACY_NOMENCLATURE = {"UGT1A1", "SLCO1B1"} # GeT-RM uses pre-2023 names (*60, SLCO1B1 *1B/*17/*21) not in CPIC 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) a = ap.parse_args() 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"): 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 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} summary = {"individuals": len(seen), "genes": {k: dict(v) for k, v in sorted(per_gene.items())}, "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 "") + (" [legacy nomenclature]" if g in LEGACY_NOMENCLATURE else "")) # fmt: skip print("total (excluding legacy-nomenclature genes)", summary["total"]) return 0 if __name__ == "__main__": raise SystemExit(main())