PGx star alleles: pgx star-allele call (PGX-011..021), probe masks vs 1000G, GeT-RM validation (K-029)
- simurg/pgx/star.py: solver-free diplotype enumeration incl. CN forms (deletion, xN), wildcards, array-equivalence classes, raw_score/score with CPIC frequency priors; calibrated error model (spec PGX-012 amended). - pipeline/pgx_star.py: DB gene models, SNV + PGx CN VCF evidence, masks, diplotype_summary.tsv + pgx_calls.json. - scripts/pgx: derive_masks.py (array vs 1000G 30x over HTTPS, GS calibration), getrm_truth.py, validate_star_getrm.py. - GSA v4 PGx demo vs GeT-RM (16 CEPH): 100/109 gene calls; CYP2D6/2C19/2C9 41/45 (docs/validation/V2_pgx_star.md).
This commit is contained in:
kemal committed
1 parent
7e9ff28501
commit
194c087b24
15 files changed
+1168
-5
No files matched your search
@@ -838,6 +838,37 @@ def pgx_db_build_cmd(
|
||||
run_guarded(rep, lambda: run_db_build(o, rep))
|
||||
|
||||
|
||||
pgx_star_app = typer.Typer(help="Star-allele calling (PGX-011).")
|
||||
pgx_app.add_typer(pgx_star_app, name="star-allele")
|
||||
|
||||
|
||||
@pgx_star_app.command("call")
|
||||
def pgx_star_call_cmd(
|
||||
vcfs: Annotated[list[Path], typer.Option("--vcfs", help="Folder(s) with SNV + PGx CN VCFs (repeatable).")],
|
||||
db: Annotated[Path, typer.Option("--db", help="PGx DB from pgx db build (required).")],
|
||||
out: Annotated[Path, typer.Option("--out", help="Output folder (required).")],
|
||||
population: Annotated[str, typer.Option("--population", help="Frequency prior population.")] = "global",
|
||||
max_candidates: Annotated[int, typer.Option("--max-candidates", min=1)] = 10,
|
||||
feasible_ratio: Annotated[float, typer.Option("--feasible-ratio", min=0.0, max=1.0)] = 0.5,
|
||||
low_gs: Annotated[float, typer.Option("--low-gs")] = 0.15,
|
||||
no_cn: Annotated[bool, typer.Option("--no-cn", help="Run without PGx CN VCFs.")] = False,
|
||||
allow_foreign_vcf: Annotated[bool, typer.Option("--allow-foreign-vcf")] = False,
|
||||
gene: Annotated[list[str] | None, typer.Option("--gene", help="Restrict to gene(s).")] = None,
|
||||
force: ForceOpt = False,
|
||||
log_json: LogJson = False,
|
||||
debug: Debug = False,
|
||||
quiet: Quiet = False,
|
||||
) -> None:
|
||||
"""Call star-allele diplotypes (CLI-033, PGX-012 ... PGX-021)."""
|
||||
from simurg.pipeline.pgx_star import StarOptions, run_star_call
|
||||
|
||||
rep = _reporter("pgx star-allele call", out, log_json, debug, quiet)
|
||||
o = StarOptions(vcfs=vcfs, db=db, out=out, population=population, max_candidates=max_candidates,
|
||||
feasible_ratio=feasible_ratio, low_gs=low_gs, no_cn=no_cn, allow_foreign_vcf=allow_foreign_vcf,
|
||||
genes=gene or [], force=force) # fmt: skip
|
||||
run_guarded(rep, lambda: run_star_call(o, rep))
|
||||
|
||||
|
||||
@pgx_cn_app.command("train")
|
||||
def pgx_cn_train_cmd(
|
||||
gtcs: Annotated[Path, typer.Option("--gtcs", help="Folder searched recursively for training GTCs (required).")],
|
||||
|
||||
@@ -171,6 +171,9 @@ CATALOGUE: dict[str, MessageDef] = dict(
|
||||
_m("W_PGX_NO_ALT_CONTIG", "warning", "FASTA lacks a contig carrying a PGx region; REF written as N", "IN-062"),
|
||||
_m("W_PGX_BATCH_COMPOSITION", "warning", "Batch spans several lots/batches or > 14 days of scans", "PCN-006"),
|
||||
_m("E_PGX_CN_FAILED", "error", "PGx CN calling failed for the sample", "PCN-002"),
|
||||
_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_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"),
|
||||
|
||||
@@ -0,0 +1,172 @@
|
||||
"""Star-allele / haplotype calling from array genotypes — spec PGX-012 … PGX-020. Pure functions, no I/O.
|
||||
|
||||
Evidence model (PGX-012) [OURS]: each defining variant v of a gene is ``called`` (alt dosage g, GenCall score GS,
|
||||
BAF), ``no_call``, ``masked`` or ``not_covered``; only called variants contribute. With diploid CN 2 the genotype
|
||||
likelihood is P(g | d) = 1 − ε for the expected alt dosage d and ε/2 otherwise, with ε = 1e-3 for GS ≥ 0.15 and
|
||||
0.05 below (calibrated: after probe masking, 0 of 9740 GSA v4 PGx calls at DB variants disagreed with 1000G 30x in any
|
||||
GS bin — errors are probe-specific and handled by masks, not by GS; K-029). When the gene's copy number n ≠ 2
|
||||
(targeted CN), the BAF likelihood N(b; k/n, σ_b) replaces the
|
||||
genotype term (k = expected alt copies). Haplotype positions with an IUPAC value (CPIC core definitions) are wildcards:
|
||||
the best-fitting dosage is used.
|
||||
|
||||
Candidates (PGX-013) [DIFF: solver-free]: every unordered pair of defined haplotypes, and for n ≠ 2 the structural
|
||||
forms of CYP2D6-like genes: n = 1 → one whole-gene deletion haplotype; n = 0 → two; n ≥ 3 → one haplotype duplicated
|
||||
(``*1x2``, n − 1 copies). Diplotypes that no called variant separates are array-equivalent (PGX-017) and form one
|
||||
solution, named after the member with the highest prior. Scores (PGX-014): raw_score = softmax(LL); score =
|
||||
softmax(LL + log π), π = Σ over members of f(h1)·f(h2) (×2 for h1 ≠ h2), f floored at 1e-5.
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import math
|
||||
from dataclasses import dataclass, field
|
||||
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
|
||||
F64 = npt.NDArray[np.float64]
|
||||
FREQ_FLOOR = 1e-5
|
||||
CALLED, NO_CALL, MASKED, NOT_COVERED = 0, 1, 2, 3
|
||||
STATE_NAMES = {CALLED: "called", NO_CALL: "no_call", MASKED: "masked", NOT_COVERED: "not_on_array"}
|
||||
|
||||
|
||||
@dataclass
|
||||
class GeneModel:
|
||||
gene: str
|
||||
variants: list[tuple[str, int, str, str]] # (chrom, pos, ref, alt), one per defining variant
|
||||
variant_ids: list[tuple[str, str | None]] # (hgvs, rsid)
|
||||
haplotypes: list[str]
|
||||
carries: npt.NDArray[np.int8] # (H, V): 1 carries alt, 0 reference, −1 wildcard
|
||||
freq: F64 # (H,) population frequency of each haplotype (prior)
|
||||
deletion_allele: str | None = None # named whole-gene deletion (CYP2D6 *5), or "*0"
|
||||
reference: str = ""
|
||||
|
||||
|
||||
@dataclass
|
||||
class Observations:
|
||||
state: npt.NDArray[np.int8] # (V,)
|
||||
dosage: npt.NDArray[np.int8] # alt dosage (0 … ploidy), −1 when not called
|
||||
gs: F64
|
||||
baf: F64
|
||||
|
||||
|
||||
@dataclass
|
||||
class Solution:
|
||||
haplotypes: tuple[str, ...] # representative, e.g. ("*1", "*4") or ("*2x2", "*4")
|
||||
ll: float
|
||||
prior: float
|
||||
members: list[tuple[str, ...]] = field(default_factory=list)
|
||||
raw_score: float = 0.0
|
||||
score: float = 0.0
|
||||
copy_number: int = 2
|
||||
|
||||
@property
|
||||
def genotype(self) -> str:
|
||||
return "/".join(self.haplotypes)
|
||||
|
||||
|
||||
ERR_CALLED = 1e-3
|
||||
ERR_LOW_GS = 0.05
|
||||
LOW_GS = 0.15
|
||||
|
||||
|
||||
def gs_error(gs: F64) -> F64:
|
||||
"""PGX-012: per-call error rate from the GenCall score (calibrated, see the module docstring)."""
|
||||
return np.where(np.nan_to_num(gs, nan=0.0) >= LOW_GS, ERR_CALLED, ERR_LOW_GS)
|
||||
|
||||
|
||||
def _dosage_ll(obs: Observations, expected: npt.NDArray[np.int16], wild: npt.NDArray[np.bool_]) -> F64:
|
||||
"""(D, V) log-likelihood of each candidate's expected dosage (CN 2 genotype model)."""
|
||||
eps = gs_error(obs.gs)[None, :]
|
||||
g = obs.dosage[None, :].astype(np.int16)
|
||||
hit = (expected == g) | wild
|
||||
return np.where(hit, np.log(1 - eps), np.log(eps / 2))
|
||||
|
||||
|
||||
def _baf_ll(obs: Observations, frac_lo: F64, frac_hi: F64, sigma: float) -> F64:
|
||||
"""(D, V) BAF log-likelihood for expected alt fractions (a wildcard spans [lo, hi]: the closest value is used)."""
|
||||
b = np.nan_to_num(obs.baf, nan=0.5)[None, :]
|
||||
near = np.clip(b, frac_lo, frac_hi)
|
||||
return np.log(0.98 * np.exp(-0.5 * ((b - near) / sigma) ** 2) / (sigma * math.sqrt(2 * math.pi)) + 0.02)
|
||||
|
||||
|
||||
def candidates(m: GeneModel, cn: int | None) -> list[tuple[tuple[int, ...], tuple[str, ...]]]:
|
||||
"""(haplotype index per gene copy, display names) for every candidate diplotype (PGX-013)."""
|
||||
h = len(m.haplotypes)
|
||||
out: list[tuple[tuple[int, ...], tuple[str, ...]]] = []
|
||||
if cn is None or cn == 2:
|
||||
for i in range(h):
|
||||
for j in range(i, h):
|
||||
out.append(((i, j), (m.haplotypes[i], m.haplotypes[j])))
|
||||
return out
|
||||
dele = m.deletion_allele or "*0"
|
||||
if cn == 0:
|
||||
return [((), (dele, dele))]
|
||||
if cn == 1:
|
||||
return [((i,), (m.haplotypes[i], dele)) for i in range(h)]
|
||||
extra = cn - 2 # copies beyond two: one haplotype duplicated (e.g. CN 3 → *2x2/*4)
|
||||
for i in range(h):
|
||||
for j in range(h):
|
||||
out.append(((i,) * (extra + 1) + (j,), (f"{m.haplotypes[i]}x{extra + 1}", m.haplotypes[j])))
|
||||
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,
|
||||
) -> 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)
|
||||
if not cands:
|
||||
return [], []
|
||||
carries = m.carries.astype(np.int16)
|
||||
lo = np.zeros((len(cands), carries.shape[1]))
|
||||
hi = np.zeros_like(lo)
|
||||
for k, (ix, _) in enumerate(cands):
|
||||
for i in ix:
|
||||
c = carries[i]
|
||||
lo[k] += np.where(c < 0, 0, c)
|
||||
hi[k] += np.where(c < 0, 1, c)
|
||||
ncp = np.array([max(1, len(ix)) for ix, _ in cands], dtype=np.float64)[:, None]
|
||||
if ploidy == 1 or cn is None or cn == 2:
|
||||
wild = hi != lo
|
||||
ll_v = _dosage_ll(obs, lo.astype(np.int16), wild)
|
||||
if wild.any(): # a wildcard may also take the higher dosage
|
||||
ll_v = np.maximum(ll_v, _dosage_ll(obs, hi.astype(np.int16), wild))
|
||||
else:
|
||||
ll_v = _baf_ll(obs, lo / ncp, hi / ncp, baf_sigma)
|
||||
ll = (ll_v * called[None, :]).sum(axis=1)
|
||||
f = np.maximum(m.freq, FREQ_FLOOR)
|
||||
# array-equivalence classes: identical expected dosages at every called variant (PGX-017)
|
||||
keys: dict[bytes, list[int]] = {}
|
||||
for k in range(len(cands)):
|
||||
sig = np.round(np.concatenate([lo[k, called], hi[k, called]]) / ncp[k, 0], 6).tobytes()
|
||||
keys.setdefault(sig + str(len(cands[k][0])).encode(), []).append(k)
|
||||
sols = []
|
||||
for members in keys.values():
|
||||
pri = []
|
||||
for k in members:
|
||||
ix = cands[k][0]
|
||||
if not ix:
|
||||
pri.append(1e-6)
|
||||
continue
|
||||
p = float(np.prod([f[i] for i in ix]))
|
||||
if len(set(ix)) > 1:
|
||||
p *= 2
|
||||
pri.append(p)
|
||||
best = members[int(np.argmax(pri))]
|
||||
sols.append(Solution(cands[best][1], float(ll[best]), float(sum(pri)),
|
||||
[cands[k][1] for k in members], copy_number=cn if cn is not None else 2)) # fmt: skip
|
||||
lls = np.array([s.ll for s in sols])
|
||||
lp = lls + np.log([s.prior for s in sols])
|
||||
raw = np.exp(lls - lls.max())
|
||||
sc = np.exp(lp - lp.max())
|
||||
raw, sc = raw / raw.sum(), sc / sc.sum()
|
||||
for s, r, v in zip(sols, raw, sc, strict=True):
|
||||
s.raw_score, s.score = float(r), float(v)
|
||||
sols.sort(key=lambda s: (-s.score, -s.raw_score, s.genotype))
|
||||
top = sols[0].score
|
||||
feasible = [s for s in sols if s.score >= feasible_ratio * top]
|
||||
return sols[: max(max_candidates, len(feasible))], feasible
|
||||
@@ -0,0 +1,256 @@
|
||||
"""``simurg pgx star-allele call`` — star alleles from SNV VCFs, PGx DB and PGx CN VCFs (spec PGX-011 … PGX-021,
|
||||
CLI-033). I/O and orchestration; the scoring lives in ``simurg.pgx.star``.
|
||||
|
||||
Scope of this version: genes whose result basis is ``star_allele`` in the DB. Named-haplotype (G6PD), variant-based
|
||||
and tagging genes, sub-allele collapsing and hybrid CYP2D6 structures follow (K-029).
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import gzip
|
||||
import json
|
||||
import math
|
||||
import sqlite3
|
||||
from collections import defaultdict
|
||||
from dataclasses import dataclass, field
|
||||
from pathlib import Path
|
||||
from typing import Any
|
||||
|
||||
import numpy as np
|
||||
|
||||
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
|
||||
|
||||
CN_REGION = {"CYP2D6": "CYP2D6.p5"} # PGX-013: total gene copies from the 5′ region
|
||||
DELETION_ALLELE = {"CYP2D6": "*5"}
|
||||
|
||||
|
||||
@dataclass
|
||||
class StarOptions:
|
||||
vcfs: list[Path]
|
||||
db: Path
|
||||
out: Path
|
||||
population: str = "global"
|
||||
max_candidates: int = 10
|
||||
feasible_ratio: float = 0.5
|
||||
low_gs: float = 0.15
|
||||
no_cn: bool = False
|
||||
allow_foreign_vcf: bool = False
|
||||
genes: list[str] = field(default_factory=list)
|
||||
force: bool = False
|
||||
|
||||
|
||||
def load_models(db: Path, population: str, genes: list[str] | None = None) -> dict[str, GeneModel]:
|
||||
con = sqlite3.connect(db)
|
||||
q = "SELECT gene FROM genes WHERE result_basis = 'star_allele' ORDER BY gene"
|
||||
out = {}
|
||||
for (gene,) in con.execute(q).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 "
|
||||
"ORDER BY allele", (gene,)).fetchall() # fmt: skip
|
||||
names = [a for a, _ in alle]
|
||||
ref = next((a for a, k in alle if k == "reference"), names[0] if names else "")
|
||||
rows = con.execute("SELECT allele, chrom, pos, ref, alt, hgvs, rsid, location_id FROM allele_variants "
|
||||
"WHERE gene = ?", (gene,)).fetchall() # fmt: skip
|
||||
keys = sorted({(r[1], r[2], r[3], r[4]) for r in rows}, key=lambda k: (k[0], k[1], k[2], k[3]))
|
||||
kidx = {k: i for i, k in enumerate(keys)}
|
||||
ids: list[tuple[str, str | None]] = [("", None)] * len(keys)
|
||||
loc_keys: dict[int, list[int]] = defaultdict(list)
|
||||
carries = np.zeros((len(names), len(keys)), dtype=np.int8)
|
||||
hidx = {n: i for i, n in enumerate(names)}
|
||||
for allele, chrom, pos, r_, a_, hgvs, rsid, lid in rows:
|
||||
k = kidx[(chrom, pos, r_, a_)]
|
||||
ids[k] = (hgvs, rsid)
|
||||
loc_keys[lid].append(k)
|
||||
if allele in hidx:
|
||||
carries[hidx[allele], k] = 1
|
||||
for allele, lid in con.execute("SELECT allele, location_id FROM allele_wildcards WHERE gene = ?", (gene,)):
|
||||
if allele in hidx:
|
||||
for k in loc_keys.get(lid, []):
|
||||
carries[hidx[allele], k] = -1
|
||||
freq = np.full(len(names), np.nan)
|
||||
fr: dict[str, list[float]] = defaultdict(list)
|
||||
for allele, pop, f in con.execute("SELECT allele, population, freq FROM frequencies WHERE gene = ?", (gene,)):
|
||||
if population in ("global", pop):
|
||||
fr[allele].append(float(f))
|
||||
for i, n in enumerate(names):
|
||||
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)
|
||||
con.close()
|
||||
return out
|
||||
|
||||
|
||||
def read_vcf_sites(
|
||||
path: Path, wanted: dict[str, set[int]]
|
||||
) -> tuple[dict[tuple[str, int], list[dict[str, Any]]], dict[str, str]]:
|
||||
"""Records at the wanted positions, and the ``##`` header meta (``##key=value``)."""
|
||||
recs: dict[tuple[str, int], list[dict[str, Any]]] = defaultdict(list)
|
||||
meta: dict[str, str] = {}
|
||||
with gzip.open(path, "rt") as fh:
|
||||
for ln in fh:
|
||||
if ln.startswith("##"):
|
||||
k, _, v = ln[2:].rstrip("\n").partition("=")
|
||||
meta.setdefault(k, v)
|
||||
continue
|
||||
if ln.startswith("#"):
|
||||
continue
|
||||
c, p, _ = ln.split("\t", 2)
|
||||
ps = wanted.get(c)
|
||||
if ps is None or int(p) not in ps:
|
||||
continue
|
||||
f = ln.rstrip("\n").split("\t")
|
||||
fmt = f[8].split(":")
|
||||
val = dict(zip(fmt, f[9].split(":"), strict=False))
|
||||
recs[(c, int(p))].append({"id": f[2], "ref": f[3], "alts": f[4].split(","), "gt": val.get("GT", "./."),
|
||||
"gs": val.get("GS"), "baf": val.get("BAF")}) # fmt: skip
|
||||
return recs, meta
|
||||
|
||||
|
||||
def observations(
|
||||
m: GeneModel, recs: dict[tuple[str, int], list[dict[str, Any]]], masked: set[str] | None = None
|
||||
) -> tuple[Observations, int]:
|
||||
"""PGX-012 evidence per defining variant; returns the observations and the ploidy seen in the GTs. A record with
|
||||
any masked probe (PGX-010) makes the variant ``masked``."""
|
||||
v = len(m.variants)
|
||||
state = np.full(v, NOT_COVERED, dtype=np.int8)
|
||||
dosage = np.full(v, -1, dtype=np.int8)
|
||||
gs = np.full(v, np.nan)
|
||||
baf = np.full(v, np.nan)
|
||||
ploidy = 2
|
||||
for k, (c, p, r_, a_) in enumerate(m.variants):
|
||||
for rec in recs.get((c, p), []):
|
||||
if rec["ref"] != r_ or a_ not in rec["alts"]:
|
||||
continue
|
||||
ai = rec["alts"].index(a_) + 1
|
||||
if masked and any(pn in masked for pn in rec["id"].replace(";", ",").split(",")):
|
||||
state[k] = MASKED
|
||||
break
|
||||
alleles = rec["gt"].replace("|", "/").split("/")
|
||||
ploidy = len(alleles)
|
||||
if any(a in (".", "") for a in alleles):
|
||||
state[k] = NO_CALL
|
||||
break
|
||||
state[k] = CALLED
|
||||
dosage[k] = sum(int(a) == ai for a in alleles)
|
||||
gs[k] = float(rec["gs"]) if rec["gs"] not in (None, ".") else math.nan
|
||||
baf[k] = float(rec["baf"]) if rec["baf"] not in (None, ".") else math.nan
|
||||
if len(rec["alts"]) > 1 and not math.isnan(baf[k]):
|
||||
baf[k] = math.nan # BAF of a merged multi-allelic record is not per-ALT
|
||||
break
|
||||
return Observations(state, dosage, gs, baf), ploidy
|
||||
|
||||
|
||||
def read_cn_vcf(path: Path) -> dict[str, int | None]:
|
||||
out: dict[str, int | None] = {}
|
||||
with gzip.open(path, "rt") as fh:
|
||||
for ln in fh:
|
||||
if ln.startswith("#"):
|
||||
continue
|
||||
f = ln.rstrip("\n").split("\t")
|
||||
region = f[7].split("REGION=")[1].split(";")[0]
|
||||
cn = f[9].split(":")[0]
|
||||
out[region] = None if cn == "." or "NoCall" in f[6] else int(cn)
|
||||
return out
|
||||
|
||||
|
||||
def _discover(dirs: list[Path]) -> tuple[dict[str, Path], dict[str, Path]]:
|
||||
snv, cn = {}, {}
|
||||
for d in dirs:
|
||||
for p in sorted(d.rglob("*.vcf.gz")):
|
||||
if p.name.endswith(".pgx_cn.vcf.gz"):
|
||||
cn[p.name[: -len(".pgx_cn.vcf.gz")]] = p
|
||||
elif p.name.endswith(".snv.vcf.gz"):
|
||||
snv[p.name[: -len(".snv.vcf.gz")]] = p
|
||||
return snv, cn
|
||||
|
||||
|
||||
def run_star_call(o: StarOptions, rep: Reporter) -> int:
|
||||
"""``simurg pgx star-allele call`` (CLI-033)."""
|
||||
from simurg.core.atomic import check_outputs_absent
|
||||
from simurg.pipeline.qc import write_tsv
|
||||
|
||||
snv, cnv = _discover(o.vcfs)
|
||||
if not snv:
|
||||
raise SimurgError("E_NO_INPUT", f"no *.snv.vcf.gz under {', '.join(map(str, o.vcfs))}")
|
||||
missing_cn = sorted(set(snv) - set(cnv))
|
||||
if missing_cn and not o.no_cn:
|
||||
raise SimurgError("E_PGX_CN_MISSING", f"{len(missing_cn)} sample(s) have no .pgx_cn.vcf.gz (e.g. "
|
||||
f"{missing_cn[0]})", hint="run pgx copy-number call, or give --no-cn") # fmt: skip
|
||||
for vp in snv.values():
|
||||
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)
|
||||
mcon = sqlite3.connect(o.db)
|
||||
masks: dict[str, set[str]] = defaultdict(set)
|
||||
for man, probe in mcon.execute("SELECT manifest, probe_name FROM masks"):
|
||||
masks[man].add(probe)
|
||||
mcon.close()
|
||||
wanted: dict[str, set[int]] = defaultdict(set)
|
||||
for m in models.values():
|
||||
for c, p, _, _ in m.variants:
|
||||
wanted[c].add(p)
|
||||
rows = []
|
||||
calls_json = []
|
||||
for sid, path in sorted(snv.items()):
|
||||
recs, meta = read_vcf_sites(path, wanted)
|
||||
opts = meta.get("simurg_options", "")
|
||||
if not opts and not o.allow_foreign_vcf:
|
||||
raise SimurgError("E_PGX_VCF_UNSUITABLE", f"{path}: not a SimurgArray VCF (no ##simurg_options)",
|
||||
hint="give --allow-foreign-vcf to use it anyway") # fmt: skip
|
||||
kv = dict(x.partition("=")[::2] for x in opts.split(";") if x)
|
||||
if (
|
||||
kv.get("unsquash_duplicates") == "true"
|
||||
or kv.get("exclude_loci", "none").lower() != "none"
|
||||
or kv.get("include_loci", "none").lower() != "none"
|
||||
): # PGX-011
|
||||
raise SimurgError("E_PGX_VCF_UNSUITABLE", f"{path}: made with unsquashed duplicates or a locus filter")
|
||||
cn = {} if o.no_cn or sid not in cnv else read_cn_vcf(cnv[sid])
|
||||
masked = masks.get(meta.get("simurg_manifest", "").split(";")[0], set())
|
||||
sample_name = sid
|
||||
genes_out = []
|
||||
for gene, m in sorted(models.items()):
|
||||
obs, ploidy = observations(m, recs, masked)
|
||||
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
|
||||
flags = []
|
||||
low = [k for k in range(len(m.variants)) if obs.state[k] == CALLED and obs.gs[k] < o.low_gs]
|
||||
if low:
|
||||
flags.append("review_alternatives")
|
||||
no_calls = [k for k in range(len(m.variants)) if obs.state[k] == NO_CALL]
|
||||
genotype = ";".join(s.genotype for s in feasible)
|
||||
for rank, s in enumerate(sols, start=1):
|
||||
rows.append({"sample_id": sid, "sample_name": sample_name, "solution_id": rank, "gene": gene,
|
||||
"call_type": "star_allele", "solution": s.genotype,
|
||||
"array_equivalent_alleles_1": ";".join(sorted({mm[0] for mm in s.members})),
|
||||
"array_equivalent_alleles_2": ";".join(sorted({mm[-1] for mm in s.members})),
|
||||
"no_call_variants": ";".join(m.variant_ids[k][0] for k in no_calls),
|
||||
"alleles_impacted_by_no_call": "", "background_alleles": "", "variant_change": "",
|
||||
"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 "",
|
||||
"flags": ";".join(flags + (["feasible"] if s in feasible else []))}) # fmt: skip
|
||||
genes_out.append({"gene": gene, "genotype": genotype, "feasible_solutions": [s.genotype for s in feasible],
|
||||
"candidates": [_sol_json(s) for s in sols], "copy_number": gcn,
|
||||
"flags": flags}) # fmt: skip
|
||||
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",
|
||||
"variant_change", "score", "raw_score", "copy_number_solution", "flags"] # fmt: skip
|
||||
o.out.mkdir(parents=True, exist_ok=True)
|
||||
write_tsv(o.out / "diplotype_summary.tsv", cols, rows)
|
||||
(o.out / "pgx_calls.json").write_text(json.dumps({"software": f"SimurgArray {__version__}", "db": o.db.name,
|
||||
"samples": calls_json}, indent=1)) # fmt: skip
|
||||
rep.info(f"star alleles: {len(calls_json)} sample(s), {len(models)} gene(s)")
|
||||
return ExitCode.OK
|
||||
|
||||
|
||||
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
|
||||
Reference in new issue
Block a user