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

scripts/pgx/wgs_star_truth.py

@@ -0,0 +1,130 @@
+#!/usr/bin/env python3
+"""Star-allele truth from 1000 Genomes 30x phased genotypes with the PGx DB definitions (SNV/indel alleles only).
+
+    python scripts/pgx/wgs_star_truth.py --db DB.sqlite --individuals ids.txt --out wgs_truth.tsv [--genes G,G]
+
+For every star-allele gene of the DB the defining variants are fetched from the NYGC phased panel (pysam over HTTPS,
+one window per gene; ``uv pip install pysam``). Each phased haplotype gets the allele whose defining variants it
+carries completely and that explains most of the haplotype's carried variants — its defining ones plus those at its
+wildcard (IUPAC) positions, so a *4 haplotype carrying the *10 variants stays *4 (ties: more defining variants,
+higher CPIC/PharmVar frequency, natural order); a haplotype carrying none is the reference allele. A variant absent
+from the filtered panel counts as reference only when rare (max allele frequency < 5 %); a gene with a common
+defining variant missing from the panel is written with ``incomplete`` (not comparable). Structural alleles (CYP2D6
+*5, xN, hybrids; CYP2A6 *4) are invisible to this truth: the diplotypes are SNV-level only.
+"""
+
+from __future__ import annotations
+
+import argparse
+import os
+import re
+import sqlite3
+from collections import defaultdict
+from pathlib import Path
+
+PANEL = (
+    "https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/"
+    "20220422_3202_phased_SNV_INDEL_SV/1kGP_high_coverage_Illumina.{c}.filtered.SNV_INDEL_SV_phased_panel.vcf.gz"
+)
+
+
+def natural(a: str) -> tuple[object, ...]:
+    return tuple(int(x) if x.isdigit() else x for x in re.split(r"(\d+)", a))
+
+
+def main() -> int:
+    import pysam
+
+    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
+    ap.add_argument("--db", required=True, type=Path)
+    ap.add_argument("--individuals", required=True, type=Path)
+    ap.add_argument("--out", required=True, type=Path)
+    ap.add_argument("--genes", default="")
+    ap.add_argument("--cache", type=Path, default=Path("data/pgx/cache"))
+    a = ap.parse_args()
+    inds = sorted(set(a.individuals.read_text().split()))
+    out_path = a.out.resolve()
+    con = sqlite3.connect(a.db)
+    genes = [g for (g,) in con.execute("SELECT gene FROM genes WHERE result_basis = 'star_allele' ORDER BY gene")]
+    if a.genes:
+        genes = [g for g in genes if g in a.genes.split(",")]
+    a.cache.mkdir(parents=True, exist_ok=True)
+    os.chdir(a.cache)
+    rows = []
+    for gene in genes:
+        q_alle = "SELECT allele FROM alleles WHERE gene = ? AND defined = 1 AND structural = 0"
+        alle = [x for (x,) in con.execute(q_alle, (gene,))]
+        q_ref = "SELECT allele FROM alleles WHERE gene = ? AND kind = 'reference'"
+        ref_allele = next((x for (x,) in con.execute(q_ref, (gene,))), "*1")
+        defs: dict[str, set[tuple[int, str, str]]] = defaultdict(set)
+        chrom = None
+        loc_keys: dict[int, set[tuple[int, str, str]]] = defaultdict(set)
+        q_var = "SELECT allele, chrom, pos, ref, alt, location_id FROM allele_variants WHERE gene = ?"
+        for allele, c, p, r, al, lid in con.execute(q_var, (gene,)):
+            defs[allele].add((p, r, al))
+            loc_keys[lid].add((p, r, al))
+            chrom = c
+        wild: dict[str, set[tuple[int, str, str]]] = defaultdict(set)
+        for allele, lid in con.execute("SELECT allele, location_id FROM allele_wildcards WHERE gene = ?", (gene,)):
+            wild[allele] |= loc_keys.get(lid, set())
+        if chrom is None:
+            continue
+        freq: dict[str, float] = defaultdict(float)
+        for allele, f in con.execute("SELECT allele, AVG(freq) FROM frequencies WHERE gene = ? GROUP BY allele",
+                                     (gene,)):  # fmt: skip
+            freq[allele] = float(f or 0)
+        vmax: dict[tuple[int, str, str], float] = defaultdict(float)
+        for allele, vs in defs.items():
+            for v in vs:
+                vmax[v] = max(vmax[v], freq.get(allele, 0.0))
+        keys = set().union(*defs.values()) if defs else set()
+        core = chrom.removeprefix("chr")
+        vf = pysam.VariantFile(PANEL.format(c=f"chr{core}"))
+        present = [s for s in inds if s in vf.header.samples]
+        vf.subset_samples(present)
+        lo = min(p for p, _, _ in keys) - 10
+        hi = max(p + len(r) for p, r, _ in keys) + 10
+        hap_alts: dict[tuple[str, int], set[tuple[int, str, str]]] = defaultdict(set)
+        seen = set()
+        for rec in vf.fetch(f"chr{core}", lo, hi):
+            for ai, alt in enumerate(rec.alts or (), start=1):
+                k = (rec.pos, rec.ref, alt)
+                if k not in keys:
+                    continue
+                seen.add(k)
+                for s in present:
+                    gt = rec.samples[s]["GT"]
+                    for h, x in enumerate(gt[:2]):
+                        if x == ai:
+                            hap_alts[(s, h)].add(k)
+        missing_common = [k for k in keys - seen if vmax[k] >= 0.05]
+        for s in present:
+            haps = []
+            for h in (0, 1):
+                carried = hap_alts[(s, h)]
+                # the allele that explains most carried variants (its defining ones plus those at its wildcard
+                # positions), then the most specific, then the most frequent
+                best = (0, 0, 0.0)
+                pick = ref_allele
+                for x in alle:
+                    d = defs.get(x, set())
+                    if not d or not d <= carried:
+                        continue
+                    score = (len(carried & (d | wild.get(x, set()))), len(d), freq.get(x, 0.0))
+                    if score > best or (score == best and natural(x) < natural(pick)):
+                        best, pick = score, x
+                haps.append(pick)
+            dip = "/".join(sorted(haps, key=natural))
+            rows.append((s, gene, dip, "incomplete" if missing_common else ""))
+        print(f"{gene}: {len(keys)} defining variant(s), {len(seen)} in the panel, "
+              f"{len(missing_common)} common missing", flush=True)  # fmt: skip
+    with out_path.open("w") as fh:
+        fh.write("individual\tgene\tdiplotype\tflag\n")
+        for r in sorted(rows):
+            fh.write("\t".join(r) + "\n")
+    print(f"{len(rows)} gene call(s) -> {out_path}")
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())