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

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

docs/validation/V2_pgx_star.md

@@ -37,7 +37,11 @@ Status: **implemented for CPIC star-allele genes, interim validation** (2026-09-
 | CYP4F2 | 8/9 | 8/9 |
 | NAT2 | 8/9 | 8/9 |
 | NUDT15 | 1/1 | 1/1 |
-| **total** | **100/109 (91.7 %)** | **103/109 (94.5 %)** |
+| VKORC1 (rs9923231, c.-1639G>A) | 13/15 | 13/15 |
+| GSTM1 (CN: DEL/DEL vs present) | 9/9 | 9/9 |
+| GSTT1 (CN 0/1/2) | 7/9 | 7/9 |
+| UGT2B17 (CN) | 2/2 | 2/2 |
+| **total** | **131/144 (91.0 %)** | **134/144 (93.1 %)** |
 
 PGX-032 target (≥ 90 % for CYP2D6/CYP2C19/CYP2C9 against truth): **41/45 = 91 %** (top solution).
 Details and every discordant call: `V2_pgx_star_getrm16.json`.
@@ -70,5 +74,10 @@ so the per-call error is 1e-3 (GS ≥ 0.15) and errors are handled by the probe
   `variants_not_covered.tsv`, per-allele supporting/unresolved sites for **both** alleles of every candidate,
   `variant_change` vs the top solution (`distant_alternative` above 3 changes), `unexplained_variants`, full HGVS on
   GRCh38 RefSeq accessions.
-* Not yet: named-haplotype (G6PD), variant-based and tagging genes; CYP2D6 hybrids (*13, *36, *68); PharmVar-only
+* All CPIC gene kinds are called: star alleles, G6PD (named haplotypes; males hemizygous from the VCF ploidy),
+  variant-based genes (ABCG2, CACNA1S, CFTR, DPYD, IFNL3, MT-RNR1, RYR1, VKORC1; genotype as `rs… REF/ALT` or
+  `Reference`, phenotype diplotype from the two lowest-activity carried variants), HLA-A/B tag variants (PGX-007,
+  positions resolved from the manifest's probes, flag `tag_caveat`), and GSTM1/GSTT1/UGT2B17 `*0` genotypes from
+  the targeted CN (PGX-006).
+* Not yet: CYP2D6 hybrids (*13, *36, *68); PharmVar-only
   genes (CYP2A6, CYP3A4, CYP1A2, …); DPWG tables.

docs/validation/V2_pgx_star_getrm16.json

@@ -32,6 +32,16 @@
    "top": 8,
    "equivalent": 8
   },
+  "GSTM1": {
+   "n": 9,
+   "top": 9,
+   "equivalent": 9
+  },
+  "GSTT1": {
+   "n": 9,
+   "top": 7,
+   "equivalent": 7
+  },
   "NAT2": {
    "n": 9,
    "top": 8,
@@ -56,6 +66,16 @@
    "n": 9,
    "top": 0,
    "equivalent": 1
+  },
+  "UGT2B17": {
+   "n": 2,
+   "top": 2,
+   "equivalent": 2
+  },
+  "VKORC1": {
+   "n": 15,
+   "top": 13,
+   "equivalent": 13
   }
  },
  "legacy_nomenclature_genes": [
@@ -63,9 +83,9 @@
   "UGT1A1"
  ],
  "total": {
-  "n": 109,
-  "top": 100,
-  "equivalent": 103
+  "n": 144,
+  "top": 131,
+  "equivalent": 134
  },
  "discordant": [
   {
@@ -110,6 +130,13 @@
    "call": "*1/*1",
    "within_equivalents": false
   },
+  {
+   "individual": "NA06991",
+   "gene": "VKORC1",
+   "truth": "AA",
+   "call": "rs9923231 T/T",
+   "within_equivalents": false
+  },
   {
    "individual": "NA12878",
    "gene": "CYP2D6",
@@ -145,6 +172,20 @@
    "call": "*28/*80",
    "within_equivalents": false
   },
+  {
+   "individual": "NA12003",
+   "gene": "VKORC1",
+   "truth": "AA",
+   "call": "rs9923231 T/T",
+   "within_equivalents": false
+  },
+  {
+   "individual": "NA12003",
+   "gene": "GSTT1",
+   "truth": "(*A/*A)",
+   "call": "*1/*0",
+   "within_equivalents": false
+  },
   {
    "individual": "NA07029",
    "gene": "CYP2B6",
@@ -229,6 +270,13 @@
    "call": "*80/*80",
    "within_equivalents": false
   },
+  {
+   "individual": "NA12813",
+   "gene": "GSTT1",
+   "truth": "(*A/*A)",
+   "call": "*1/*0",
+   "within_equivalents": false
+  },
   {
    "individual": "NA11881",
    "gene": "SLCO1B1",

scripts/pgx/getrm_truth.py

@@ -31,6 +31,10 @@ def main() -> int:
     rows = list(wb[wb.sheetnames[0]].iter_rows(values_only=True))
     head = [str(h or "").split("\n")[0].strip() for h in rows[1]]
     col = {g: head.index(g) for g in GENES if g in head}
+    vk = next((i for i, h in enumerate(rows[1]) if str(h or "").startswith("VKORC1") and "rs9923231" in str(h)
+               and "Reference" not in str(h)), None)  # fmt: skip
+    if vk is not None:
+        col["VKORC1"] = vk  # c.-1639G>A genotype, e.g. "GA"
     out = []
     for r in rows[2:]:
         if not r or not r[1] or str(r[1]).strip() not in inds:

scripts/pgx/validate_star_getrm.py

@@ -43,6 +43,20 @@ def truth_options(text: str) -> list[set[str]] | None:
     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
 
 
@@ -84,8 +98,30 @@ def main() -> int:
             t = truth.get((ind, g["gene"]))
             if t is None:
                 continue
-            opts = truth_alternatives(t)
             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

src/simurg/pgx/phenotype.py

@@ -48,16 +48,24 @@ def build_table(gene: str, rows: list[tuple[str, str | None, str]], alleles: lis
     by_key: dict[str, tuple[str, str | None]] = {}
     totals: dict[float, str] = {}
     for key, total, pheno in rows:
-        k = json.dumps(json.loads(key), sort_keys=True)
+        k = _canon(json.loads(key))
         by_key.setdefault(k, (pheno, total))
         t = _num(total)
         if t is not None:
             totals.setdefault(t, pheno)
-    func = {a: f for a, f, _ in alleles}
+    func: dict[str, str | None] = {a: f.lower() for a, f, _ in alleles if f}
     act = {a: _num(v) for a, _, v in alleles}
     return GuidelineTable(gene, by_key, sorted(totals.items()), func, act)
 
 
+def _canon(counts: dict[str, int]) -> str:
+    """Case-insensitive lookup key (CPIC mixes "Normal risk …" and "normal risk …", MT-RNR1)."""
+    merged: dict[str, int] = {}
+    for k, v in counts.items():
+        merged[k.lower()] = merged.get(k.lower(), 0) + v
+    return json.dumps(merged, sort_keys=True)
+
+
 def _split(allele: str) -> tuple[str, int]:
     m = re.fullmatch(r"(.+?)x(\d+|N)", allele)
     if not m:
@@ -80,7 +88,7 @@ def phenotype(t: GuidelineTable, haplotypes: tuple[str, ...]) -> PhenotypeResult
         for v in vals:
             s = f"{v:.1f}" if abs(v * 10 - round(v * 10)) < 1e-9 else f"{v:g}"
             counts[s] = counts.get(s, 0) + 1
-        key = json.dumps(counts, sort_keys=True)
+        key = _canon(counts)
         total = sum(vals)
         hit = t.by_key.get(key)
         if hit is not None:
@@ -92,13 +100,13 @@ def phenotype(t: GuidelineTable, haplotypes: tuple[str, ...]) -> PhenotypeResult
     funcs: dict[str, int] = {}
     for h in haplotypes:
         core, n = _split(h)
-        f = "No function" if core in deletion else t.function.get(core)
+        f = "no function" if core in deletion else t.function.get(core)
         if f is None:
             return PhenotypeResult("Indeterminate", None, ["phenotype_not_in_guideline"])
         if n > 1:
             return PhenotypeResult("Indeterminate", None, ["phenotype_not_in_guideline"])  # no CPIC rule for xN here
         funcs[f] = funcs.get(f, 0) + 1
-    hit = t.by_key.get(json.dumps(funcs, sort_keys=True))
+    hit = t.by_key.get(_canon(funcs))
     if hit is None:
         return PhenotypeResult("Indeterminate", None, ["phenotype_not_in_guideline"])
     return PhenotypeResult(hit[0], None)

src/simurg/pgx/star.py

@@ -40,6 +40,9 @@ class GeneModel:
     freq: F64  # (H,) population frequency of each haplotype (prior)
     deletion_allele: str | None = None  # named whole-gene deletion (CYP2D6 *5), or "*0"
     reference: str = ""
+    basis: str = "star_allele"  # PGX-002 result basis: star_allele | haplotype_named | variant_based | tagging_variant
+    activity: dict[str, float] = field(default_factory=dict)  # allele → activity value (variant-based ordering)
+    function: dict[str, str] = field(default_factory=dict)
 
 
 @dataclass

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]] = {}

src/simurg/pipeline/pgx_pheno.py

@@ -97,6 +97,12 @@ def run_phenotype(o: PhenoOptions, rep: Reporter) -> int:
                 res = None
             else:
                 res, per = gene_phenotype(t, feas)
+            if g.get("call_type") == "tagging_variant":  # PGX-007: the CPIC result is the tag status itself
+                from simurg.pgx.phenotype import PhenotypeResult
+
+                res, per = PhenotypeResult(g["genotype"], None), []
+            elif g.get("call_type") == "copy_number":
+                res, per = None, []
             flags = list(g.get("flags", [])) + (res.flags if res else [])
             top = g["candidates"][0] if g["candidates"] else {}
             top_alleles = top.get("alleles", [])

src/simurg/pipeline/pgx_star.py

@@ -7,6 +7,7 @@ and tagging genes, sub-allele collapsing and hybrid CYP2D6 structures follow (K-
 
 from __future__ import annotations
 
+import contextlib
 import gzip
 import json
 import math
@@ -27,6 +28,9 @@ from simurg.pgx.star import CALLED, MASKED, NO_CALL, NOT_COVERED, GeneModel, Obs
 MAX_EDIT_DISTANCE = 3  # PGX-019
 
 CN_REGION = {"CYP2D6": "CYP2D6.p5"}  # PGX-013: total gene copies from the 5′ region
+CN_ONLY_GENES = {"GSTM1": "GSTM1", "GSTT1": "GSTT1", "UGT2B17": "UGT2B17"}  # PGX-006: *0 from whole-gene CN
+ALL_BASES = ("star_allele", "haplotype_named", "variant_based", "tagging_variant")
+TAG_CAVEAT = "tag reliability varies by population"
 DELETION_ALLELE = {"CYP2D6": "*5"}
 
 
@@ -45,11 +49,13 @@ class StarOptions:
     force: bool = False
 
 
-def load_models(db: Path, population: str, genes: list[str] | None = None) -> dict[str, GeneModel]:
+def load_models(
+    db: Path, population: str, genes: list[str] | None = None, bases: tuple[str, ...] = ("star_allele",)
+) -> dict[str, GeneModel]:
     con = sqlite3.connect(db)
-    q = "SELECT gene FROM genes WHERE result_basis = 'star_allele' ORDER BY gene"
+    q = f"SELECT gene, result_basis FROM genes WHERE result_basis IN ({','.join('?' * len(bases))}) ORDER BY gene"
     out = {}
-    for (gene,) in con.execute(q).fetchall():
+    for gene, basis in con.execute(q, bases).fetchall():
         if genes and gene not in genes:
             continue
         alle = con.execute("SELECT allele, kind FROM alleles WHERE gene = ? AND defined = 1 AND structural = 0 "
@@ -83,7 +89,14 @@ def load_models(db: Path, population: str, genes: list[str] | None = None) -> di
             if fr.get(n):
                 freq[i] = float(np.mean(fr[n]))
         freq = np.nan_to_num(freq, nan=0.0)
-        out[gene] = GeneModel(gene, keys, ids, names, carries, freq, DELETION_ALLELE.get(gene), ref)
+        fa = {a: (f, v) for a, f, v in con.execute("SELECT allele, function, activity FROM alleles WHERE gene = ?",
+                                                    (gene,))}  # fmt: skip
+        act = {}
+        for a, (_, v) in fa.items():
+            with contextlib.suppress(TypeError, ValueError):
+                act[a] = float(v)
+        out[gene] = GeneModel(gene, keys, ids, names, carries, freq, DELETION_ALLELE.get(gene), ref, basis, act,
+                              {a: f for a, (f, _) in fa.items() if f})  # fmt: skip
     con.close()
     return out
 
@@ -254,7 +267,7 @@ def run_star_call(o: StarOptions, rep: Reporter) -> int:
         if not Path(str(vp) + ".tbi").exists() and not Path(str(vp) + ".csi").exists():
             raise SimurgError("E_PGX_VCF_UNINDEXED", f"{vp} is not indexed")
     check_outputs_absent([o.out / "diplotype_summary.tsv"], force=o.force)
-    models = load_models(o.db, o.population, o.genes or None)
+    models = load_models(o.db, o.population, o.genes or None, ALL_BASES)
     mcon = sqlite3.connect(o.db)
     masks: dict[str, set[str]] = defaultdict(set)
     for man, probe in mcon.execute("SELECT manifest, probe_name FROM masks"):
@@ -296,6 +309,13 @@ def run_star_call(o: StarOptions, rep: Reporter) -> int:
             unres_rows += g_unres
             for x in g_ncov:
                 not_cov_rows[(gene, str(x["hgvs"]))] = x
+        if not o.no_cn:
+            for gene, region in CN_ONLY_GENES.items():
+                if o.genes and gene not in o.genes:
+                    continue
+                gj, row = cn_only_gene(sid, sample_name, gene, cn.get(region))
+                genes_out.append(gj)
+                rows.append(row)
         calls_json.append({"sample_id": sid, "genes": genes_out})
     cols = ["sample_id", "sample_name", "solution_id", "gene", "call_type", "solution", "array_equivalent_alleles_1",
             "array_equivalent_alleles_2", "no_call_variants", "alleles_impacted_by_no_call", "background_alleles",
@@ -323,6 +343,8 @@ def call_one_gene(
 ) -> tuple[dict[str, Any] | None, list[dict[str, Any]], list[dict[str, Any]], list[dict[str, Any]],
            list[dict[str, Any]]]:  # fmt: skip
     """One gene of one sample: call, details, and the rows of the four TSVs (PGX-012 … PGX-031)."""
+    if m.basis in ("variant_based", "tagging_variant"):
+        return call_variant_gene(sid, sample_name, m, recs, masked, o)
     gene = m.gene
     obs, ploidy = observations(m, recs, masked)
     vrec = variant_records(m, recs)
@@ -367,7 +389,7 @@ def call_one_gene(
     for rank, (s_, d) in enumerate(zip(sols, details, strict=True), start=1):
         sflags = flags + d["flags"] + (["feasible"] if s_ in feasible else [])
         rows.append({"sample_id": sid, "sample_name": sample_name, "solution_id": rank, "gene": gene,
-                     "call_type": "star_allele", "solution": s_.genotype,
+                     "call_type": m.basis, "solution": s_.genotype,
                      "array_equivalent_alleles_1": ";".join(d["alleles"][0]["array_equivalent_alleles"]),
                      "array_equivalent_alleles_2": ";".join(d["alleles"][-1]["array_equivalent_alleles"]),
                      "no_call_variants": ";".join(str(x["hgvs"]) for x in no_call_list),
@@ -398,13 +420,122 @@ def call_one_gene(
         cj.update({"alleles": d["alleles"], "variant_change": d["variant_change"], "flags": d["flags"],
                    "variant_no_calls": no_call_list, "alleles_impacted_by_no_calls": imp_nc})  # fmt: skip
         cands.append(cj)
-    gene_json = {"gene": gene, "call_type": "star_allele", "genotype": ";".join(s_.genotype for s_ in feasible),
+    gene_json = {"gene": gene, "call_type": m.basis, "genotype": ";".join(s_.genotype for s_ in feasible),
                  "feasible_solutions": [s_.genotype for s_ in feasible], "candidates": cands, "copy_number": gcn,
                  "flags": flags, "variant_no_calls": no_call_list, "variants_not_covered": not_cov,
                  "unexplained_variants": unexplained, "alleles_tested": tested}  # fmt: skip
     return gene_json, rows, sup_rows, unres_rows, ncov_rows
 
 
+def call_variant_gene(
+    sid: str, sample_name: str, m: GeneModel, recs: dict[tuple[str, int], list[dict[str, Any]]], masked: set[str],
+    o: StarOptions,
+) -> tuple[dict[str, Any] | None, list[dict[str, Any]], list[dict[str, Any]], list[dict[str, Any]],
+           list[dict[str, Any]]]:  # fmt: skip
+    """PGX-007 / PGX-021: variant-based and tagging genes are reported per variant (``rs9923231 C/T`` or
+    ``Reference``); the diplotype for the phenotype lookup is built from the carried variant alleles — with more than
+    two, the two of lowest activity (CPIC DPYD rule), non-normal function before normal [DB: CPIC DPYD guideline]."""
+    obs, ploidy = observations(m, recs, masked)
+    vrec = variant_records(m, recs)
+    carried: list[str] = []
+    variant_text: list[str] = []
+    sup_rows: list[dict[str, Any]] = []
+    for k, (c, p, r_, a_) in enumerate(m.variants):
+        if obs.state[k] != CALLED or obs.dosage[k] <= 0:
+            continue
+        rec = vrec[k]
+        hg, rs = m.variant_ids[k]
+        gt = "/".join([r_] * (ploidy - int(obs.dosage[k])) + [a_] * int(obs.dosage[k]))
+        variant_text.append(f"{rs or hg} {gt}")
+        owners = [m.haplotypes[h] for h in np.flatnonzero(m.carries[:, k] == 1).tolist()]
+        carried += owners[:1] * int(obs.dosage[k])
+        sup_rows.append({"sample_id": sid, "solution_id": 1, "gene": m.gene, "allele": owners[0] if owners else "",
+                         "hgvs": hg, "rsid": rs, "chrom": c, "pos": p, "ref": r_, "alts": ";".join(rec["alts"]) if rec
+                         else a_, "genotype": rec["gt"] if rec else "", "variant_copy_number": int(obs.dosage[k]),
+                         "gs": None if math.isnan(obs.gs[k]) else round(float(obs.gs[k]), 4), "baf": None
+                         if math.isnan(obs.baf[k]) else round(float(obs.baf[k]), 4), "probe_ids": rec["id"] if rec
+                         else "", "low_gs": str(bool(obs.gs[k] < o.low_gs)).lower()})  # fmt: skip
+
+    def severity(a: str) -> tuple[float, int]:
+        f = (m.function.get(a) or "").lower()
+        return (m.activity.get(a, 1.0), 0 if f and "normal" not in f else 1)
+
+    carried.sort(key=severity)
+    ref = m.reference or "Reference"
+    dip = (carried + [ref] * ploidy)[:ploidy]
+    if m.basis == "tagging_variant":
+        genotype = ";".join(f"{a} {'positive' if a in carried else 'negative'}" for a in m.haplotypes)
+    else:
+        genotype = "; ".join(variant_text) if variant_text else "Reference"
+    v_nc = [k for k in range(len(m.variants)) if obs.state[k] in (NO_CALL, MASKED)]
+    v_ncov = [k for k in range(len(m.variants)) if obs.state[k] == NOT_COVERED]
+    flags = ["unresolved_due_to_no_call"] if v_nc else []
+    if m.basis == "tagging_variant":
+        flags.append("tag_caveat")
+
+    def probe_ids(k: int) -> list[str]:
+        rec = vrec[k]
+        return str(rec["id"]).split(",") if rec else []
+
+    no_call_list = [{"hgvs": m.variant_ids[k][0], "rsid": m.variant_ids[k][1], "probe_ids": probe_ids(k),
+                     "alleles_impacted": [m.haplotypes[h] for h in np.flatnonzero(m.carries[:, k] == 1).tolist()]}
+                    for k in v_nc]  # fmt: skip
+    not_cov = [{"hgvs": m.variant_ids[k][0], "rsid": m.variant_ids[k][1],
+                "alleles_impacted": [m.haplotypes[h] for h in np.flatnonzero(m.carries[:, k] == 1).tolist()]}
+               for k in v_ncov]  # fmt: skip
+    diplotype = "/".join(dip)
+    imp_nc = sorted({m.haplotypes[h] for k in v_nc for h in np.flatnonzero(m.carries[:, k] == 1).tolist()})
+    cand = {"genotype": diplotype, "score": 1.0, "raw_score": 1.0, "equivalents": [diplotype], "copy_number": 2,
+            "alleles": [{"allele": a, "supporting_variants": [], "unresolved_variant_sites": [],
+                         "array_equivalent_alleles": [a], "background_allele": None} for a in dip],
+            "variant_change": None, "flags": [], "variant_no_calls": no_call_list,
+            "alleles_impacted_by_no_calls": imp_nc}  # fmt: skip
+    row = {"sample_id": sid, "sample_name": sample_name, "solution_id": 1, "gene": m.gene,
+           "call_type": "variant" if m.basis == "variant_based" else m.basis, "solution": genotype,
+           "array_equivalent_alleles_1": "", "array_equivalent_alleles_2": "",
+           "no_call_variants": ";".join(str(x["hgvs"]) for x in no_call_list),
+           "alleles_impacted_by_no_call": ";".join(imp_nc), "background_alleles": "",
+           "variant_change": "", "score": "1.0000", "raw_score": "1.0000", "copy_number_solution": "",
+           "flags": ";".join([*flags, "feasible"])}  # fmt: skip
+    gene_json = {"gene": m.gene, "call_type": "variant" if m.basis == "variant_based" else m.basis,
+                 "genotype": genotype, "feasible_solutions": [diplotype], "candidates": [cand], "copy_number": None,
+                 "flags": flags, "variant_no_calls": no_call_list, "variants_not_covered": not_cov,
+                 "unexplained_variants": [], "alleles_tested": sorted(m.haplotypes, key=_natural),
+                 "caveat": TAG_CAVEAT if m.basis == "tagging_variant" else None}  # fmt: skip
+    ncov_rows = [{"gene": m.gene, "hgvs": m.variant_ids[k][0], "rsid": m.variant_ids[k][1],
+                  "alleles_impacted": ";".join(m.haplotypes[h] for h in np.flatnonzero(m.carries[:, k] == 1).tolist())}
+                 for k in v_ncov]  # fmt: skip
+    return gene_json, [row], sup_rows, [], ncov_rows
+
+
+def cn_only_gene(sid: str, sample_name: str, gene: str, cn: int | None) -> tuple[dict[str, Any], dict[str, Any]]:
+    """PGX-006: genes with a whole-gene CN region and no named deletion allele — CN 0 → *0/*0, 1 → *1/*0,
+    2 → *1/*1, n ≥ 3 → *1x(n−1)/*1 (N when capped)."""
+    if cn is None:
+        geno = ""
+    elif cn == 0:
+        geno = "*0/*0"
+    elif cn == 1:
+        geno = "*1/*0"
+    elif cn == 2:
+        geno = "*1/*1"
+    else:
+        geno = f"*1x{'N' if cn >= 5 else cn - 1}/*1"
+    flags = [] if geno else ["no_call"]
+    cand = {"genotype": geno, "score": 1.0, "raw_score": 1.0, "equivalents": [geno] if geno else [], "copy_number": cn,
+            "alleles": [], "variant_change": None, "flags": [], "variant_no_calls": [],
+            "alleles_impacted_by_no_calls": []}  # fmt: skip
+    gj = {"gene": gene, "call_type": "copy_number", "genotype": geno, "feasible_solutions": [geno] if geno else [],
+          "candidates": [cand] if geno else [], "copy_number": cn, "flags": flags, "variant_no_calls": [],
+          "variants_not_covered": [], "unexplained_variants": [], "alleles_tested": ["*0", "*1"]}  # fmt: skip
+    row = {"sample_id": sid, "sample_name": sample_name, "solution_id": 1, "gene": gene, "call_type": "copy_number",
+           "solution": geno, "array_equivalent_alleles_1": "", "array_equivalent_alleles_2": "", "no_call_variants": "",
+           "alleles_impacted_by_no_call": "", "background_alleles": "", "variant_change": "", "score": "1.0000",
+           "raw_score": "1.0000", "copy_number_solution": f"{gene}={'' if cn is None else cn}",
+           "flags": ";".join([*flags, "feasible"])}  # fmt: skip
+    return gj, row
+
+
 def _sol_json(s: Solution) -> dict[str, Any]:
     return {"genotype": s.genotype, "score": round(s.score, 4), "raw_score": round(s.raw_score, 4),
             "equivalents": ["/".join(m) for m in s.members][:20], "copy_number": s.copy_number}  # fmt: skip

tests/unit/test_pgx_star.py

@@ -2,6 +2,8 @@
 
 from __future__ import annotations
 
+from pathlib import Path
+
 import numpy as np
 
 from simurg.pgx.star import CALLED, MASKED, NOT_COVERED, GeneModel, Observations, call_gene, gs_error
@@ -66,3 +68,32 @@ def test_masked_record() -> None:  # PGX-010
     obs, ploidy = observations(m, recs, masked={"p1b"})
     assert ploidy == 2 and obs.state[0] == MASKED and obs.state[1] == CALLED and obs.dosage[1] == 1
     assert obs.state[2] == NOT_COVERED
+
+
+def test_cn_only_genes() -> None:  # PGX-006
+    from simurg.pipeline.pgx_star import cn_only_gene
+
+    assert [cn_only_gene("s", "s", "GSTM1", c)[0]["genotype"] for c in (0, 1, 2, 3, 5, None)] == [
+        "*0/*0", "*1/*0", "*1/*1", "*1x2/*1", "*1xN/*1", ""]  # fmt: skip
+
+
+def test_variant_based_gene() -> None:  # PGX-021: per-variant notation, diplotype by lowest activity
+    from simurg.pipeline.pgx_star import StarOptions, call_variant_gene
+
+    variants = [("1", 100, "C", "T"), ("1", 200, "G", "A"), ("1", 300, "A", "G")]
+    carries = np.array([[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1]], dtype=np.int8)
+    m = GeneModel("DPYD", variants, [("h0", "rs1"), ("h1", "rs2"), ("h2", "rs3")],
+                  ["Reference", "c.A", "c.B", "c.C"], carries, np.array([0.9, 0.05, 0.03, 0.02]), None, "Reference",
+                  "variant_based", {"Reference": 1.0, "c.A": 1.0, "c.B": 0.0, "c.C": 0.5},
+                  {"Reference": "Normal function", "c.A": "Normal function", "c.B": "No function",
+                   "c.C": "Decreased function"})  # fmt: skip
+    recs = {("1", 100): [{"id": "p1", "ref": "C", "alts": ["T"], "gt": "0/1", "gs": "0.5", "baf": "0.5"}],
+            ("1", 200): [{"id": "p2", "ref": "G", "alts": ["A"], "gt": "0/1", "gs": "0.5", "baf": "0.5"}],
+            ("1", 300): [{"id": "p3", "ref": "A", "alts": ["G"], "gt": "0/1", "gs": "0.5", "baf": "0.5"}]}  # fmt: skip
+    gj, rows, sup, _, _ = call_variant_gene("s", "s", m, recs, set(), StarOptions([], Path("x"), Path("y")))
+    assert gj is not None and gj["genotype"] == "rs1 C/T; rs2 G/A; rs3 A/G"
+    assert gj["feasible_solutions"] == ["c.B/c.C"]  # three carried: the two of lowest activity
+    assert rows[0]["call_type"] == "variant" and len(sup) == 3
+    recs = {("1", 100): [{"id": "p1", "ref": "C", "alts": ["T"], "gt": "0/0", "gs": "0.5", "baf": "0"}]}
+    gj, *_ = call_variant_gene("s", "s", m, recs, set(), StarOptions([], Path("x"), Path("y")))
+    assert gj is not None and gj["genotype"] == "Reference" and gj["feasible_solutions"] == ["Reference/Reference"]