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

docs/known_unknowns.md

@@ -4,7 +4,7 @@ Plan §20 lists the open questions; this file records what the implementation ha
 
 | ID | Date | Topic | Finding | Status |
 |---|---|---|---|---|
-| K-029 | 2026-09-29 | PGx star alleles (v2) | (1) PGX-012 error model: ε(GS) = 0.5·(1 − GS)² was replaced by ε = 1e-3 (GS ≥ 0.15) / 0.05 — after probe masking 0 of 9 740 GSA v4 PGx calls disagreed with 1000G 30x at any GS, and the old model made correct low-GS calls lose to the prior (spec amended). (2) PGX-010 masks derived from 14 individuals (spec: ≥ 100); a variant absent from the filtered 1000G panel counts as hom-ref only if rare (CYP2B6 rs2279343 is filtered out, not absent). (3) CPIC core-allele definitions use IUPAC wildcards; alleles without definitions (NUDT15 *2, structural CYP2D6) are not callable from SNVs. (4) GeT-RM uses legacy UGT1A1/SLCO1B1 names. Result vs GeT-RM (16 CEPH): 100/109 gene calls, CYP2D6/2C19/2C9 41/45. (5) DPWG: no public machine-readable genotype→phenotype mapping; CPIC mappings are used for genes CPIC covers (as ClinPGx does), DPWG-only genes (CYP1A2, CYP3A4, F5) get no phenotype. (6) GeT-RM definitions are older than CPIC/PharmVar for UGT1A1, SLCO1B1, CYP1A2 and CYP2C8 (e.g. PharmVar CYP2C8 *19 = rs2071426 alone). Result vs GeT-RM (16 CEPH, 21 genes): 154/168 comparable gene calls. | Open: larger/independent truth set; KNMP DPWG tables if they become available (`docs/validation/V2_pgx_star.md`). |
+| K-029 | 2026-09-29 | PGx star alleles (v2) | (1) PGX-012 error model: ε(GS) = 0.5·(1 − GS)² was replaced by ε = 1e-3 (GS ≥ 0.15) / 0.05 — after probe masking 0 of 9 740 GSA v4 PGx calls disagreed with 1000G 30x at any GS, and the old model made correct low-GS calls lose to the prior (spec amended). (2) PGX-010 masks derived from 14 individuals (spec: ≥ 100); a variant absent from the filtered 1000G panel counts as hom-ref only if rare (CYP2B6 rs2279343 is filtered out, not absent). (3) CPIC core-allele definitions use IUPAC wildcards; alleles without definitions (NUDT15 *2, structural CYP2D6) are not callable from SNVs. (4) GeT-RM uses legacy UGT1A1/SLCO1B1 names. Result vs GeT-RM (16 CEPH): 100/109 gene calls, CYP2D6/2C19/2C9 41/45. (5) DPWG: no public machine-readable genotype→phenotype mapping; CPIC mappings are used for genes CPIC covers (as ClinPGx does), DPWG-only genes (CYP1A2, CYP3A4, F5) get no phenotype. (6) GeT-RM definitions are older than CPIC/PharmVar for UGT1A1, SLCO1B1, CYP1A2 and CYP2C8 (e.g. PharmVar CYP2C8 *19 = rs2071426 alone). Result vs GeT-RM (16 CEPH, 21 genes): 154/168 comparable gene calls. Independent check on GSE165845 (GSA v1 C2, 120 individuals × 3, WGS-derived truth, 2-fold masks): 95.7 % array-resolvable concordance, 67.9 % top-solution (the chip lacks many defining variants), replicate concordance 99.2 % without the failed array. | Open: KNMP DPWG tables if they become available (`docs/validation/V2_pgx_star.md`). |
 | K-028 | 2026-09-29 | PGx CN (v2) | Deviations from PCN-009: the per-probe reference is the training cohort's median stored in the model (a calling-batch median is not CN 2 where deletions are common, e.g. GSTM1 null ≈ 50 %); no GC-correction step; the CN mixture is ours: µ_c = δ + γ·log2((c+1)/3) per region with a free lower CN-0 component and Hardy–Weinberg-constrained weights (probe compression differs by region: γ ≈ 0.19 CYP2D6.p5 … 0.56 GSTT1). Interim validation on the GSA v4 PGx demo — whose 36 arrays turned out to be 30 CEPH/1000G individuals — against NYGC 30x SV RD_CN: 224/240 region calls (GSTM1, UGT2B17 30/30), but in-sample (trained on the same arrays, < 96 samples, priors chosen with this truth in view). | Open: validate on an independent GSA-PGx cohort; CYP2D6 via GeT-RM star alleles (`docs/validation/V2_pgx_cn.md`). |
 | K-027 | 2026-09-29 | CN/LOH (v2) | CYTO-014 noise gate: the summary `logr_dev` also measures cluster-file mismatch (GSE165845: 0.246, all samples `SampleQuality`) while the LRR used after the batch reference has SD ≈ 0.1; the gate now uses the SD of the normalised LRR (the JSON still reports the summary value). Spec: built-in mask/cytobands from UCSC; the annotation bundle uses UCSC ncbiRefSeq instead of RefSeq GFF + Ensembl GTF + HGNC. Real-data result: constitutional DEL/DUP calls 74–77 % in NYGC truth, mosaic calls mostly false positives. Pi tuning: mosaic prior / switch distance barely matter; CYTO-041 (noise-normalised BAF split, complete-loss promotion) cuts PASS mosaic calls by ~30 % without losing implanted ≥ 15 % mosaics. Two limits remain: (1) without phasing, < 20 % mosaic losses are not separable from regional BAF noise; (2) common CNVs (carrier frequency > 20 %) are almost never detected (1/236 DEL) because the median-LRR batch reference absorbs them; rare DEL (< 5 %) 38 %. | Open: rare-event sensitivity as primary metric; phasing-based mosaic BAF (`docs/validation/V2_cn.md`). |
 | K-026 | 2026-09-29 | Methylation QC (v2) | Spec deviations: vendor EPIC manifests put a 5th field on some control lines (≥ 4 fields accepted); detection p follows sesame's pOOBAH (1 − max(F(M), F(U))) rather than the ambiguous "1 − ECDF(M + U)"; minfi-style `Sentrix_ID`/`Sentrix_Position` sheet keys accepted. Against sesame 1.28.1 on minfiDataEPIC: β r = 0.99998 (max \|Δβ\| 0.011), detection calls 99.1 % identical, noob signals within 0.5 % (median) / 5 % (max). The systematic remainder is unexplained; it is not the negative controls (sesame did not attach the supplied control table). | Open; `docs/validation/V2_methylation.md`. Next step: dump sesame's background vectors and noob parameters to find the difference. |

docs/validation/V2_pgx_star.md

@@ -101,3 +101,54 @@ so the per-call error is 1e-3 (GS ≥ 0.15) and errors are handled by the probe
   for genes CPIC covers, so those genes reuse the CPIC tables (flag `dpwg_via_cpic_mapping`); DPWG-only genes
   (CYP1A2, CYP3A4, F5) have no public machine-readable mapping and get phenotype `null` with the flag
   `dpwg_mapping_not_available` (K-029).
+
+## Independent validation: GSE165845 (GSA v1 C2, 120 1000 Genomes individuals × 3 arrays)
+
+A second cohort, independent of the demo used for development (2026-09-29, all on the Pi, $0):
+
+* **Arrays:** the 360 GSE165845 GTCs (GSA-MD v1 chips, GSA v1 C2 manifest; K-024) → SNV VCFs; no targeted CN probes
+  on this chip, so star alleles are called with `--no-cn`. The chip carries little PGx content: 292 probe-to-variant
+  links (GSA v4 PGx: 5 280); `pgx coverage`: 495 alleles callable, 136 array-equivalent, 1 129 not covered — e.g.
+  no CYP2D6 *2/*4/*10 probes (rs16947, rs3892097, rs1065852 are absent from the manifest).
+* **Truth:** star alleles derived from the 1000G 30x phased panel with the same DB definitions
+  (`scripts/pgx/wgs_star_truth.py`; SNV/indel level — CYP2D6 *5, xN and hybrids invisible; CYP2B6 left out because
+  its common core variant rs2279343 is filtered from the panel). The derivation was checked against GeT-RM on the
+  demo's CEPH individuals (CYP2D6 14/15, CYP2C9/CYP3A4/CYP3A5 15/15, TPMT 16/16; the rest are nomenclature changes).
+* **Masks:** 2-fold — individuals split in halves; each half is called with the mask derived from the other half
+  (A → 0 probes masked, B → 1). After masking the calls agree with 1000G at 99 % or better in every GS bin.
+* **Metrics:** *top* = the top solution equals the WGS diplotype; *array-resolvable* = both diplotypes projected
+  onto the alleles this array can tell apart (`pgx coverage`: callable → itself, equivalent_to:X → X, not covered →
+  the reference's class) — the fair measure for a chip that lacks a gene's defining variants.
+
+| Gene | Top = WGS | Truth within feasible equivalents | Array-resolvable |
+|---|---|---|---|
+| CYP1A2 | 18/120 | 72/120 | 120/120 |
+| CYP2A13 | 97/120 | 119/120 | 117/120 |
+| CYP2A6 | 73/120 | 112/120 | 117/120 |
+| CYP2C19 | 100/120 | 116/120 | 110/120 |
+| CYP2C8 | 114/120 | 120/120 | 115/120 |
+| CYP2C9 | 95/120 | 108/120 | 116/120 |
+| CYP2D6 (SNV level) | 32/120 | 54/120 | 111/120 |
+| CYP3A4 | 112/120 | 117/120 | 117/120 |
+| CYP3A5 | 120/120 | 120/120 | 120/120 |
+| CYP4F2 | 43/120 | 86/120 | 96/120 |
+| NAT2 | 95/120 | 105/120 | 120/120 |
+| NUDT15 | 120/120 | 120/120 | 120/120 |
+| SLCO1B1 | 31/120 | 87/120 | 109/120 |
+| TPMT | 107/120 | 115/120 | 102/103 |
+| UGT1A1 | 66/120 | 120/120 | 47/47 |
+| **total** | **1223/1800 (67.9 %)** | | **1637/1710 (95.7 %)** |
+
+* Top-solution concordance is low where the chip lacks the defining variants: the call is then the most probable
+  member of an equivalence class the array cannot split (CYP2D6 *1/*2/*10, SLCO1B1 *1/*37, CYP1A2), so the prior
+  decides. Projected onto what the array can see, 95.7 % agree (half A 815/854, half B 822/856).
+* CYP4F2 (96/120) is the weakest: the array genotypes at its covered sites match 1000G exactly (rs3093200 180/180),
+  so the residue is definitional — *5 and *23 share the only covered variant, and the projection cannot separate
+  alleles whose other defining variants are off-array.
+* **GeT-RM** (7 of the 120 individuals, older consensus): 28/35 top, array-resolvable 29/30.
+* **Replicates:** the top genotype of all 70 genes agrees between replicate arrays of one individual in 24 916 of
+  25 200 pairs (98.9 %); 87 of the 284 discordant pairs come from NA20362's failed replicate (call rate 0.76, K-024);
+  without it 99.2 %.
+
+Files: `docs/validation/V2_pgx_star_gse165845.json`; scripts `wgs_star_truth.py`, `validate_star_getrm.py`
+(`--all-genes --individuals --coverage`), `replicate_concordance.py`.

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

scripts/pgx/validate_star_getrm.py

@@ -79,13 +79,42 @@ def main() -> int:
     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))
@@ -135,11 +164,21 @@ def main() -> int:
                                   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
@@ -147,8 +186,12 @@ def main() -> int:
     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
 
 

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