From 9a02f39d01991447bb32375765f2865d71ad04e4 Mon Sep 17 00:00:00 2001 From: Kemal Yaylali Date: Tue, 29 Sep 2026 16:24:33 +0100 Subject: [PATCH] =?UTF-8?q?PGx:=20all=20CPIC=20gene=20kinds=20=E2=80=94=20?= =?UTF-8?q?G6PD,=20variant-based,=20HLA=20tags,=20*0=20CN=20genes=20(PGX-0?= =?UTF-8?q?06/007/021)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - 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). --- docs/validation/V2_pgx_star.md | 13 +- docs/validation/V2_pgx_star_getrm16.json | 54 ++++++++- scripts/pgx/getrm_truth.py | 4 + scripts/pgx/validate_star_getrm.py | 38 +++++- src/simurg/pgx/phenotype.py | 18 ++- src/simurg/pgx/star.py | 3 + src/simurg/pipeline/pgx_db.py | 47 +++++++- src/simurg/pipeline/pgx_pheno.py | 6 + src/simurg/pipeline/pgx_star.py | 145 +++++++++++++++++++++-- tests/unit/test_pgx_star.py | 31 +++++ 10 files changed, 336 insertions(+), 23 deletions(-) diff --git a/docs/validation/V2_pgx_star.md b/docs/validation/V2_pgx_star.md index 234f8e9..f6aff98 100644 --- a/docs/validation/V2_pgx_star.md +++ b/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. diff --git a/docs/validation/V2_pgx_star_getrm16.json b/docs/validation/V2_pgx_star_getrm16.json index 00899e5..12ac5b4 100644 --- a/docs/validation/V2_pgx_star_getrm16.json +++ b/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", diff --git a/scripts/pgx/getrm_truth.py b/scripts/pgx/getrm_truth.py index 74c5e10..e5cd35f 100644 --- a/scripts/pgx/getrm_truth.py +++ b/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: diff --git a/scripts/pgx/validate_star_getrm.py b/scripts/pgx/validate_star_getrm.py index 8821719..f8a627b 100644 --- a/scripts/pgx/validate_star_getrm.py +++ b/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 diff --git a/src/simurg/pgx/phenotype.py b/src/simurg/pgx/phenotype.py index e16dd69..1e91c30 100644 --- a/src/simurg/pgx/phenotype.py +++ b/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) diff --git a/src/simurg/pgx/star.py b/src/simurg/pgx/star.py index dfe851f..2ccb9f1 100644 --- a/src/simurg/pgx/star.py +++ b/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 diff --git a/src/simurg/pipeline/pgx_db.py b/src/simurg/pipeline/pgx_db.py index 2acf9a4..acb3d72 100644 --- a/src/simurg/pipeline/pgx_db.py +++ b/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]] = {} diff --git a/src/simurg/pipeline/pgx_pheno.py b/src/simurg/pipeline/pgx_pheno.py index 6e17a22..8e5224b 100644 --- a/src/simurg/pipeline/pgx_pheno.py +++ b/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", []) diff --git a/src/simurg/pipeline/pgx_star.py b/src/simurg/pipeline/pgx_star.py index 4a85460..f28e530 100644 --- a/src/simurg/pipeline/pgx_star.py +++ b/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 diff --git a/tests/unit/test_pgx_star.py b/tests/unit/test_pgx_star.py index 1cf6ec7..d41ec8a 100644 --- a/tests/unit/test_pgx_star.py +++ b/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"]