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

@@ -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 |

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 `<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`).

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).")],

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"),

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]

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))

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 == []