From 46422b8de4a9eaf31e987577d01a371e3c377c77 Mon Sep 17 00:00:00 2001 From: Kemal Yaylali Date: Tue, 29 Sep 2026 16:29:54 +0100 Subject: [PATCH] =?UTF-8?q?PGx:=20CYP2D6=20structures=20from=20sub-region?= =?UTF-8?q?=20CN=20=E2=80=94=20*68+*4,=20*36+*10,=20*13=20(PGX-013/018)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Sub-region CNs (p5, intron 2, exon 9) resolve full copies and *36/*68/*13 units; candidates carry the tandem names, background alleles, and CPIC activity sums. GeT-RM: CYP2D6 13/15 top (NA12878 *3/*68+*4 and NA11832 *1/*68+*4 exact), 15/15 within feasible; all genes 132/144. --- docs/validation/V2_pgx_star.md | 17 ++++-- docs/validation/V2_pgx_star_getrm16.json | 15 ++--- scripts/pgx/validate_star_getrm.py | 2 +- src/simurg/pgx/phenotype.py | 15 +++-- src/simurg/pgx/star.py | 73 +++++++++++++++++++++++- src/simurg/pipeline/pgx_star.py | 23 ++++++-- tests/unit/test_pgx_star.py | 15 +++++ 7 files changed, 129 insertions(+), 31 deletions(-) diff --git a/docs/validation/V2_pgx_star.md b/docs/validation/V2_pgx_star.md index f6aff98..6fef22a 100644 --- a/docs/validation/V2_pgx_star.md +++ b/docs/validation/V2_pgx_star.md @@ -30,7 +30,7 @@ Status: **implemented for CPIC star-allele genes, interim validation** (2026-09- |---|---|---| | CYP2C9 | **15/15** | 15/15 | | CYP2C19 | **14/15** | 15/15 | -| CYP2D6 | **12/15** | 14/15 | +| CYP2D6 | **13/15** | **15/15** | | CYP3A5 | 15/15 | 15/15 | | TPMT | 15/16 | 15/16 | | CYP2B6 | 12/14 | 12/14 | @@ -41,16 +41,21 @@ Status: **implemented for CPIC star-allele genes, interim validation** (2026-09- | 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 %)** | +| **total** | **132/144 (91.7 %)** | **135/144 (93.8 %)** | -PGX-032 target (≥ 90 % for CYP2D6/CYP2C19/CYP2C9 against truth): **41/45 = 91 %** (top solution). +PGX-032 target (≥ 90 % for CYP2D6/CYP2C19/CYP2C9 against truth): **42/45 = 93 %** (top solution), 45/45 within +the array-equivalent feasible solutions. + +CYP2D6 structures (PGX-013) come from the three sub-region CNs (`simurg.pgx.star.cyp2d6_structure`): NA12878 +`*3/*68+*4` and NA11832 `*1/*68+*4` are called exactly as GeT-RM reports them; one sample is called `*36+*10/*4` +(p5 = intron 2 = 3, exon 9 = 2; no GeT-RM truth). The two remaining top-solution misses are *1 vs *2, which this +array cannot separate (both within the feasible set). Details and every discordant call: `V2_pgx_star_getrm16.json`. Not in the total (GeT-RM uses pre-2023 names the CPIC definitions no longer carry): UGT1A1 (GeT-RM `*60`; 0/9 top) and SLCO1B1 (`*1B/*17/*21` vs CPIC `*37/*15/*20`; 9/15 top). -Remaining discordances: CYP2D6 `*3/(*68)+*4` (NA12878; *68 hybrids not modelled, *3 frameshift probe), CYP2D6 *2 vs -*1 where the array cannot separate them (prior picks *1), CYP2C19 *8, TPMT *21, CYP2B6 *27, CYP4F2 *2 — rare +Remaining discordances: CYP2D6 *2 vs *1 where the array cannot separate them (prior picks *1), CYP2C19 *8, TPMT *21, CYP2B6 *27, CYP4F2 *2 — rare alleles whose defining variants are absent or masked on this array. ## Error model (K-029) @@ -79,5 +84,5 @@ so the per-call error is 1e-3 (GS ≥ 0.15) and errors are handled by the probe `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 +* Not yet: PharmVar-only genes (CYP2A6, CYP3A4, CYP1A2, …); DPWG tables. diff --git a/docs/validation/V2_pgx_star_getrm16.json b/docs/validation/V2_pgx_star_getrm16.json index 12ac5b4..e77370d 100644 --- a/docs/validation/V2_pgx_star_getrm16.json +++ b/docs/validation/V2_pgx_star_getrm16.json @@ -19,8 +19,8 @@ }, "CYP2D6": { "n": 15, - "top": 12, - "equivalent": 14 + "top": 13, + "equivalent": 15 }, "CYP3A5": { "n": 15, @@ -84,8 +84,8 @@ ], "total": { "n": 144, - "top": 131, - "equivalent": 134 + "top": 132, + "equivalent": 135 }, "discordant": [ { @@ -137,13 +137,6 @@ "call": "rs9923231 T/T", "within_equivalents": false }, - { - "individual": "NA12878", - "gene": "CYP2D6", - "truth": "*3/(*68)+*4", - "call": "*4x2/*1;*4x2/*82", - "within_equivalents": false - }, { "individual": "NA12878", "gene": "UGT1A1", diff --git a/scripts/pgx/validate_star_getrm.py b/scripts/pgx/validate_star_getrm.py index f8a627b..39bbeee 100644 --- a/scripts/pgx/validate_star_getrm.py +++ b/scripts/pgx/validate_star_getrm.py @@ -67,7 +67,7 @@ def truth_alternatives(text: str) -> list[list[set[str]]]: def matches(call: tuple[str, str], truth: list[list[set[str]]]) -> bool: - a, b = (re.sub(r"x\d+$", "", x) for x in call) + a, b = (re.sub(r"x\d+$", "", x.split("+")[-1]) for x in call) # tandem "*68+*4" → core "*4" return any((a in t[0] and b in t[1]) or (a in t[1] and b in t[0]) for t in truth) diff --git a/src/simurg/pgx/phenotype.py b/src/simurg/pgx/phenotype.py index 1e91c30..41273c2 100644 --- a/src/simurg/pgx/phenotype.py +++ b/src/simurg/pgx/phenotype.py @@ -78,12 +78,15 @@ def phenotype(t: GuidelineTable, haplotypes: tuple[str, ...]) -> PhenotypeResult deletion = {"*5", "*0"} if t.activity_based: vals: list[float] = [] - for h in haplotypes: - core, n = _split(h) - a = 0.0 if core in deletion else t.activity.get(core) - if a is None: - return PhenotypeResult("Indeterminate", None, ["phenotype_not_in_guideline"]) - vals.append(a * n) + for h in haplotypes: # a tandem "*36+*10" or "*68+*4" sums its units (CPIC activity-score rule) + hap_act = 0.0 + for unit in h.split("+"): + core, n = _split(unit) + a = 0.0 if core in deletion else t.activity.get(core) + if a is None: + return PhenotypeResult("Indeterminate", None, ["phenotype_not_in_guideline"]) + hap_act += a * n + vals.append(hap_act) counts: dict[str, int] = {} for v in vals: s = f"{v:.1f}" if abs(v * 10 - round(v * 10)) < 1e-9 else f"{v:g}" diff --git a/src/simurg/pgx/star.py b/src/simurg/pgx/star.py index 2ccb9f1..ba4769c 100644 --- a/src/simurg/pgx/star.py +++ b/src/simurg/pgx/star.py @@ -115,14 +115,81 @@ def candidates(m: GeneModel, cn: int | None) -> list[tuple[tuple[int, ...], tupl return out +@dataclass(frozen=True) +class Cyp2d6Structure: + """CYP2D6 copies from the sub-region CNs (PGX-013) [OURS; DOC: PharmVar CYP2D6 structural variation]: + ``full`` complete CYP2D6 copies, ``u36`` *36-type units (CYP2D7 exon 9), ``u68`` *68-type units (CYP2D6 exon 1 + only), ``u13`` *13-type units (CYP2D7 5′ part). CN(p5) = full + u36 + u68, CN(intron 2) = full + u36 + u13, + CN(exon 9) = full + u13; with u68 · u13 = 0 this has one solution.""" + + full: int + u36: int + u68: int + u13: int + + +def cyp2d6_structure(p5: int | None, intron2: int | None, exon9: int | None) -> Cyp2d6Structure | None: + if p5 is None: + return None + ni = p5 if intron2 is None else intron2 + n9 = ni if exon9 is None else exon9 + u36 = max(0, ni - n9) + x = p5 - u36 - n9 + if x >= 0: + full, u68, u13 = n9, x, 0 + else: + full, u68, u13 = p5 - u36, 0, -x + if full < 0: + return Cyp2d6Structure(max(0, p5), 0, 0, 0) + return Cyp2d6Structure(full, u36, u68, u13) + + +def structure_candidates(m: GeneModel, st: Cyp2d6Structure) -> list[tuple[tuple[int, ...], tuple[str, ...]]]: + """Candidates for a CYP2D6 structure: SNV copies = full + u36 (a *36 unit carries the *10 sequence up to exon 8); + *36 units attach to a *10 copy (``*36+*10``), *68 units to *4 when present (``*68+*4``), *13 units are separate + haplotypes.""" + snv_copies = st.full + st.u36 + base = candidates(m, snv_copies) + i10 = m.haplotypes.index("*10") if "*10" in m.haplotypes else None + out = [] + for ix, names in base: + names_l = list(names) + if st.u36: + if i10 is None or ix.count(i10) < st.u36 + 1: # the *36 unit(s) plus the *10 they precede + continue + k = next((n for n, x in enumerate(names_l) if x.startswith("*10")), None) + if k is None: + continue + unit = "*36" if st.u36 == 1 else f"*36x{st.u36}" + names_l[k] = f"{unit}+*10" + if st.u68: + unit = "*68" if st.u68 == 1 else f"*68x{st.u68}" + k = next((n for n, x in enumerate(names_l) if x.split("+")[-1].split("x")[0] == "*4"), 0) + names_l[k] = f"{unit}+{names_l[k]}" + if st.u13: + if len(names_l) == 2 and names_l[1] == (m.deletion_allele or "*0"): + names_l[1] = "*13" if st.u13 == 1 else f"*13x{st.u13}" + else: + names_l[0] = f"*13+{names_l[0]}" + out.append((ix, tuple(names_l))) + return out + + def call_gene( m: GeneModel, obs: Observations, cn: int | None = None, baf_sigma: float = 0.06, ploidy: int = 2, - max_candidates: int = 10, feasible_ratio: float = 0.5, + max_candidates: int = 10, feasible_ratio: float = 0.5, structure: Cyp2d6Structure | None = None, ) -> tuple[list[Solution], list[Solution]]: # fmt: skip """All scored solutions (ranked) and the feasible ones (PGX-014/015).""" called = obs.state == CALLED - haploid = [((i,), (m.haplotypes[i],)) for i in range(len(m.haplotypes))] - cands = haploid if ploidy == 1 else candidates(m, cn) + haploid: list[tuple[tuple[int, ...], tuple[str, ...]]] = [((i,), (m.haplotypes[i],)) for i in + range(len(m.haplotypes))] # fmt: skip + if ploidy == 1: + cands = haploid + elif structure is not None and (structure.u36 or structure.u68 or structure.u13): + cands = structure_candidates(m, structure) + cn = structure.full + structure.u36 + else: + cands = candidates(m, cn) if not cands: return [], [] carries = m.carries.astype(np.int16) diff --git a/src/simurg/pipeline/pgx_star.py b/src/simurg/pipeline/pgx_star.py index f28e530..96999a4 100644 --- a/src/simurg/pipeline/pgx_star.py +++ b/src/simurg/pipeline/pgx_star.py @@ -23,11 +23,22 @@ import numpy.typing as npt from simurg import __version__ from simurg.core.errors import ExitCode, SimurgError from simurg.core.log import Reporter -from simurg.pgx.star import CALLED, MASKED, NO_CALL, NOT_COVERED, GeneModel, Observations, Solution, call_gene +from simurg.pgx.star import ( + CALLED, + MASKED, + NO_CALL, + NOT_COVERED, + GeneModel, + Observations, + Solution, + call_gene, + cyp2d6_structure, +) 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_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" @@ -350,8 +361,11 @@ def call_one_gene( vrec = variant_records(m, recs) region = CN_REGION.get(gene) gcn = cn.get(region) if region else None - sols, feasible = call_gene(m, obs, cn=gcn if region else None, ploidy=ploidy, - max_candidates=o.max_candidates, feasible_ratio=o.feasible_ratio) # fmt: skip + structure = None + if gene == "CYP2D6" and gcn is not None: + structure = cyp2d6_structure(gcn, cn.get("CYP2D6.intron.2"), cn.get("CYP2D6.exon.9")) + sols, feasible = call_gene(m, obs, cn=gcn if region else None, ploidy=ploidy, max_candidates=o.max_candidates, + feasible_ratio=o.feasible_ratio, structure=structure) # fmt: skip if not sols: return None, [], [], [], [] top = sols[0] @@ -398,7 +412,8 @@ def call_one_gene( if al["background_allele"]), "variant_change": ";".join(d["variant_change"] or []), "score": f"{s_.score:.4f}", "raw_score": f"{s_.raw_score:.4f}", - "copy_number_solution": f"{region}={gcn}" if region and gcn is not None else "", + "copy_number_solution": ";".join(f"{r}={cn.get(r)}" for r in CN_SOLUTION_REGIONS.get(gene, ()) + if cn.get(r) is not None), "flags": ";".join(sflags)}) # fmt: skip for al in d["alleles"]: for sv in al["supporting_variants"]: diff --git a/tests/unit/test_pgx_star.py b/tests/unit/test_pgx_star.py index d41ec8a..74a9992 100644 --- a/tests/unit/test_pgx_star.py +++ b/tests/unit/test_pgx_star.py @@ -97,3 +97,18 @@ def test_variant_based_gene() -> None: # PGX-021: per-variant notation, diploty 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"] + + +def test_cyp2d6_structures() -> None: # PGX-013 CN rules + from simurg.pgx.star import Cyp2d6Structure, cyp2d6_structure, structure_candidates + + assert cyp2d6_structure(3, 2, 2) == Cyp2d6Structure(2, 0, 1, 0) # *68 unit (NA12878 *3/*68+*4) + assert cyp2d6_structure(3, 3, 2) == Cyp2d6Structure(2, 1, 0, 0) # *36 unit + assert cyp2d6_structure(1, 2, 2) == Cyp2d6Structure(1, 0, 0, 1) # *13 hybrid + assert cyp2d6_structure(2, 2, 2) == Cyp2d6Structure(2, 0, 0, 0) + m = _model() + m.haplotypes = ["*1", "*4", "*10", "*2", "*5"] + names = {n for _, n in structure_candidates(m, Cyp2d6Structure(2, 0, 1, 0))} + assert ("*1", "*68+*4") in names and ("*68+*1", "*2") in names + names = {n for _, n in structure_candidates(m, Cyp2d6Structure(2, 1, 0, 0))} + assert ("*36+*10", "*4") in names and all("*36+*10" in n for n in names)