PGx validation on GSE165845: WGS-derived star truth, 2-fold masks, array-resolvable concordance

- 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).

Co-Authored-By: Claude Opus 5.5 <[email protected]>
Claude-Session: https://claude.ai/code/session_01814GKVzZVL3FuF3DeiZXhv

scripts/pgx/replicate_concordance.py

@@ -0,0 +1,53 @@
+#!/usr/bin/env python3
+"""Replicate concordance of `pgx star-allele call`: the top genotype of every gene compared between the replicate
+arrays of one individual (PGX-029 reproducibility on real replicates).
+
+    python scripts/pgx/replicate_concordance.py --calls STAR_DIR/pgx_calls.json --identity identity.tsv --out rep.json
+"""
+
+from __future__ import annotations
+
+import argparse
+import csv
+import itertools
+import json
+from collections import defaultdict
+from pathlib import Path
+
+
+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("--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
+    calls: dict[str, dict[str, dict[str, str]]] = defaultdict(dict)
+    for s in json.loads(a.calls.read_text())["samples"]:
+        ind = ident.get(s["sample_id"])
+        if ind:
+            calls[ind][s["sample_id"]] = {g["gene"]: (g["candidates"][0]["genotype"] if g["candidates"] else "")
+                                          for g in s["genes"]}  # fmt: skip
+    per_gene: dict[str, list[int]] = defaultdict(lambda: [0, 0])
+    discord = []
+    for ind, reps in sorted(calls.items()):
+        for (_s1, c1), (_s2, c2) in itertools.combinations(sorted(reps.items()), 2):
+            for gene in sorted(set(c1) & set(c2)):
+                ok = c1[gene] == c2[gene]
+                per_gene[gene][0] += ok
+                per_gene[gene][1] += 1
+                if not ok:
+                    discord.append({"individual": ind, "gene": gene, "a": c1[gene], "b": c2[gene]})
+    tot = [sum(v[0] for v in per_gene.values()), sum(v[1] for v in per_gene.values())]
+    out = {"pairs_by_gene": {g: {"concordant": v[0], "pairs": v[1]} for g, v in sorted(per_gene.items())},
+           "total": {"concordant": tot[0], "pairs": tot[1], "rate": round(tot[0] / tot[1], 4) if tot[1] else None},
+           "discordant": discord}  # fmt: skip
+    a.out.write_text(json.dumps(out, indent=1))
+    print(json.dumps(out["total"]), f"{len(discord)} discordant gene pair(s)")
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())