kemal
f885593e32
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.
133 lines
6.3 KiB
Python
133 lines
6.3 KiB
Python
"""``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))
|