From 9d3242828e42c5651a9184e5267bd1edfead32d8 Mon Sep 17 00:00:00 2001 From: Kemal Yaylali Date: Tue, 29 Sep 2026 16:03:27 +0100 Subject: [PATCH] PGx phenotypes: pgx phenotype|annotate (PGX-024..028, CLI-034) CPIC lookups by allele function counts or activity values (duplications N x activity, *5/*0 = 0, total-score fallback), Indeterminate when feasible solutions disagree or an allele is not in the guideline; per-sample simurg.pgx/1 JSON and diplotype_annotations.tsv with PGX-026 metadata lines. --- docs/messages.md | 1 + docs/validation/V2_pgx_star.md | 10 ++- src/simurg/cli/main.py | 41 ++++++++++ src/simurg/core/errors.py | 1 + src/simurg/pgx/phenotype.py | 118 +++++++++++++++++++++++++++ src/simurg/pipeline/pgx_pheno.py | 132 +++++++++++++++++++++++++++++++ tests/unit/test_pgx_phenotype.py | 36 +++++++++ 7 files changed, 338 insertions(+), 1 deletion(-) create mode 100644 src/simurg/pgx/phenotype.py create mode 100644 src/simurg/pipeline/pgx_pheno.py create mode 100644 tests/unit/test_pgx_phenotype.py diff --git a/docs/messages.md b/docs/messages.md index 2f9f672..3dce3ca 100644 --- a/docs/messages.md +++ b/docs/messages.md @@ -45,6 +45,7 @@ Generated from `src/simurg/core/errors.py` — do not edit by hand. | `E_PGX_BATCH_TOO_SMALL` | error | Too few (passing) samples in the batch for PGx CN calling | PCN-004 | | `E_PGX_CN_FAILED` | error | PGx CN calling failed for the sample | PCN-002 | | `E_PGX_CN_MISSING` | error | A sample has no PGx CN VCF (give --no-cn to run without) | PGX-011 | +| `E_PGX_GUIDELINE` | error | Guideline tables not present in the PGx DB | PGX-025 | | `E_PGX_SOURCE` | error | PGx source tables missing; run scripts/pgx/fetch_sources.sh | PGX-001 | | `E_PGX_TRAIN_TOO_SMALL` | error | Too few (passing) samples to train a PGx CN model | PCN-005 | | `E_PGX_VCF_UNINDEXED` | error | SNV VCFs for PGx must be bgzipped and indexed | PGX-011 | diff --git a/docs/validation/V2_pgx_star.md b/docs/validation/V2_pgx_star.md index 2522d96..8d30cd5 100644 --- a/docs/validation/V2_pgx_star.md +++ b/docs/validation/V2_pgx_star.md @@ -6,6 +6,10 @@ Status: **implemented for CPIC star-allele genes, interim validation** (2026-09- * `simurg pgx db build` → deterministic SQLite: 22 genes, 1354 alleles, 1897 normalised variant rows; GSA v4 PGx probe map 5002 links; two builds byte-identical. * `simurg pgx star-allele call` → `diplotype_summary.tsv`, `pgx_calls.json` (11 star-allele genes). +* `simurg pgx phenotype` (alias `annotate`) → per-sample `.pgx.json` (`simurg.pgx/1`) and + `diplotype_annotations.tsv` with `#` metadata lines (PGX-024 … PGX-028); CPIC activity scores (xN = N × activity, + *5/*0 = 0) and function-count lookups; disagreeing feasible solutions → `Indeterminate` plus per-solution + phenotypes. `tests/unit/test_pgx_phenotype.py`. * Unit tests: `tests/unit/test_pgx_db.py`, `tests/unit/test_pgx_star.py`. ## Data and truth @@ -54,5 +58,9 @@ so the per-call error is 1e-3 (GS ≥ 0.15) and errors are handled by the probe ## Caveats * 16 individuals, one array type, European (CEPH) only; masks from 14 individuals (spec: ≥ 100 samples). +* CPIC now defines the GRCh38 reference as its own allele where it differs from *1 (CYP2C19 *38): a sample + heterozygous at the *1-defining rs3758581 fits *1/*38 best (raw score 0.98), but *38 has no CPIC frequency, so the + prior keeps *1/*1 on top — which is what legacy GeT-RM reports. * Not yet: named-haplotype (G6PD), variant-based and tagging genes; CYP2D6 hybrids (*13, *36, *68); PharmVar-only - genes (CYP2A6, CYP3A4, CYP1A2, …); phenotype annotation (PGX-025) and the JSON schema (PGX-028). + genes (CYP2A6, CYP3A4, CYP1A2, …); DPWG tables; the full per-allele supporting/unresolved-variant detail of + PGX-027/028 (`supporting_variants.tsv`, `unresolved_variants.tsv`, `variants_not_covered.tsv`). diff --git a/src/simurg/cli/main.py b/src/simurg/cli/main.py index 394192f..41da1a8 100644 --- a/src/simurg/cli/main.py +++ b/src/simurg/cli/main.py @@ -869,6 +869,47 @@ def pgx_star_call_cmd( run_guarded(rep, lambda: run_star_call(o, rep)) +def _pgx_phenotype(calls: Path, db: Path, out: Path | None, guideline: str, force: bool, log_json: bool, + debug: bool, quiet: bool, name: str) -> None: # fmt: skip + from simurg.pipeline.pgx_pheno import PhenoOptions, run_phenotype + + if guideline not in ("cpic", "dpwg"): + raise typer.BadParameter("--guideline cpic|dpwg") + rep = _reporter(name, out or calls, log_json, debug, quiet) + o = PhenoOptions(calls=calls, db=db, out=out, guideline=guideline, force=force) + run_guarded(rep, lambda: run_phenotype(o, rep)) + + +@pgx_app.command("phenotype") +def pgx_phenotype_cmd( + calls: Annotated[Path, typer.Option("--calls", help="Output folder of pgx star-allele call (required).")], + db: Annotated[Path, typer.Option("--db", help="PGx DB (required).")], + out: Annotated[Path | None, typer.Option("--out", help="Output folder (default: --calls).")] = None, + guideline: Annotated[str, typer.Option("--guideline", help="cpic|dpwg")] = "cpic", + force: ForceOpt = False, + log_json: LogJson = False, + debug: Debug = False, + quiet: Quiet = False, +) -> None: + """Phenotypes from star-allele calls: per-sample PGx JSON and diplotype_annotations.tsv (CLI-034).""" + _pgx_phenotype(calls, db, out, guideline, force, log_json, debug, quiet, "pgx phenotype") + + +@pgx_app.command("annotate", hidden=True) +def pgx_annotate_cmd( + calls: Annotated[Path, typer.Option("--calls")], + db: Annotated[Path, typer.Option("--db")], + out: Annotated[Path | None, typer.Option("--out")] = None, + guideline: Annotated[str, typer.Option("--guideline")] = "cpic", + force: ForceOpt = False, + log_json: LogJson = False, + debug: Debug = False, + quiet: Quiet = False, +) -> None: + """Alias of pgx phenotype (CLI-034).""" + _pgx_phenotype(calls, db, out, guideline, force, log_json, debug, quiet, "pgx annotate") + + @pgx_cn_app.command("train") def pgx_cn_train_cmd( gtcs: Annotated[Path, typer.Option("--gtcs", help="Folder searched recursively for training GTCs (required).")], diff --git a/src/simurg/core/errors.py b/src/simurg/core/errors.py index 6758677..9b62008 100644 --- a/src/simurg/core/errors.py +++ b/src/simurg/core/errors.py @@ -174,6 +174,7 @@ CATALOGUE: dict[str, MessageDef] = dict( _m("E_PGX_CN_MISSING", "error", "A sample has no PGx CN VCF (give --no-cn to run without)", "PGX-011"), _m("E_PGX_VCF_UNINDEXED", "error", "SNV VCFs for PGx must be bgzipped and indexed", "PGX-011"), _m("E_PGX_VCF_UNSUITABLE", "error", "SNV VCF unsuitable for PGx (squashing/locus filters/foreign)", "PGX-011"), + _m("E_PGX_GUIDELINE", "error", "Guideline tables not present in the PGx DB", "PGX-025"), _m("E_PGX_SOURCE", "error", "PGx source tables missing; run scripts/pgx/fetch_sources.sh", "PGX-001"), # --- methylation control QC (MET-001 .. MET-018) ------------------------------------------------- _m("E_METH_NO_ASSAY", "error", "Probe section not found", "MET-002"), diff --git a/src/simurg/pgx/phenotype.py b/src/simurg/pgx/phenotype.py new file mode 100644 index 0000000..e16dd69 --- /dev/null +++ b/src/simurg/pgx/phenotype.py @@ -0,0 +1,118 @@ +"""Diplotype → phenotype lookup — spec PGX-024, PGX-025. Pure functions, no I/O. + +[DB: CPIC gene_result_lookup — the diplotype's allele functions (or, for activity-score genes, allele activity +values) as a count key, e.g. {"Normal function": 1, "No function": 1} or {"1.0": 1, "0.25": 1}.] A duplicated allele +``*Nxk`` contributes k × its activity value (CPIC activity-score rule for gene duplications); a whole-gene deletion +(``*5``, ``*0``) has activity 0 and "No function". When a key is absent from the table, activity-score genes fall back +to the total score (the phenotype of the largest tabulated total ≤ the score); anything else is ``Indeterminate`` with +the flag ``phenotype_not_in_guideline`` (PGX-024). +""" + +from __future__ import annotations + +import json +import re +from dataclasses import dataclass, field + + +@dataclass +class GuidelineTable: + gene: str + by_key: dict[str, tuple[str, str | None]] # canonical key JSON → (phenotype, total activity) + by_total: list[tuple[float, str]] # (total activity, phenotype), sorted + function: dict[str, str | None] # allele → clinical function + activity: dict[str, float | None] # allele → activity value + + @property + def activity_based(self) -> bool: + return bool(self.by_total) + + +@dataclass +class PhenotypeResult: + phenotype: str | None + activity_score: float | None + flags: list[str] = field(default_factory=list) + + +def _num(v: object) -> float | None: + try: + return float(str(v)) + except (TypeError, ValueError): + return None + + +def build_table(gene: str, rows: list[tuple[str, str | None, str]], alleles: list[tuple[str, str | None, str | None]] + ) -> GuidelineTable: # fmt: skip + """``rows``: (lookup_key JSON, total activity, phenotype); ``alleles``: (allele, function, activity).""" + 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) + 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} + act = {a: _num(v) for a, _, v in alleles} + return GuidelineTable(gene, by_key, sorted(totals.items()), func, act) + + +def _split(allele: str) -> tuple[str, int]: + m = re.fullmatch(r"(.+?)x(\d+|N)", allele) + if not m: + return allele, 1 + return m.group(1), (3 if m.group(2) == "N" else int(m.group(2))) + + +def phenotype(t: GuidelineTable, haplotypes: tuple[str, ...]) -> PhenotypeResult: + """Phenotype of one diplotype (haplotypes as reported, e.g. ("*2x2", "*4") or ("*1", "*5")).""" + 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) + 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}" + counts[s] = counts.get(s, 0) + 1 + key = json.dumps(counts, sort_keys=True) + total = sum(vals) + hit = t.by_key.get(key) + if hit is not None: + return PhenotypeResult(hit[0], total) + below = [p for tot, p in t.by_total if tot <= total + 1e-9] + if below: + return PhenotypeResult(below[-1], total) + return PhenotypeResult("Indeterminate", total, ["phenotype_not_in_guideline"]) + funcs: dict[str, int] = {} + for h in haplotypes: + core, n = _split(h) + 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)) + if hit is None: + return PhenotypeResult("Indeterminate", None, ["phenotype_not_in_guideline"]) + return PhenotypeResult(hit[0], None) + + +def gene_phenotype(t: GuidelineTable, feasible: list[tuple[str, ...]]) -> tuple[PhenotypeResult, list[str]]: + """PGX-025: a phenotype only when every feasible solution agrees; otherwise ``Indeterminate`` plus the + per-solution phenotypes.""" + res = [phenotype(t, s) for s in feasible] + if not res: + return PhenotypeResult(None, None, ["no_call"]), [] + phenos = sorted({r.phenotype or "" for r in res}) + if len(phenos) == 1: + acts = {r.activity_score for r in res} + return PhenotypeResult(res[0].phenotype, res[0].activity_score if len(acts) == 1 else None, + sorted({f for r in res for f in r.flags})), [] # fmt: skip + return PhenotypeResult("Indeterminate", None, ["feasible_solutions_disagree"]), [r.phenotype or "" for r in res] diff --git a/src/simurg/pipeline/pgx_pheno.py b/src/simurg/pipeline/pgx_pheno.py new file mode 100644 index 0000000..8fd95f2 --- /dev/null +++ b/src/simurg/pipeline/pgx_pheno.py @@ -0,0 +1,132 @@ +"""``simurg pgx phenotype`` (alias ``annotate``) — phenotypes and the per-sample PGx JSON (spec PGX-024 … PGX-028, +CLI-034). Reads the ``pgx star-allele call`` output folder and the PGx DB; can be re-run with another guideline. +""" + +from __future__ import annotations + +import hashlib +import json +import sqlite3 +from dataclasses import dataclass +from pathlib import Path +from typing import Any + +from simurg import __version__ +from simurg.core.errors import ExitCode, SimurgError +from simurg.core.log import Reporter +from simurg.pgx.phenotype import build_table, gene_phenotype + +JSON_SCHEMA = "simurg.pgx/1" +TSV_SCHEMA = "simurg.pgx_tsv/1" + + +@dataclass +class PhenoOptions: + calls: Path # folder written by pgx star-allele call + db: Path + out: Path | None = None + guideline: str = "cpic" + force: bool = False + + +def db_meta(db: Path) -> dict[str, str]: + con = sqlite3.connect(db) + meta = dict(con.execute("SELECT key, value FROM meta").fetchall()) + con.close() + return meta + + +def tsv_header(db: Path, meta: dict[str, str], guideline: str) -> list[str]: + """PGX-026: ``#`` metadata lines at the top of every PGx TSV.""" + sha = meta.get("content_sha256") or hashlib.sha256(db.read_bytes()).hexdigest() + return [f"#simurg_schema={TSV_SCHEMA}", f"#software_version={__version__}", "#genome_build=GRCh38", + f"#pgx_db={db.name};sha256={sha}", f"#allele_sources={meta.get('allele_sources', '')}", + f"#phenotype_sources={meta.get('phenotype_sources', '')}", f"#guideline={guideline}"] # fmt: skip + + +def write_meta_tsv(path: Path, header: list[str], cols: list[str], rows: list[dict[str, Any]]) -> None: + from simurg.core.atomic import atomic_write_text + + def cell(v: Any) -> str: + if v is None: + return "" + if isinstance(v, float): + if v != v: + return "" + txt = f"{v:.4f}".rstrip("0") + return txt + "0" if txt.endswith(".") else txt + return str(v).replace("\t", " ") + + lines = header + ["\t".join(cols)] + ["\t".join(cell(r.get(c)) for c in cols) for r in rows] + atomic_write_text(path, "\n".join(lines) + "\n") + + +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 + 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") + out = o.out or o.calls + out.mkdir(parents=True, exist_ok=True) + calls = json.loads(src.read_text()) + 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,)): + rows = con.execute("SELECT lookup_key, activity_score, phenotype FROM phenotypes WHERE gene = ? AND " + "guideline = ?", (gene, o.guideline)).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,))] + for (g,) in con.execute("SELECT gene FROM genes")} # fmt: skip + con.close() + sha = meta.get("content_sha256", "") + ann_rows = [] + for s in calls["samples"]: + genes = [] + for g in s["genes"]: + t = tables.get(g["gene"]) + feas = [tuple(x.split("/")) for x in g["feasible_solutions"]] + per: list[str] = [] + if t is None: + res = None + else: + res, per = gene_phenotype(t, feas) + flags = list(g.get("flags", [])) + (res.flags if res else []) + entry = {"gene": g["gene"], "call_type": "star_allele", "genotype": g["genotype"], + "feasible_solutions": g["feasible_solutions"], + "activity_score": res.activity_score if res else None, + "phenotype": res.phenotype if res else None, "per_solution_phenotypes": per, + "score": g["candidates"][0]["score"] if g["candidates"] else None, + "raw_score": g["candidates"][0]["raw_score"] if g["candidates"] else None, + "flags": sorted(set(flags)), "copy_number": g.get("copy_number"), + "candidate_solutions": [{"rank": i + 1, **c} for i, c in enumerate(g["candidates"])], + "alleles_tested": sorted(tested.get(g["gene"], []), key=_natural)} # fmt: skip + genes.append(entry) + ann_rows.append({"sample_id": s["sample_id"], "sample_name": s.get("sample_name", s["sample_id"]), + "gene": g["gene"], "guideline": o.guideline, "genotype": g["genotype"], + "activity_score": entry["activity_score"], "phenotype": entry["phenotype"], + "flags": ";".join(entry["flags"])}) # fmt: skip + doc = {"schema": JSON_SCHEMA, "software_version": __version__, "genome_build": "GRCh38", + "pgx_db": {"file": o.db.name, "sha256": sha}, + "allele_sources": json.loads(meta.get("allele_sources", "[]")), + "phenotype_sources": json.loads(meta.get("phenotype_sources", "[]")), "guideline": o.guideline, + "sample_id": s["sample_id"], "sample_name": s.get("sample_name", s["sample_id"]), + "genes": genes} # fmt: skip + atomic_write_text(out / f"{s['sample_id']}.pgx.json", json.dumps(doc, indent=1, ensure_ascii=False)) + cols = ["sample_id", "sample_name", "gene", "guideline", "genotype", "activity_score", "phenotype", "flags"] + ann_rows.sort(key=lambda r: (r["sample_id"], r["gene"])) + write_meta_tsv(out / "diplotype_annotations.tsv", tsv_header(o.db, meta, o.guideline), cols, ann_rows) + rep.info(f"phenotypes ({o.guideline}): {len(calls['samples'])} sample(s) -> {out}") + return ExitCode.OK + + +def _natural(a: str) -> tuple[Any, ...]: + import re + + return tuple(int(x) if x.isdigit() else x for x in re.split(r"(\d+)", a)) diff --git a/tests/unit/test_pgx_phenotype.py b/tests/unit/test_pgx_phenotype.py new file mode 100644 index 0000000..de95d65 --- /dev/null +++ b/tests/unit/test_pgx_phenotype.py @@ -0,0 +1,36 @@ +"""Diplotype → phenotype (spec PGX-024, PGX-025).""" + +from __future__ import annotations + +import json + +from simurg.pgx.phenotype import build_table, gene_phenotype, phenotype + +ACT_ROWS = [(json.dumps({"1.0": 2}), "2.0", "Normal Metabolizer"), (json.dumps({"0.0": 1, "1.0": 1}), "1.0", + "Intermediate Metabolizer"), (json.dumps({"0.0": 2}), "0.0", "Poor Metabolizer"), + (json.dumps({"1.0": 1, "2.0": 1}), "3.0", "Ultrarapid Metabolizer")] # fmt: skip +ACT_ALLELES = [("*1", "Normal function", "1.0"), ("*4", "No function", "0.0"), ("*10", "Decreased function", "0.25")] +FUNC_ROWS = [(json.dumps({"Normal function": 2}), "n/a", "Normal Metabolizer"), + (json.dumps({"Normal function": 1, "No function": 1}), "n/a", "Intermediate Metabolizer")] # fmt: skip +FUNC_ALLELES = [("*1", "Normal function", None), ("*2", "No function", None)] + + +def test_activity_genes() -> None: + t = build_table("CYP2D6", ACT_ROWS, ACT_ALLELES) + assert phenotype(t, ("*1", "*4")).phenotype == "Intermediate Metabolizer" + r = phenotype(t, ("*1x2", "*1")) # duplication: 2 × activity + assert r.phenotype == "Ultrarapid Metabolizer" and r.activity_score == 3.0 + assert phenotype(t, ("*1", "*5")).phenotype == "Intermediate Metabolizer" # *5 = whole-gene deletion + r = phenotype(t, ("*1", "*10")) # key absent → largest tabulated total ≤ 1.25 + assert r.phenotype == "Intermediate Metabolizer" and r.activity_score == 1.25 + r = phenotype(t, ("*1", "*99")) # PGX-024: unknown allele never crashes + assert r.phenotype == "Indeterminate" and "phenotype_not_in_guideline" in r.flags + + +def test_function_genes_and_disagreement() -> None: # PGX-025 + t = build_table("CYP2C19", FUNC_ROWS, FUNC_ALLELES) + assert phenotype(t, ("*2", "*1")).phenotype == "Intermediate Metabolizer" + res, per = gene_phenotype(t, [("*1", "*1"), ("*1", "*2")]) + assert res.phenotype == "Indeterminate" and per == ["Normal Metabolizer", "Intermediate Metabolizer"] + res, per = gene_phenotype(t, [("*1", "*2"), ("*2", "*1")]) + assert res.phenotype == "Intermediate Metabolizer" and per == []