From ed5e92b8593b82d2e61531bfd7c18570d8331c9c Mon Sep 17 00:00:00 2001 From: Kemal Yaylali Date: Tue, 29 Sep 2026 17:01:36 +0100 Subject: [PATCH] PGx: ClinPGx variant-based genes (PGX-008) and --guideline dpwg (PGX-025) - fetch_sources.sh downloads ClinPGx clinicalVariants (CC BY-SA 4.0); pgx db build adds genotype-only variant-based genes (levels 1A-2B, level 3 for the parity list), rsIDs resolved from the manifest's probes: 70 genes in the DB. - pgx phenotype --guideline dpwg: CPIC mappings for CPIC genes (as ClinPGx curates DPWG), no phenotype for DPWG-only genes (flags dpwg_via_cpic_mapping / dpwg_mapping_not_available). K-029 updated. --- THIRD_PARTY_NOTICES.md | 1 + docs/known_unknowns.md | 2 +- docs/validation/V2_pgx_star.md | 10 ++++++- scripts/pgx/fetch_sources.sh | 17 +++++++++++ src/simurg/pgx/db.py | 29 ++++++++++++++++++ src/simurg/pipeline/pgx_db.py | 51 ++++++++++++++++++++++++++++++-- src/simurg/pipeline/pgx_pheno.py | 18 +++++++---- tests/unit/test_pgx_db.py | 43 +++++++++++++++++++++++++++ 8 files changed, 161 insertions(+), 10 deletions(-) diff --git a/THIRD_PARTY_NOTICES.md b/THIRD_PARTY_NOTICES.md index 2410019..69102e0 100644 --- a/THIRD_PARTY_NOTICES.md +++ b/THIRD_PARTY_NOTICES.md @@ -18,6 +18,7 @@ used only as external subprocess oracles and are not copied or linked. | [CPIC](https://cpicpgx.org) database via the CPIC API (allele definitions, allele function, diplotype→phenotype tables) | CC0 1.0 | downloaded by `scripts/pgx/fetch_sources.sh` into `data/` (never committed); `simurg pgx db build` derives the local PGx DB from it; source version and file hashes are recorded in the DB `meta` table | | [PharmVar](https://www.pharmvar.org) haplotype definitions (Gaedigk et al. 2018 CPT 103:399) | PharmVar terms of use (attribution) | per-gene zips downloaded by `scripts/pgx/fetch_sources.sh` into `data/` (never committed); `simurg pgx db build` derives core-allele definitions for genes CPIC does not define | | [GeT-RM](https://www.cdc.gov/lab-quality/php/get-rm/reference-materials.html) consolidated PGx/HLA table (CDC) | US government work | validation truth only (`scripts/pgx/getrm_truth.py`), not redistributed | +| [ClinPGx](https://www.clinpgx.org) (formerly PharmGKB) clinical-annotation variants | CC BY-SA 4.0 (attribution: ClinPGx, https://www.clinpgx.org/page/citingClinpgx) | `clinicalVariants.zip` downloaded by `scripts/pgx/fetch_sources.sh` into `data/` (never committed); rsIDs become genotype-only variant-based genes in the local PGx DB | ## BSD-2-Clause (BeadArrayFiles, verbatim from its README) diff --git a/docs/known_unknowns.md b/docs/known_unknowns.md index 6c30ee3..4389c57 100644 --- a/docs/known_unknowns.md +++ b/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. | Open: hybrids, G6PD/variant-based genes, PharmVar genes, phenotypes; larger/independent truth (`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. | Open: larger/independent truth set; 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. | diff --git a/docs/validation/V2_pgx_star.md b/docs/validation/V2_pgx_star.md index 2164381..ff3e449 100644 --- a/docs/validation/V2_pgx_star.md +++ b/docs/validation/V2_pgx_star.md @@ -92,4 +92,12 @@ so the per-call error is 1e-3 (GS ≥ 0.15) and errors are handled by the probe from the GRCh38 haplotype VCFs, `*1` = reference, prior reference 0.8 (no population frequencies); genotype only (no CPIC phenotype). The probe mask was re-derived with these variants (one more record: CYP2A6 19:40843845 C>A, concordance 0.47, which had produced false *5 calls). -* Not yet: DPWG tables. +* ClinPGx clinical-annotation variants (CC BY-SA 4.0; levels 1A–2B, and level 3 for the PGX-008 parity genes) add + genotype-only variant-based genes — OPRM1, COMT, DRD2, ANKK1, MTHFR, F5, … (43 genes, 171 of 187 rsIDs resolved on + GSA v4 PGx from the manifest's probes). The DB now holds 70 genes (16 star-allele, G6PD, 2 HLA tag, 51 + variant-based; spec parity list: 58). `pgx coverage`: 1225 alleles callable, 93 array-equivalent, 437 not + covered, 5 masked. +* `--guideline dpwg` (PGX-025): ClinPGx curates DPWG guidelines with the CPIC allele-function and phenotype mappings + 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). diff --git a/scripts/pgx/fetch_sources.sh b/scripts/pgx/fetch_sources.sh index be0f2a4..4ca1f8f 100644 --- a/scripts/pgx/fetch_sources.sh +++ b/scripts/pgx/fetch_sources.sh @@ -56,3 +56,20 @@ files = {p.name: hashlib.sha256(p.read_bytes()).hexdigest() for p in sorted(d.gl "files": files}, indent=1)) print("PharmVar version", ",".join(sorted(vers))) PY + +# ---- ClinPGx (formerly PharmGKB; CC BY-SA 4.0): clinical-annotation variants (variant-based genes, genotype only) +CP="$OUT/clinpgx-$DATE" +mkdir -p "$CP" +curl -sSfL --retry 3 -o "$CP/clinicalVariants.zip" "https://api.clinpgx.org/v1/download/file/data/clinicalVariants.zip" +python3 - "$CP" <<'PY' +import hashlib, json, pathlib, sys, zipfile +d = pathlib.Path(sys.argv[1]) +z = zipfile.ZipFile(d / "clinicalVariants.zip") +created = next((n.split("_", 1)[1].rsplit(".", 1)[0] for n in z.namelist() if n.startswith("CREATED_")), "") +(d / "source.json").write_text(json.dumps({"name": "ClinPGx clinical variants", "version": created, + "retrieved": d.name.split("-")[1], "url": "https://www.clinpgx.org", + "license": "CC BY-SA 4.0", + "files": {"clinicalVariants.zip": hashlib.sha256( + (d / "clinicalVariants.zip").read_bytes()).hexdigest()}}, indent=1)) +print("ClinPGx clinical variants", created) +PY diff --git a/src/simurg/pgx/db.py b/src/simurg/pgx/db.py index e32076e..6f0b6b7 100644 --- a/src/simurg/pgx/db.py +++ b/src/simurg/pgx/db.py @@ -258,3 +258,32 @@ def pharmvar_core_alleles(gene: str, vcfs: dict[str, str]) -> dict[str, list[tup recs.append((f[0], int(f[1]), f[3].upper(), alt.upper(), rsid)) out[f"*{m.group(1)}"] = recs return out + + +CLINPGX_LEVELS = ("1A", "1B", "2A", "2B") +# PGX-008 parity list of variant-based genes: level-3 annotations are kept for these too +PARITY_VARIANT_GENES = {"DPYD", "RYR1", "OPRM1", "CFTR", "DRD2", "COMT", "OPRD1", "CACNA1C", "HTR2A", "OPRK1", + "MT-RNR1", "TNF", "BDNF", "IL6", "HTR2C", "GRIN2B", "MTHFR", "ITGB3", "FKBP5", "ADH1B", + "ANKK1", "HMGCR", "ATM", "KIF6", "F2", "CACNA1S", "LPA", "F13A1", "GRIK1", "GRIK4", "YEATS4", + "ABCG2", "F5", "ALDH2", "ANK3", "VKORC1"} # fmt: skip + + +def clinpgx_variant_genes(tsv: str, exclude: set[str], levels: tuple[str, ...] = CLINPGX_LEVELS) -> dict[str, set[str]]: + """ClinPGx ``clinicalVariants.tsv`` → gene → rsIDs with a clinical annotation at the given evidence levels, for + single-gene rows of genes that no allele-definition source covers (level 3 as well for the PGX-008 parity genes) + [DB: ClinPGx clinical annotations, CC BY-SA 4.0]. + Positions and alleles are resolved later from the manifest's probes (genotype-only variant-based genes).""" + out: dict[str, set[str]] = {} + lines = tsv.splitlines() + head = lines[0].split("\t") + iv, ig, il = head.index("variant"), head.index("gene"), head.index("level of evidence") + for ln in lines[1:]: + f = ln.split("\t") + if len(f) <= max(iv, ig, il): + continue + gene, var, lev = f[ig].strip(), f[iv].strip(), f[il].strip() + ok_level = lev in levels or (lev == "3" and gene in PARITY_VARIANT_GENES) + if "," in gene or not gene or gene in exclude or not ok_level or not re.fullmatch(r"rs\d+", var): + continue + out.setdefault(gene, set()).add(var) + return out diff --git a/src/simurg/pipeline/pgx_db.py b/src/simurg/pipeline/pgx_db.py index 644710d..6988471 100644 --- a/src/simurg/pipeline/pgx_db.py +++ b/src/simurg/pipeline/pgx_db.py @@ -254,6 +254,47 @@ def run_db_build(o: PgxDbOptions, rep: Reporter) -> int: 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)) + # ClinPGx clinical-annotation variants of genes no definition source covers: genotype-only variant-based genes + cp_src = None + cdir2 = _latest(o.sources, "clinpgx") + if cdir2 is not None: + import zipfile + + from simurg.pgx.db import clinpgx_variant_genes + + cp_src = json.loads((cdir2 / "source.json").read_text()) + with zipfile.ZipFile(cdir2 / "clinicalVariants.zip") as z: + tsv = z.read("clinicalVariants.tsv").decode() + known = {g for g, *_ in rows["genes"]} | {g for g, *_ in TAGS} + cp_genes = clinpgx_variant_genes(tsv, known) + resolved: dict[tuple[str, str], tuple[Any, ...]] = {} + want = {rs for v in cp_genes.values() for rs in v} + for bpm_path in o.manifests: + name = read_bpm(bpm_path).manifest_name + for ps in _probe_sites(bpm_path, csvs.get(bpm_path), genome, cmap, {}, rep, want): + alts = [a for a in ps.alleles if a != ps.ref] + if len(alts) != 1 or ps.ref not in ps.alleles: + continue + for rs in re.findall(r"rs\d+", ps.name): + if rs in want: + hg = full_hgvs(ps.contig, f"g.{ps.pos}{ps.ref}>{alts[0]}") + resolved.setdefault((rs, ps.contig), (ps.contig, ps.pos, ps.ref, alts[0], hg)) + rows["probe_map"].append((name, ps.name, ps.contig, ps.pos, ps.ref, alts[0])) + by_rs = {rs: v for (rs, _), v in resolved.items()} + lid2 = 950_000_000 + for gene in sorted(cp_genes): + chrom = None + rows["alleles"].append((gene, "Reference", "reference", None, None, 0, 1)) + for rs in sorted(cp_genes[gene], key=lambda r: int(r[2:])): + v = by_rs.get(rs) + rows["alleles"].append((gene, rs, "rsid_variant_allele", None, None, 0, int(v is not None))) + if v is not None: + lid2 += 1 + chrom = v[0] + rows["allele_variants"].append((gene, rs, "GRCh38", v[0], v[1], v[2], v[3], v[4], rs, lid2, v[3])) + rows["genes"].append((gene, "variant_based", chrom)) + rep.info(f"ClinPGx: {len(cp_genes)} variant-based gene(s), {len(by_rs)} of {len(want)} rsID(s) resolved " + f"on the manifest(s)") # fmt: skip 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 @@ -270,7 +311,10 @@ def run_db_build(o: PgxDbOptions, rep: Reporter) -> int: "license": csrc.get("license")}] + ([{"name": "PharmVar", "version": pv_src["version"], "retrieved": pv_src["retrieved"], "license": pv_src.get("license")}] - if pv_src else [])), + if pv_src else []) + + ([{"name": "ClinPGx clinical variants", "version": cp_src["version"], + "retrieved": cp_src["retrieved"], "license": cp_src.get("license")}] + if cp_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 @@ -301,8 +345,9 @@ def run_db_build(o: PgxDbOptions, rep: Reporter) -> int: rep.info(f"identical PGx DB already present: {final}") return ExitCode.OK tmp.replace(final) - rep.info(f"PGx DB: {len(defs)} genes, {len(rows['alleles'])} alleles, {len(rows['allele_variants'])} variant " - f"rows, {len(rows['probe_map'])} probe links → {final} (content sha256 {sha})") # fmt: skip + n_g, n_a, n_v, n_p = (len(rows[k]) for k in ("genes", "alleles", "allele_variants", "probe_map")) + rep.info(f"PGx DB: {n_g} genes, {n_a} alleles, {n_v} variant rows, {n_p} probe links → {final} " + f"(content sha256 {sha})") # fmt: skip return ExitCode.OK diff --git a/src/simurg/pipeline/pgx_pheno.py b/src/simurg/pipeline/pgx_pheno.py index 8e5224b..dd61007 100644 --- a/src/simurg/pipeline/pgx_pheno.py +++ b/src/simurg/pipeline/pgx_pheno.py @@ -17,6 +17,11 @@ from simurg.core.log import Reporter from simurg.pgx.phenotype import build_table, gene_phenotype JSON_SCHEMA = "simurg.pgx/1" +# PGX-025 DPWG gene coverage. ClinPGx curates DPWG guidelines with the CPIC allele-function and phenotype mappings for +# genes CPIC covers [DOC: ClinPGx, "Genotype-phenotype mapping for curation of DPWG guidelines"], so those genes use the +# CPIC tables (flag dpwg_via_cpic_mapping); DPWG-only genes have no public machine-readable mapping here. +DPWG_GENES = ("CYP1A2", "CYP2B6", "CYP2C19", "CYP2C9", "CYP2D6", "CYP3A4", "CYP3A5", "DPYD", "NUDT15", "SLCO1B1", + "TPMT", "UGT1A1", "VKORC1", "F5") # fmt: skip TSV_SCHEMA = "simurg.pgx_tsv/1" @@ -65,9 +70,8 @@ def run_phenotype(o: PhenoOptions, rep: Reporter) -> int: """``simurg pgx phenotype`` (CLI-034).""" from simurg.core.atomic import atomic_write_text - if o.guideline != "cpic": - raise SimurgError("E_PGX_GUIDELINE", f"guideline {o.guideline!r} is not in this PGx DB (CPIC only)", - hint="DPWG tables are not built yet (K-029)") # fmt: skip + if o.guideline not in ("cpic", "dpwg"): + raise SimurgError("E_PGX_GUIDELINE", f"unknown guideline {o.guideline!r}", hint="cpic or dpwg") src = o.calls / "pgx_calls.json" if not src.exists(): raise SimurgError("E_NO_INPUT", f"{src} not found", hint="run pgx star-allele call first") @@ -77,9 +81,11 @@ def run_phenotype(o: PhenoOptions, rep: Reporter) -> int: meta = db_meta(o.db) con = sqlite3.connect(o.db) tables = {} - for (gene,) in con.execute("SELECT DISTINCT gene FROM phenotypes WHERE guideline = ?", (o.guideline,)): + for (gene,) in con.execute("SELECT DISTINCT gene FROM phenotypes WHERE guideline = 'cpic'"): + if o.guideline == "dpwg" and gene not in DPWG_GENES: + continue rows = con.execute("SELECT lookup_key, activity_score, phenotype FROM phenotypes WHERE gene = ? AND " - "guideline = ?", (gene, o.guideline)).fetchall() # fmt: skip + "guideline = 'cpic'", (gene,)).fetchall() # fmt: skip alle = con.execute("SELECT allele, function, activity FROM alleles WHERE gene = ?", (gene,)).fetchall() tables[gene] = build_table(gene, rows, alle) tested = {g: [a for (a,) in con.execute("SELECT allele FROM alleles WHERE gene = ? AND defined = 1", (g,))] @@ -104,6 +110,8 @@ def run_phenotype(o: PhenoOptions, rep: Reporter) -> int: elif g.get("call_type") == "copy_number": res, per = None, [] flags = list(g.get("flags", [])) + (res.flags if res else []) + if o.guideline == "dpwg" and g["gene"] in DPWG_GENES: + flags.append("dpwg_via_cpic_mapping" if t is not None else "dpwg_mapping_not_available") top = g["candidates"][0] if g["candidates"] else {} top_alleles = top.get("alleles", []) entry = {"gene": g["gene"], "call_type": g.get("call_type", "star_allele"), "genotype": g["genotype"], diff --git a/tests/unit/test_pgx_db.py b/tests/unit/test_pgx_db.py index 08b7598..e672f12 100644 --- a/tests/unit/test_pgx_db.py +++ b/tests/unit/test_pgx_db.py @@ -156,3 +156,46 @@ def test_pharmvar_core_alleles() -> None: # PGX-001 PharmVar source core = pharmvar_core_alleles("CYP2A6", {"CYP2A6-6.2/GRCh38/CYP2A6_9.vcf": vcf, "CYP2A6-6.2/GRCh38/CYP2A6_9.001.vcf": vcf}) # fmt: skip assert core == {"*9": [("chr19", 100, "A", "C", "rs1")]} # sub-alleles (*9.001) are not core definitions + + +def test_clinpgx_variant_genes() -> None: # PGX-008 variant-based genes from ClinPGx + from simurg.pgx.db import clinpgx_variant_genes + + tsv = ("variant\tgene\ttype\tlevel of evidence\tchemicals\tphenotypes\n" + "rs1\tOPRM1\tEfficacy\t3\tx\t\nrs2\tFOO\tEfficacy\t3\tx\t\nrs3\tFOO\tEfficacy\t2A\tx\t\n" + "CYP2C9*3\tCYP2C9\tPK\t1A\tx\t\nrs4\tCYP2C9\tPK\t1A\tx\t\nrs5\tA,B\tPK\t1A\tx\t\n") # fmt: skip + assert clinpgx_variant_genes(tsv, {"CYP2C9"}) == {"OPRM1": {"rs1"}, "FOO": {"rs3"}} + + +def test_dpwg_mode(tmp_path: Path) -> None: # PGX-025 + from simurg.core.log import Reporter + from simurg.pipeline.pgx_pheno import PhenoOptions, run_phenotype + + src = tmp_path / "sources" / "cpic-20260101" + src.mkdir(parents=True) + t = _tables() + t["gene"][0]["symbol"] = "CYP2D6" + for tbl in ("sequence_location", "allele_definition", "allele", "gene_result"): + for r in t[tbl]: + r["genesymbol"] = "CYP2D6" + for k, v in t.items(): + (src / f"{k}.json").write_text(json.dumps(v)) + (src / "source.json").write_text(json.dumps({"name": "CPIC", "version": "2026-01-01", "retrieved": "20260101"})) + fa = tmp_path / "g.fa" + fa.write_text(">chr1\n" + SEQ + "\n") + r = CliRunner().invoke(app, ["pgx", "db", "build", "--sources", str(tmp_path / "sources"), "--fasta", str(fa), + "--out", str(tmp_path / "db"), "--quiet"]) # fmt: skip + assert r.exit_code == 0, r.output + db = next((tmp_path / "db").glob("pgx_db-1-*.sqlite")) + calls = tmp_path / "calls" + calls.mkdir() + (calls / "pgx_calls.json").write_text(json.dumps({"samples": [{"sample_id": "s1", "genes": [ + {"gene": "CYP2D6", "genotype": "*2/*2", "feasible_solutions": ["*2/*2"], "candidates": [ + {"genotype": "*2/*2", "score": 1.0, "raw_score": 1.0}]}, + {"gene": "CYP3A4", "genotype": "*1/*22", "feasible_solutions": ["*1/*22"], "candidates": []}]}]})) # fmt: skip + rep = Reporter(command="t", quiet=True) + assert run_phenotype(PhenoOptions(calls, db, guideline="dpwg"), rep) == 0 + doc = json.loads((calls / "s1.pgx.json").read_text()) + by = {g["gene"]: g for g in doc["genes"]} + assert by["CYP2D6"]["phenotype"] == "Poor Metabolizer" and "dpwg_via_cpic_mapping" in by["CYP2D6"]["flags"] + assert by["CYP3A4"]["phenotype"] is None and "dpwg_mapping_not_available" in by["CYP3A4"]["flags"]