PGx: PharmVar genes (CYP1A2, CYP2A6, CYP2A13, CYP2C8, CYP3A4), CYP2A6 *4 from CN (PGX-001/006)

fetch_sources.sh downloads PharmVar per-gene zips; pgx db build adds core alleles for genes CPIC does not
define (normalised against the FASTA, *1 = reference). Masks re-derived with the new variants. GeT-RM:
CYP3A4 15/15, CYP2A6 7/9; all comparable genes 154/168 (CYP1A2/CYP2C8 definitions differ from GeT-RM's).

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

src/simurg/pipeline/pgx_db.py

@@ -168,6 +168,51 @@ def run_db_build(o: PgxDbOptions, rep: Reporter) -> int:
         rows["sites"].append((s.gene, lid, s.chrom, full_hgvs(s.chrom, s.hgvs), s.rsid, s.ref_value))
         for val, (pos, ref, alt) in sorted(s.alleles.items()):
             rows["site_values"].append((s.gene, lid, val, pos, ref, alt))
+    # PharmVar genes that CPIC does not define (PGX-001): core alleles from the GRCh38 haplotype VCFs
+    pv_src = None
+    pdir = _latest(o.sources, "pharmvar")
+    if pdir is not None:
+        import zipfile
+
+        from simurg.pgx.db import normalize, pharmvar_core_alleles
+
+        pv_src = json.loads((pdir / "source.json").read_text())
+        lid = 900_000_000
+        for zpath in sorted(pdir.glob("*.zip")):
+            gene = zpath.stem
+            if gene in defs:
+                continue
+            with zipfile.ZipFile(zpath) as z:
+                vcfs = {n: z.read(n).decode() for n in z.namelist() if "/GRCh38/" in n and n.endswith(".vcf")}
+            core = pharmvar_core_alleles(gene, vcfs)
+            if not core:
+                continue
+            core.setdefault("*1", [])  # the reference haplotype has no core VCF of its own
+            chrom = None
+            n_var = 0
+            for allele in sorted(core, key=lambda a: int(a[1:])):
+                bad = False
+                norm: list[tuple[str, int, str, str, str | None, str]] = []
+                for c, vpos, r, alt_v, rsid in core[allele]:
+                    contig = cmap.get(c)
+                    if contig is None or genome.fetch(contig, vpos - 1, vpos - 1 + len(r)) != r:
+                        bad = True
+                        break
+                    npos, nref, nalt = normalize(genome.fetch, contig, vpos, vpos + len(r) - 1, alt_v)
+                    norm.append((contig, npos, nref, nalt, rsid, f"g.{vpos}{r}>{alt_v}"))
+                    chrom = c
+                if bad:
+                    rep.info(f"definition note: PharmVar {gene} {allele}: REF does not match the FASTA, skipped")
+                    continue
+                kind = "reference" if allele == "*1" else "star_allele"
+                rows["alleles"].append((gene, allele, kind, None, None, 0, int(bool(norm) or allele == "*1")))
+                for contig, pos, ref, alt, rsid, hg in norm:
+                    lid += 1
+                    rows["allele_variants"].append((gene, allele, "GRCh38", contig, pos, ref, alt,
+                                                    full_hgvs(contig, hg), rsid, lid, alt))  # fmt: skip
+                    n_var += 1
+            rows["genes"].append((gene, "star_allele", chrom))
+            rep.info(f"PharmVar {gene}: {len(core)} core allele(s), {n_var} variant row(s)")
     results = {r["id"]: r for r in t["gene_result"]}
     for lk in t["gene_result_lookup"]:
         res = results.get(lk["phenotypeid"])
@@ -222,7 +267,10 @@ def run_db_build(o: PgxDbOptions, rep: Reporter) -> int:
     meta = {"schema": "simurg.pgx_db", "schema_version": str(DB_SCHEMA), "builder": f"SimurgArray {__version__}",
             "genome_build": "GRCh38", "fasta_sha256": genome.sha256, "sqlite_version": sqlite3.sqlite_version,
             "allele_sources": json.dumps([{"name": "CPIC", "version": csrc["version"], "retrieved": csrc["retrieved"],
-                                           "license": csrc.get("license")}]),
+                                           "license": csrc.get("license")}]
+                                         + ([{"name": "PharmVar", "version": pv_src["version"],
+                                              "retrieved": pv_src["retrieved"], "license": pv_src.get("license")}]
+                                            if pv_src else [])),
             "phenotype_sources": json.dumps([{"name": "CPIC", "version": csrc["version"]}]),
             "cpic_files_sha256": json.dumps(csrc.get("files", {}), sort_keys=True),
             "manifests": json.dumps(sorted(manifests)), "definition_notes": json.dumps(problems)}  # fmt: skip

src/simurg/pipeline/pgx_star.py

@@ -37,12 +37,15 @@ from simurg.pgx.star import (
 
 MAX_EDIT_DISTANCE = 3  # PGX-019
 
-CN_REGION = {"CYP2D6": "CYP2D6.p5"}  # PGX-013: total gene copies from the 5′ region
-CN_SOLUTION_REGIONS = {"CYP2D6": ("CYP2D6.p5", "CYP2D6.intron.2", "CYP2D6.exon.9")}  # PGX-027 copy_number_solution
+CN_REGION = {"CYP2D6": "CYP2D6.p5", "CYP2A6": "CYP2A6.exon.1"}  # PGX-013: gene copies from the 5′ region
+CN_SOLUTION_REGIONS = {
+    "CYP2D6": ("CYP2D6.p5", "CYP2D6.intron.2", "CYP2D6.exon.9"),
+    "CYP2A6": ("CYP2A6.exon.1", "CYP2A6.intron.7"),
+}  # PGX-027 copy_number_solution
 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"}
+DELETION_ALLELE = {"CYP2D6": "*5", "CYP2A6": "*4"}  # named whole-gene deletions (PGX-006)
 
 
 @dataclass
@@ -100,6 +103,10 @@ def load_models(
             if fr.get(n):
                 freq[i] = float(np.mean(fr[n]))
         freq = np.nan_to_num(freq, nan=0.0)
+        if not freq.any() and names:  # no population frequencies (PharmVar genes): reference 0.8, others share 0.2
+            freq = np.full(len(names), 0.2 / max(1, len(names) - 1))
+            if ref in names:
+                freq[names.index(ref)] = 0.8
         fa = {a: (f, v) for a, f, v in con.execute("SELECT allele, function, activity FROM alleles WHERE gene = ?",
                                                     (gene,))}  # fmt: skip
         act = {}