PGx: all CPIC gene kinds — G6PD, variant-based, HLA tags, *0 CN genes (PGX-006/007/021)

- DB build resolves the PGX-007 tag variants from the manifest's probes; coverage over all gene kinds.
- Variant-based genes reported per variant (rs.. REF/ALT | Reference) with the CPIC lowest-activity diplotype;
  tagging genes as positive/negative with the population caveat; GSTM1/GSTT1/UGT2B17 *0 from targeted CN.
- Case-insensitive phenotype keys (MT-RNR1). GeT-RM check extended to VKORC1 and the CN genes:
  131/144 gene calls (docs/validation/V2_pgx_star.md).

src/simurg/pipeline/pgx_db.py

@@ -9,6 +9,7 @@ from __future__ import annotations
 
 import hashlib
 import json
+import re
 import sqlite3
 from dataclasses import dataclass, field
 from pathlib import Path
@@ -70,6 +71,12 @@ def _load_cpic(d: Path) -> dict[str, list[dict[str, Any]]]:
     return t
 
 
+# PGX-007 tagging variants (dbSNP plus-strand alleles); positions are resolved from the manifest's probes
+TAGS = [("HLA-A", "*31:01", "rs1061235", "A", "T"), ("HLA-B", "*15:02", "rs144012689", "T", "A"),
+        ("HLA-B", "*57:01", "rs2395029", "T", "G"), ("HLA-B", "*58:01", "rs9263726", "G", "A")]  # fmt: skip
+TAG_CAVEAT = "tag reliability varies by population"
+
+
 def _probe_sites(
     bpm_path: Path,
     csv_path: Path | None,
@@ -77,6 +84,7 @@ def _probe_sites(
     cmap: Any,
     windows: dict[str, list[tuple[int, int]]],
     rep: Reporter,
+    rsids: set[str] | None = None,
 ) -> list[Any]:
     """Plus-strand, normalised sites of the manifest probes inside the gene windows (the SNV VCF representation)."""
     from simurg.formats.bpm import read_bpm
@@ -96,13 +104,18 @@ def _probe_sites(
                       csvm.source_seqs if csvm else [""] * bpm.num_loci, bpm.assay_types.tolist(),
                       bpm.intensity_only.tolist())  # fmt: skip
     out = []
+    tok = re.compile(r"rs\d+")
     for i in range(bpm.num_loci):
         c = cmap.get(bpm.chroms[i])
-        if c is None or c not in windows:
-            continue
-        p = int(bpm.map_infos[i])
-        if not any(a - GENE_PAD <= p <= b + GENE_PAD for a, b in windows[c]):
+        if rsids is not None:
+            if c is None or not (set(tok.findall(bpm.names[i])) & rsids):
+                continue
+        elif c is None or c not in windows:
             continue
+        else:
+            p = int(bpm.map_infos[i])
+            if not any(a - GENE_PAD <= p <= b + GENE_PAD for a, b in windows[c]):
+                continue
         site, _ = map_locus(i, mv, genome.fetch, cmap.get, genome.build or "GRCh38")
         if site is not None:
             out.append(site)
@@ -168,6 +181,7 @@ def run_db_build(o: PgxDbOptions, rep: Reporter) -> int:
         windows.setdefault(av[3], []).append((av[4], av[4] + len(av[5])))
     var_keys = {(av[3], av[4], av[5], av[6]) for av in rows["allele_variants"]}
     manifests = []
+    tag_rows: dict[tuple[str, str], tuple[Any, ...]] = {}
     csvs = dict(zip(o.manifests, o.csvs, strict=False))
     for bpm_path in o.manifests:
         name = read_bpm(bpm_path).manifest_name
@@ -179,6 +193,27 @@ def run_db_build(o: PgxDbOptions, rep: Reporter) -> int:
                     rows["probe_map"].append((name, ps.name, ps.contig, ps.pos, ps.ref, alt))
                     n += 1
         rep.info(f"{name}: {n} probe-to-variant link(s)")
+        comp = str.maketrans("ACGT", "TGCA")
+        for ps in _probe_sites(bpm_path, csvs.get(bpm_path), genome, cmap, {}, rep, {t[2] for t in TAGS}):
+            for gene, allele, rsid, t_ref, t_alt in TAGS:
+                if rsid not in re.findall(r"rs\d+", ps.name) or len(ps.ref) != 1:
+                    continue
+                if ps.ref == t_ref:
+                    alt = t_alt
+                elif ps.ref == t_ref.translate(comp):
+                    alt = t_alt.translate(comp)
+                else:
+                    continue
+                if alt not in ps.alleles:
+                    continue
+                hg = full_hgvs(ps.contig, f"g.{ps.pos}{ps.ref}>{alt}")
+                tag_rows[(gene, allele)] = (gene, allele, "GRCh38", ps.contig, ps.pos, ps.ref, alt, hg, rsid, -1, alt)
+                rows["probe_map"].append((name, ps.name, ps.contig, ps.pos, ps.ref, alt))
+    if tag_rows:  # PGX-007: tag alleles become defined single-variant alleles
+        rows["allele_variants"] += list(tag_rows.values())
+        rows["alleles"] = [(g, a, k, f, act, sv, 1 if (g, a) in tag_rows else d) for g, a, k, f, act, sv, d in
+                           rows["alleles"]]  # fmt: skip
+        rep.info(f"tagging variants resolved: {', '.join(f'{g} {a}' for g, a in sorted(tag_rows))}")
     for mp in o.masks:
         for ln in mp.read_text(encoding="utf-8").splitlines():
             if ln.strip() and not ln.startswith("#"):
@@ -240,7 +275,9 @@ def coverage_rows(db: Path, manifest: str) -> list[dict[str, Any]]:
     kinds = {(g, a): (k, d) for g, a, k, d in con.execute("SELECT gene, allele, kind, defined FROM alleles")}
     con.close()
     out = []
-    for gene, m in sorted(load_models(db, "global").items()):
+    from simurg.pipeline.pgx_star import ALL_BASES
+
+    for gene, m in sorted(load_models(db, "global", None, ALL_BASES).items()):
         usable = np.array([bool(links.get(v)) and not all(pr in masked for pr in links[v]) for v in m.variants])
         masked_v = np.array([bool(links.get(v)) and all(pr in masked for pr in links[v]) for v in m.variants])
         sig: dict[bytes, list[int]] = {}