From 194c087b240562cc6276396de74a0113b5842679 Mon Sep 17 00:00:00 2001 From: Kemal Yaylali Date: Tue, 29 Sep 2026 14:08:03 +0100 Subject: [PATCH] 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). --- SIMURGARRAY_FUNCTIONAL_SPEC.md | 5 +- docs/known_unknowns.md | 1 + docs/messages.md | 3 + docs/validation/V2_pgx_star.md | 58 +++++ docs/validation/V2_pgx_star_getrm16.json | 240 +++++++++++++++++++++ scripts/pgx/build_cn_regions.py | 4 +- scripts/pgx/derive_masks.py | 157 ++++++++++++++ scripts/pgx/getrm_truth.py | 51 +++++ scripts/pgx/validate_cn_truth.py | 6 +- scripts/pgx/validate_star_getrm.py | 118 +++++++++++ src/simurg/cli/main.py | 31 +++ src/simurg/core/errors.py | 3 + src/simurg/pgx/star.py | 172 +++++++++++++++ src/simurg/pipeline/pgx_star.py | 256 +++++++++++++++++++++++ tests/unit/test_pgx_star.py | 68 ++++++ 15 files changed, 1168 insertions(+), 5 deletions(-) create mode 100644 docs/validation/V2_pgx_star.md create mode 100644 docs/validation/V2_pgx_star_getrm16.json create mode 100644 scripts/pgx/derive_masks.py create mode 100644 scripts/pgx/getrm_truth.py create mode 100644 scripts/pgx/validate_star_getrm.py create mode 100644 src/simurg/pgx/star.py create mode 100644 src/simurg/pipeline/pgx_star.py create mode 100644 tests/unit/test_pgx_star.py diff --git a/SIMURGARRAY_FUNCTIONAL_SPEC.md b/SIMURGARRAY_FUNCTIONAL_SPEC.md index 4a6e3ab..68ef106 100644 --- a/SIMURGARRAY_FUNCTIONAL_SPEC.md +++ b/SIMURGARRAY_FUNCTIONAL_SPEC.md @@ -2759,7 +2759,10 @@ about **algorithms or numerics**, the plan wins. Each such conflict found during In a region whose CN is n ≠ 2, the variant copy number `k ∈ {0..n}` is the maximum-likelihood value under the observed BAF ~ N(k/n, σ²), with σ estimated per batch from AB clusters. It is reported as `variant_copy_number`. - GS enters as a per-call error rate: ε(GS) = 0.5·(1 − GS)², clipped to [1e-4, 0.25] [OURS]. (src: B-30, E-41) + GS enters as a per-call error rate: ε = 1e-3 for GS ≥ 0.15 and 0.05 below [OURS, calibrated: after probe masking + (PGX-010), 0 of 9 740 GSA v4 PGx calls at DB variants disagreed with 1000G 30x in any GS bin, so errors are + probe-specific and belong in the mask; the earlier ε(GS) = 0.5·(1 − GS)² made correct low-GS calls such as + CYP2C9 rs1799853 at GS 0.16 lose to the prior; K-029]. (src: B-30, E-41) - **PGX-013** `§11 · v2` **Candidate enumeration [DIFF: solver-free].** - For each gene, SimurgArray enumerates diplotypes from the gene's haplotypes: star alleles, `*0`, duplications diff --git a/docs/known_unknowns.md b/docs/known_unknowns.md index a2ccb06..6c30ee3 100644 --- a/docs/known_unknowns.md +++ b/docs/known_unknowns.md @@ -4,6 +4,7 @@ Plan §20 lists the open questions; this file records what the implementation ha | ID | Date | Topic | Finding | Status | |---|---|---|---|---| +| K-029 | 2026-09-29 | PGx star alleles (v2) | (1) PGX-012 error model: ε(GS) = 0.5·(1 − GS)² was replaced by ε = 1e-3 (GS ≥ 0.15) / 0.05 — after probe masking 0 of 9 740 GSA v4 PGx calls disagreed with 1000G 30x at any GS, and the old model made correct low-GS calls lose to the prior (spec amended). (2) PGX-010 masks derived from 14 individuals (spec: ≥ 100); a variant absent from the filtered 1000G panel counts as hom-ref only if rare (CYP2B6 rs2279343 is filtered out, not absent). (3) CPIC core-allele definitions use IUPAC wildcards; alleles without definitions (NUDT15 *2, structural CYP2D6) are not callable from SNVs. (4) GeT-RM uses legacy UGT1A1/SLCO1B1 names. Result vs GeT-RM (16 CEPH): 100/109 gene calls, CYP2D6/2C19/2C9 41/45. | Open: hybrids, G6PD/variant-based genes, PharmVar genes, phenotypes; larger/independent truth (`docs/validation/V2_pgx_star.md`). | | K-028 | 2026-09-29 | PGx CN (v2) | Deviations from PCN-009: the per-probe reference is the training cohort's median stored in the model (a calling-batch median is not CN 2 where deletions are common, e.g. GSTM1 null ≈ 50 %); no GC-correction step; the CN mixture is ours: µ_c = δ + γ·log2((c+1)/3) per region with a free lower CN-0 component and Hardy–Weinberg-constrained weights (probe compression differs by region: γ ≈ 0.19 CYP2D6.p5 … 0.56 GSTT1). Interim validation on the GSA v4 PGx demo — whose 36 arrays turned out to be 30 CEPH/1000G individuals — against NYGC 30x SV RD_CN: 224/240 region calls (GSTM1, UGT2B17 30/30), but in-sample (trained on the same arrays, < 96 samples, priors chosen with this truth in view). | Open: validate on an independent GSA-PGx cohort; CYP2D6 via GeT-RM star alleles (`docs/validation/V2_pgx_cn.md`). | | K-027 | 2026-09-29 | CN/LOH (v2) | CYTO-014 noise gate: the summary `logr_dev` also measures cluster-file mismatch (GSE165845: 0.246, all samples `SampleQuality`) while the LRR used after the batch reference has SD ≈ 0.1; the gate now uses the SD of the normalised LRR (the JSON still reports the summary value). Spec: built-in mask/cytobands from UCSC; the annotation bundle uses UCSC ncbiRefSeq instead of RefSeq GFF + Ensembl GTF + HGNC. Real-data result: constitutional DEL/DUP calls 74–77 % in NYGC truth, mosaic calls mostly false positives. Pi tuning: mosaic prior / switch distance barely matter; CYTO-041 (noise-normalised BAF split, complete-loss promotion) cuts PASS mosaic calls by ~30 % without losing implanted ≥ 15 % mosaics. Two limits remain: (1) without phasing, < 20 % mosaic losses are not separable from regional BAF noise; (2) common CNVs (carrier frequency > 20 %) are almost never detected (1/236 DEL) because the median-LRR batch reference absorbs them; rare DEL (< 5 %) 38 %. | Open: rare-event sensitivity as primary metric; phasing-based mosaic BAF (`docs/validation/V2_cn.md`). | | K-026 | 2026-09-29 | Methylation QC (v2) | Spec deviations: vendor EPIC manifests put a 5th field on some control lines (≥ 4 fields accepted); detection p follows sesame's pOOBAH (1 − max(F(M), F(U))) rather than the ambiguous "1 − ECDF(M + U)"; minfi-style `Sentrix_ID`/`Sentrix_Position` sheet keys accepted. Against sesame 1.28.1 on minfiDataEPIC: β r = 0.99998 (max \|Δβ\| 0.011), detection calls 99.1 % identical, noob signals within 0.5 % (median) / 5 % (max). The systematic remainder is unexplained; it is not the negative controls (sesame did not attach the supplied control table). | Open; `docs/validation/V2_methylation.md`. Next step: dump sesame's background vectors and noob parameters to find the difference. | diff --git a/docs/messages.md b/docs/messages.md index f010884..2f9f672 100644 --- a/docs/messages.md +++ b/docs/messages.md @@ -44,8 +44,11 @@ Generated from `src/simurg/core/errors.py` — do not edit by hand. | `E_OUTPUT_EXISTS` | error | Outputs already exist; use --force or --resume | CLI-010 | | `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_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 | +| `E_PGX_VCF_UNSUITABLE` | error | SNV VCF unsuitable for PGx (squashing/locus filters/foreign) | PGX-011 | | `E_QC_CONFIG_EXT` | error | QC config file must be .yaml or .yml | QC-021 | | `E_QC_CONFIG_MISSING` | error | QC config file does not exist | QC-021 | | `E_QC_CONFIG_PARSE` | error | QC config file is not valid YAML | QC-021 | diff --git a/docs/validation/V2_pgx_star.md b/docs/validation/V2_pgx_star.md new file mode 100644 index 0000000..2522d96 --- /dev/null +++ b/docs/validation/V2_pgx_star.md @@ -0,0 +1,58 @@ +# V2 — PGx DB and star-allele calling (PGX-001 … PGX-021) + +Status: **implemented for CPIC star-allele genes, interim validation** (2026-09-29). + +* `scripts/pgx/fetch_sources.sh` → CPIC API tables (CC0, data version 2026-08-03). +* `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). +* Unit tests: `tests/unit/test_pgx_db.py`, `tests/unit/test_pgx_star.py`. + +## Data and truth + +* Arrays: the GSA v4 PGx demo (36 arrays = 30 CEPH/1000G individuals, `scripts/pgx/demo36_identity.tsv`), SNV VCFs + and PGx CN VCFs from SimurgArray. +* Star-allele truth: GeT-RM consolidated PGx/HLA table (CDC, 2025 V4) — 16 of the 30 individuals + (`scripts/pgx/getrm_truth.py`). +* Probe mask (PGX-010): array vs 1000G 30x (NYGC phased panel, fetched per gene over HTTPS) for the **other 14 + individuals** (`scripts/pgx/derive_masks.py`), so the mask is independent of the GeT-RM validation set: 8 of 658 + assessed records (23 probe names) below 95 % concordance — e.g. CYP2C19 rs367543002 and CYP2D6 rs1135836 called + homozygous ALT in every sample, the UGT1A1 `C>CATAT` TA-repeat probe. Regenerate with the command in the script + header (the mask lists vendor probe names and is not committed). + +## Result (16 individuals; replicate arrays counted once) + +| Gene | Top solution = GeT-RM | Truth among array-equivalent feasible solutions | +|---|---|---| +| CYP2C9 | **15/15** | 15/15 | +| CYP2C19 | **14/15** | 15/15 | +| CYP2D6 | **12/15** | 14/15 | +| CYP3A5 | 15/15 | 15/15 | +| TPMT | 15/16 | 15/16 | +| CYP2B6 | 12/14 | 12/14 | +| CYP4F2 | 8/9 | 8/9 | +| NAT2 | 8/9 | 8/9 | +| NUDT15 | 1/1 | 1/1 | +| **total** | **100/109 (91.7 %)** | **103/109 (94.5 %)** | + +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`. + +Not in the total (GeT-RM uses pre-2023 names the CPIC definitions no longer carry): UGT1A1 (GeT-RM `*60`; +0/9 top) and SLCO1B1 (`*1B/*17/*21` vs CPIC `*37/*15/*20`; 9/15 top). + +Remaining discordances: CYP2D6 `*3/(*68)+*4` (NA12878; *68 hybrids not modelled, *3 frameshift probe), CYP2D6 *2 vs +*1 where the array cannot separate them (prior picks *1), CYP2C19 *8, TPMT *21, CYP2B6 *27, CYP4F2 *2 — rare +alleles whose defining variants are absent or masked on this array. + +## Error model (K-029) + +After masking, 0 of 9 740 array calls at DB variants disagreed with 1000G 30x in any GenCall-score bin (0.15–1.0), +so the per-call error is 1e-3 (GS ≥ 0.15) and errors are handled by the probe mask. The first model, +ε = 0.5·(1 − GS)², let correct low-GS calls lose to the prior (CYP2C9 7/15 → 15/15 after the change). + +## Caveats + +* 16 individuals, one array type, European (CEPH) only; masks from 14 individuals (spec: ≥ 100 samples). +* 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). diff --git a/docs/validation/V2_pgx_star_getrm16.json b/docs/validation/V2_pgx_star_getrm16.json new file mode 100644 index 0000000..00899e5 --- /dev/null +++ b/docs/validation/V2_pgx_star_getrm16.json @@ -0,0 +1,240 @@ +{ + "individuals": 30, + "genes": { + "CYP2B6": { + "n": 14, + "top": 12, + "equivalent": 12, + "not_comparable": 1 + }, + "CYP2C19": { + "n": 15, + "top": 14, + "equivalent": 15 + }, + "CYP2C9": { + "n": 15, + "top": 15, + "equivalent": 15 + }, + "CYP2D6": { + "n": 15, + "top": 12, + "equivalent": 14 + }, + "CYP3A5": { + "n": 15, + "top": 15, + "equivalent": 15 + }, + "CYP4F2": { + "n": 9, + "top": 8, + "equivalent": 8 + }, + "NAT2": { + "n": 9, + "top": 8, + "equivalent": 8 + }, + "NUDT15": { + "n": 1, + "top": 1, + "equivalent": 1 + }, + "SLCO1B1": { + "n": 15, + "top": 9, + "equivalent": 9 + }, + "TPMT": { + "n": 16, + "top": 15, + "equivalent": 15 + }, + "UGT1A1": { + "n": 9, + "top": 0, + "equivalent": 1 + } + }, + "legacy_nomenclature_genes": [ + "SLCO1B1", + "UGT1A1" + ], + "total": { + "n": 109, + "top": 100, + "equivalent": 103 + }, + "discordant": [ + { + "individual": "NA12156", + "gene": "CYP4F2", + "truth": "(*2)/*3", + "call": "*1/*3", + "within_equivalents": false + }, + { + "individual": "NA12156", + "gene": "SLCO1B1", + "truth": "*1/*21", + "call": "*1/*37", + "within_equivalents": false + }, + { + "individual": "NA12156", + "gene": "UGT1A1", + "truth": "*60/*28", + "call": "*28/*80", + "within_equivalents": false + }, + { + "individual": "NA12044", + "gene": "TPMT", + "truth": "*1/*21", + "call": "*1/*1", + "within_equivalents": false + }, + { + "individual": "NA07019", + "gene": "UGT1A1", + "truth": "*1/*60", + "call": "*1/*1", + "within_equivalents": false + }, + { + "individual": "NA06991", + "gene": "UGT1A1", + "truth": "*60/*60", + "call": "*1/*1", + "within_equivalents": false + }, + { + "individual": "NA12878", + "gene": "CYP2D6", + "truth": "*3/(*68)+*4", + "call": "*4x2/*1;*4x2/*82", + "within_equivalents": false + }, + { + "individual": "NA12878", + "gene": "UGT1A1", + "truth": "*60/*28", + "call": "*28/*80", + "within_equivalents": false + }, + { + "individual": "NA11993", + "gene": "SLCO1B1", + "truth": "*1/*17", + "call": "*1/*15", + "within_equivalents": false + }, + { + "individual": "NA11993", + "gene": "UGT1A1", + "truth": "(*28 + *60)/*60", + "call": "*28/*80", + "within_equivalents": false + }, + { + "individual": "NA12003", + "gene": "UGT1A1", + "truth": "*60/*28", + "call": "*28/*80", + "within_equivalents": false + }, + { + "individual": "NA07029", + "gene": "CYP2B6", + "truth": "*6/(*27)", + "call": "*1/*6", + "within_equivalents": false + }, + { + "individual": "NA07029", + "gene": "CYP2C19", + "truth": "*8/*17", + "call": "*1/*17", + "within_equivalents": true + }, + { + "individual": "NA07029", + "gene": "UGT1A1", + "truth": "*1/*1", + "call": "*1/*28", + "within_equivalents": true + }, + { + "individual": "NA12892", + "gene": "CYP2B6", + "truth": "*1/*1", + "call": "*1/*5", + "within_equivalents": false + }, + { + "individual": "NA11832", + "gene": "NAT2", + "truth": "*5/(*12)", + "call": "*1/*5", + "within_equivalents": false + }, + { + "individual": "NA11832", + "gene": "SLCO1B1", + "truth": "*1/*1", + "call": "*14/*20", + "within_equivalents": false + }, + { + "individual": "NA11832", + "gene": "UGT1A1", + "truth": "(*28 + *60)/*60", + "call": "*28/*80", + "within_equivalents": false + }, + { + "individual": "NA10838", + "gene": "CYP2D6", + "truth": "*2/*4", + "call": "*1/*4;*2/*4", + "within_equivalents": true + }, + { + "individual": "NA10838", + "gene": "SLCO1B1", + "truth": "*1/*1", + "call": "*1/*14", + "within_equivalents": false + }, + { + "individual": "NA12813", + "gene": "CYP2D6", + "truth": "*2/*4", + "call": "*1/*4;*2/*4", + "within_equivalents": true + }, + { + "individual": "NA12813", + "gene": "SLCO1B1", + "truth": "*1/*21", + "call": "*1/*20", + "within_equivalents": false + }, + { + "individual": "NA12813", + "gene": "UGT1A1", + "truth": "(*28 + *60)/(*28 + *60)", + "call": "*80/*80", + "within_equivalents": false + }, + { + "individual": "NA11881", + "gene": "SLCO1B1", + "truth": "*1/(*15)", + "call": "*14/*15", + "within_equivalents": false + } + ] +} \ No newline at end of file diff --git a/scripts/pgx/build_cn_regions.py b/scripts/pgx/build_cn_regions.py index 64f3cf0..87f78bc 100644 --- a/scripts/pgx/build_cn_regions.py +++ b/scripts/pgx/build_cn_regions.py @@ -49,12 +49,12 @@ def main() -> int: rows: list[tuple[str, int, int, str, str]] = [] def exon(gene: str, k: int) -> tuple[int, int]: - chrom, strand, s, e = tx[gene] + _, strand, s, e = tx[gene] i = k - 1 if strand == "+" else len(s) - k return s[i], e[i] def intron(gene: str, k: int) -> tuple[int, int]: - chrom, strand, s, e = tx[gene] + strand = tx[gene][1] a_, b_ = exon(gene, k), exon(gene, k + 1) return (a_[1], b_[0]) if strand == "+" else (b_[1], a_[0]) diff --git a/scripts/pgx/derive_masks.py b/scripts/pgx/derive_masks.py new file mode 100644 index 0000000..0fbaf1e --- /dev/null +++ b/scripts/pgx/derive_masks.py @@ -0,0 +1,157 @@ +#!/usr/bin/env python3 +"""Derive a PGx probe mask from array-vs-WGS concordance (spec PGX-010) for the individuals of a GSA PGx batch. + + python scripts/pgx/derive_masks.py --db DB.sqlite --vcfs VCF_DIR --identity demo36_identity.tsv \ + --exclude-individuals getrm.txt --manifest NAME --out MASK.tsv [--min-concordance 0.95] + +Truth: the 1000 Genomes 30x phased panel (NYGC, 3202 samples), fetched per gene window over HTTPS with pysam (MIT; +a script-only dependency: ``uv pip install pysam``). For every PGx-DB defining variant the array record(s) at that +site are compared with the truth alt dosage; records whose concordance is below the threshold (≥ 5 comparisons) +are masked by their probe names. Individuals listed in ``--exclude-individuals`` (the validation set) are not used, +so the mask stays independent of the star-allele validation. The spec asks for ≥ 100 samples; this batch has fewer +(recorded as K-029). +""" + +from __future__ import annotations + +import argparse +import csv +import gzip +import sqlite3 +from collections import defaultdict +from pathlib import Path + +PANEL = ( + "https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/" + "20220422_3202_phased_SNV_INDEL_SV/1kGP_high_coverage_Illumina.{c}.filtered.SNV_INDEL_SV_phased_panel.vcf.gz" +) + + +def main() -> int: + import pysam + + ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("--db", required=True, type=Path) + ap.add_argument("--vcfs", required=True, type=Path) + ap.add_argument("--identity", required=True, type=Path) + ap.add_argument("--exclude-individuals", type=Path) + ap.add_argument("--manifest", required=True) + ap.add_argument("--min-concordance", type=float, default=0.95) + ap.add_argument("--out", required=True, type=Path) + a = ap.parse_args() + ident = {} + with a.identity.open() as fh: + for r in csv.DictReader((ln for ln in fh if not ln.startswith("#")), delimiter="\t"): + ident[r["sample_id"]] = r["individual"] + excl = set(a.exclude_individuals.read_text().split()) if a.exclude_individuals else set() + use = {s: i for s, i in ident.items() if i not in excl} + con = sqlite3.connect(a.db) + # the highest population frequency among the alleles that carry each variant (CPIC frequencies) + vfreq: dict[tuple[str, int, str, str], float] = defaultdict(float) + for c, p, r, al, f in con.execute( + "SELECT v.chrom, v.pos, v.ref, v.alt, MAX(f.freq) FROM allele_variants v JOIN frequencies f " + "ON f.gene = v.gene AND f.allele = v.allele GROUP BY v.chrom, v.pos, v.ref, v.alt" + ): + vfreq[(c, p, r, al)] = float(f or 0.0) + keys = sorted({(c, p, r, al, g) for g, c, p, r, al in con.execute( + "SELECT gene, chrom, pos, ref, alt FROM allele_variants")}) # fmt: skip + # array dosages per (chrom, pos, ref, alt) and record ID + wanted = defaultdict(set) + for c, p, *_ in keys: + wanted[c].add(p) + arr: dict[tuple[str, int, str, str], dict[str, dict[str, int]]] = defaultdict(lambda: defaultdict(dict)) + gsd: dict[tuple[str, int, str, str, str, str], float] = {} + for sid in sorted(use): + path = a.vcfs / f"{sid}.snv.vcf.gz" + if not path.exists(): + continue + with gzip.open(path, "rt") as fh: + for ln in fh: + if ln[0] == "#": + continue + c, p, rid, ref, alts = ln.split("\t", 5)[:5] + if int(p) not in wanted.get(c, ()): + continue + sample = ln.rstrip("\n").split("\t")[9].split(":") + gt = sample[0].replace("|", "/").split("/") + if "." in gt: + continue + gs_v = float(sample[1]) if len(sample) > 1 and sample[1] not in (".", "") else float("nan") + for ai, alt in enumerate(alts.split(","), start=1): + arr[(c, int(p), ref, alt)][rid][sid] = sum(int(x) == ai for x in gt) + gsd[(c, int(p), ref, alt, rid, sid)] = gs_v + # truth per gene window + truth: dict[tuple[str, int, str, str], dict[str, int]] = {} + by_gene = defaultdict(list) + for c, p, r, al, g in keys: + by_gene[(g, c)].append((p, r, al)) + inds = sorted(set(use.values())) + for (g, c), vs in sorted(by_gene.items()): + core = c.removeprefix("chr") + chrom = "chrM" if core in ("MT", "M") else f"chr{core}" + if chrom == "chrM": + continue + try: + vf = pysam.VariantFile(PANEL.format(c=chrom)) + except (OSError, ValueError): + print(f"{g} {chrom}: panel file not available, skipped") + continue + samples = [s for s in inds if s in vf.header.samples] + vf.subset_samples(samples) + lo, hi = min(p for p, _, _ in vs) - 10, max(p + len(r) for p, r, _ in vs) + 10 + want = {(p, r, al) for p, r, al in vs} + seen = set() + for rec in vf.fetch(chrom, lo, hi): + for ai, alt in enumerate(rec.alts or (), start=1): + k = (rec.pos, rec.ref, alt) + if k in want: + truth[(c, *k)] = {s: sum(x == ai for x in rec.samples[s]["GT"] if x is not None) for s in samples} + seen.add(k) + for p, r, al in want - seen: + # absent from the (filtered) panel: homozygous reference only for a rare variant; a common one was + # most likely filtered out (e.g. CYP2B6 rs2279343 next to CYP2B7) and is not assessed + if vfreq.get((c, p, r, al), 0.0) < 0.05: + truth[(c, p, r, al)] = dict.fromkeys(samples, 0) + print(f"{g} {chrom}: {len(want)} variant(s), {len(seen)} in the panel") + rows = [] + stats = [0, 0] + bins = [0.0, 0.15, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 1.01] + calib = [[0, 0] for _ in bins[:-1]] # [discordant, n] per GS bin over unmasked records (error-model calibration) + for k, recs in sorted(arr.items()): + t = truth.get(k) + if t is None: + continue + for rid, dos in sorted(recs.items()): + pairs = [(d, t[use[s]]) for s, d in dos.items() if use[s] in t] + if len(pairs) < 5: + continue + conc = sum(x == y for x, y in pairs) / len(pairs) + stats[1] += 1 + if conc >= a.min_concordance: + for s_, d in dos.items(): + if use[s_] not in t: + continue + gv = gsd.get((*k, rid, s_), float("nan")) + for bi in range(len(bins) - 1): + if bins[bi] <= gv < bins[bi + 1]: + calib[bi][0] += d != t[use[s_]] + calib[bi][1] += 1 + if conc < a.min_concordance: + stats[0] += 1 + for probe in rid.split(","): + rows.append((probe, f"concordance {conc:.2f} vs 1000G 30x at {k[0]}:{k[1]} {k[2]}>{k[3]} " + f"(n={len(pairs)})", a.manifest)) # fmt: skip + with a.out.open("w") as fh: + fh.write(f"# PGx probe mask for {a.manifest} (PGX-010): array vs 1000G 30x, {len(inds)} individuals, " + f"threshold {a.min_concordance}\n") # fmt: skip + for r in sorted(set(rows)): + fh.write("\t".join(r) + "\n") + print(f"masked {stats[0]} of {stats[1]} assessed record(s); {len(set(rows))} probe name(s) -> {a.out}") + for bi, (dsc, n) in enumerate(calib): + if n: + print(f"GS [{bins[bi]:.2f}, {bins[bi + 1]:.2f}): discordance {dsc}/{n} = {dsc / n:.4f} (unmasked records)") + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/scripts/pgx/getrm_truth.py b/scripts/pgx/getrm_truth.py new file mode 100644 index 0000000..74c5e10 --- /dev/null +++ b/scripts/pgx/getrm_truth.py @@ -0,0 +1,51 @@ +#!/usr/bin/env python3 +"""Extract GeT-RM consensus diplotypes (CDC Genetic Testing Reference Material program, consolidated PGx/HLA table, +public) for the individuals of an identity table. + + python scripts/pgx/getrm_truth.py --xlsx Consolidated_PGx-HLA_table_1-22-25-V4.xlsx \ + --identity scripts/pgx/demo36_identity.tsv --out getrm_truth.tsv +Writes ``individual, gene, diplotype`` (the table's text, trimmed) for the star-allele genes. +""" + +from __future__ import annotations + +import argparse +import csv +from pathlib import Path + +GENES = ("CYP2B6", "CYP2C9", "CYP2C19", "CYP2D6", "CYP3A4", "CYP3A5", "CYP4F2", "NAT2", "NUDT15", "SLCO1B1", "TPMT", + "UGT1A1", "DPYD", "CYP2A6", "CYP1A2", "CYP2C8", "GSTM1", "GSTT1", "UGT2B17") # fmt: skip + + +def main() -> int: + import openpyxl + + ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("--xlsx", required=True, type=Path) + ap.add_argument("--identity", required=True, type=Path) + ap.add_argument("--out", required=True, type=Path) + a = ap.parse_args() + with a.identity.open() as fh: + inds = {r["individual"] for r in csv.DictReader((ln for ln in fh if not ln.startswith("#")), delimiter="\t")} + wb = openpyxl.load_workbook(a.xlsx, read_only=True) + 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} + out = [] + for r in rows[2:]: + if not r or not r[1] or str(r[1]).strip() not in inds: + continue + for g, i in col.items(): + v = r[i] + if v is not None and str(v).strip(): + out.append((str(r[1]).strip(), g, " ".join(str(v).split()))) + with a.out.open("w") as fh: + fh.write("individual\tgene\tdiplotype\n") + for t in sorted(out): + fh.write("\t".join(t) + "\n") + print(f"{len({t[0] for t in out})} individual(s), {len(out)} gene call(s) -> {a.out}") + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/scripts/pgx/validate_cn_truth.py b/scripts/pgx/validate_cn_truth.py index 474049b..c9c1813 100644 --- a/scripts/pgx/validate_cn_truth.py +++ b/scripts/pgx/validate_cn_truth.py @@ -59,10 +59,12 @@ def main() -> int: tot[1] += len(pairs) res[region] = {"record": rid, "concordant": ok, "n": len(pairs), "rate": round(ok / len(pairs), 4) if pairs else None, - "truth_to_call_mismatches": dict(collections.Counter(f"{t}->{c}" for t, c in pairs if t != c))} # fmt: skip + "truth_to_call_mismatches": dict(collections.Counter( + f"{t}->{c}" for t, c in pairs if t != c))} # fmt: skip summary = {"individuals": len(calls), "regions": res, "overall": {"concordant": tot[0], "n": tot[1], "rate": round(tot[0] / tot[1], 4) if tot[1] else None}, - "gstt1_cn_distribution": dict(sorted(collections.Counter(c.get("GSTT1", ".") for c in calls.values()).items()))} # fmt: skip + "gstt1_cn_distribution": dict(sorted(collections.Counter( + c.get("GSTT1", ".") for c in calls.values()).items()))} # fmt: skip a.out.write_text(json.dumps(summary, indent=1)) print(json.dumps(summary, indent=1)) return 0 diff --git a/scripts/pgx/validate_star_getrm.py b/scripts/pgx/validate_star_getrm.py new file mode 100644 index 0000000..8821719 --- /dev/null +++ b/scripts/pgx/validate_star_getrm.py @@ -0,0 +1,118 @@ +#!/usr/bin/env python3 +"""Compare `simurg pgx star-allele call` with GeT-RM consensus diplotypes (PGX-032). + + python scripts/pgx/validate_star_getrm.py --calls STAR_DIR/pgx_calls.json --truth getrm_truth.tsv \ + --identity scripts/pgx/demo36_identity.tsv --out summary.json + +GeT-RM notation is normalised per allele: "*2(*35)" / "*2 (*35)" accept *2 or *35; a trailing sub-allele letter +("*1F") also matches the core allele ("*1") when the call has no letter; "(*68)+*4" tandems and "xN" duplications +reduce to the core allele (structural forms are reported separately); "(*28 + *60)" keeps the star alleles the DB +defines. "A/B or C/D" accepts either. UGT1A1 and SLCO1B1 are reported but left out of the total: GeT-RM uses pre-2023 +names (UGT1A1 *60, SLCO1B1 *1B/*17/*21) that the CPIC definitions no longer carry. A sample's gene counts as **top** +concordant when the top solution matches, and as **equivalent** +concordant when the truth is among the array-equivalent members of any feasible solution. Replicate arrays of one +individual are counted once (the first). +""" + +from __future__ import annotations + +import argparse +import csv +import json +import re +from collections import defaultdict +from pathlib import Path + + +def truth_options(text: str) -> list[set[str]] | None: + """Per haplotype, the set of acceptable core alleles; None when the text is not a diplotype.""" + t = text.replace(" ", "") + if "/" not in t: + return None + out = [] + for side in t.split("/", 1): + opts = set(re.findall(r"\*\d+[A-Z]?", side)) + if not opts: + return None + cores = {re.sub(r"x\d+|xN", "", o) for o in opts} + cores |= {re.sub(r"[A-Z]$", "", o) for o in cores} + if "+" in side: # tandem / hybrid: the non-parenthesised last unit is the core allele + last = re.findall(r"\*\d+[A-Z]?", side.split("+")[-1]) + cores |= set(last) + out.append(cores) + return out + + +LEGACY_NOMENCLATURE = {"UGT1A1", "SLCO1B1"} # GeT-RM uses pre-2023 names (*60, SLCO1B1 *1B/*17/*21) not in CPIC + + +def truth_alternatives(text: str) -> list[list[set[str]]]: + """GeT-RM "A/B or C/D" → one option list per alternative diplotype.""" + alts = [truth_options(x) for x in re.split(r"\bor\b", text)] + return [x for x in alts if x is not None] + + +def matches(call: tuple[str, str], truth: list[list[set[str]]]) -> bool: + a, b = (re.sub(r"x\d+$", "", x) for x in call) + return any((a in t[0] and b in t[1]) or (a in t[1] and b in t[0]) for t in truth) + + +def main() -> int: + ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("--calls", required=True, type=Path) + ap.add_argument("--truth", required=True, type=Path) + ap.add_argument("--identity", required=True, type=Path) + ap.add_argument("--out", required=True, type=Path) + a = ap.parse_args() + with a.identity.open() as fh: + ident = {r["sample_id"]: r["individual"] for r in csv.DictReader((ln for ln in fh if not ln.startswith("#")), + delimiter="\t")} # fmt: skip + truth: dict[tuple[str, str], str] = {} + with a.truth.open() as fh: + for r in csv.DictReader(fh, delimiter="\t"): + truth[(r["individual"], r["gene"])] = r["diplotype"] + calls = json.loads(a.calls.read_text()) + per_gene: dict[str, dict[str, int]] = defaultdict(lambda: defaultdict(int)) + detail = [] + seen = set() + for s in calls["samples"]: + ind = ident.get(s["sample_id"]) + if ind is None or ind in seen: + continue + seen.add(ind) + for g in s["genes"]: + t = truth.get((ind, g["gene"])) + if t is None: + continue + opts = truth_alternatives(t) + st = per_gene[g["gene"]] + if not opts: + st["not_comparable"] += 1 + continue + st["n"] += 1 + top = g["candidates"][0]["genotype"].split("/") if g["candidates"] else ["?", "?"] + ok_top = len(top) == 2 and matches((top[0], top[1]), opts) + feas = [c for c in g["candidates"] if c["genotype"] in g["feasible_solutions"]] + ok_eq = ok_top or any(matches(tuple(e.split("/")[:2]), opts) for c in feas for e in c["equivalents"] + if len(e.split("/")) >= 2) # fmt: skip + st["top"] += ok_top + st["equivalent"] += ok_eq + if not ok_top: + detail.append({"individual": ind, "gene": g["gene"], "truth": t, "call": g["genotype"], + "within_equivalents": ok_eq}) # fmt: skip + comparable = {k: v for k, v in per_gene.items() if k not in LEGACY_NOMENCLATURE} + summary = {"individuals": len(seen), "genes": {k: dict(v) for k, v in sorted(per_gene.items())}, + "legacy_nomenclature_genes": sorted(LEGACY_NOMENCLATURE), + "total": {k: sum(v.get(k, 0) for v in comparable.values()) for k in ("n", "top", "equivalent")}, + "discordant": detail} # fmt: skip + a.out.write_text(json.dumps(summary, indent=1)) + for g, v in sorted(per_gene.items()): + print(f"{g:8s} top {v['top']}/{v['n']} within equivalents {v['equivalent']}/{v['n']}" + + (f" (not comparable {v['not_comparable']})" if v.get("not_comparable") else "") + + (" [legacy nomenclature]" if g in LEGACY_NOMENCLATURE else "")) # fmt: skip + print("total (excluding legacy-nomenclature genes)", summary["total"]) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/src/simurg/cli/main.py b/src/simurg/cli/main.py index e8bf0f3..394192f 100644 --- a/src/simurg/cli/main.py +++ b/src/simurg/cli/main.py @@ -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).")], diff --git a/src/simurg/core/errors.py b/src/simurg/core/errors.py index 41ad17b..6758677 100644 --- a/src/simurg/core/errors.py +++ b/src/simurg/core/errors.py @@ -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"), diff --git a/src/simurg/pgx/star.py b/src/simurg/pgx/star.py new file mode 100644 index 0000000..fb110e8 --- /dev/null +++ b/src/simurg/pgx/star.py @@ -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 diff --git a/src/simurg/pipeline/pgx_star.py b/src/simurg/pipeline/pgx_star.py new file mode 100644 index 0000000..fc0b1bb --- /dev/null +++ b/src/simurg/pipeline/pgx_star.py @@ -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 diff --git a/tests/unit/test_pgx_star.py b/tests/unit/test_pgx_star.py new file mode 100644 index 0000000..1cf6ec7 --- /dev/null +++ b/tests/unit/test_pgx_star.py @@ -0,0 +1,68 @@ +"""Star-allele calling (spec PGX-012 … PGX-017).""" + +from __future__ import annotations + +import numpy as np + +from simurg.pgx.star import CALLED, MASKED, NOT_COVERED, GeneModel, Observations, call_gene, gs_error +from simurg.pipeline.pgx_star import observations + + +def _model() -> GeneModel: + # variants: v0 defines *2, v1 defines *3, v2 is carried by *4 and *3 wildcard, v3 not on the array + variants = [("1", 100, "A", "G"), ("1", 200, "C", "T"), ("1", 300, "G", "A"), ("1", 400, "T", "C")] + carries = np.array([[0, 0, 0, 0], # *1 + [1, 0, 0, 0], # *2 + [0, 1, -1, 0], # *3 (v2 wildcard) + [0, 0, 1, 0], # *4 + [0, 0, 0, 1]], # *5 (only v3, not covered) + dtype=np.int8) # fmt: skip + return GeneModel("G", variants, [(f"v{i}", None) for i in range(4)], ["*1", "*2", "*3", "*4", "*5"], carries, + np.array([0.7, 0.1, 0.1, 0.05, 0.05]), deletion_allele="*0", reference="*1") # fmt: skip + + +def _obs(dosage: list[int], state: list[int] | None = None) -> Observations: + st = np.array(state or [CALLED, CALLED, CALLED, NOT_COVERED], dtype=np.int8) + return Observations(st, np.array(dosage, dtype=np.int8), np.full(4, 0.2), np.full(4, np.nan)) + + +def test_error_model() -> None: # PGX-012 (calibrated, K-029) + assert gs_error(np.array([0.16, 0.9, 0.1, np.nan])).tolist() == [1e-3, 1e-3, 0.05, 0.05] + + +def test_simple_diplotypes() -> None: # PGX-013/014 + m = _model() + sols, feas = call_gene(m, _obs([1, 0, 0, -1])) + assert sols[0].genotype == "*1/*2" and sols[0].score > 0.9 and feas == [sols[0]] + sols, _ = call_gene(m, _obs([0, 1, 1, -1])) + assert sols[0].genotype == "*1/*3" # *3 with its wildcard at v2 explains the v2 het too (vs *3/*4) + sols, _ = call_gene(m, _obs([2, 0, 0, -1])) + assert sols[0].genotype == "*2/*2" + + +def test_array_equivalence() -> None: # PGX-017: *5 differs from *1 only at an uncovered variant + m = _model() + sols, _ = call_gene(m, _obs([0, 0, 0, -1])) + top = sols[0] + assert top.genotype == "*1/*1" and ("*1", "*5") in top.members and ("*5", "*5") in top.members + assert abs(sum(s.raw_score for s in sols) - 1) < 1e-9 + + +def test_copy_number_forms() -> None: # PGX-013 structural candidates + m = _model() + obs = Observations(np.array([CALLED, CALLED, CALLED, NOT_COVERED], np.int8), np.array([1, 0, 0, -1], np.int8), + np.full(4, 0.5), np.array([2 / 3, 0.0, 0.0, np.nan])) # fmt: skip + sols, _ = call_gene(m, obs, cn=3) + assert sols[0].genotype == "*2x2/*1" and sols[0].copy_number == 3 + obs.baf[0] = 1.0 + sols, _ = call_gene(m, obs, cn=1) + assert sols[0].genotype == "*2/*0" + + +def test_masked_record() -> None: # PGX-010 + m = _model() + recs = {("1", 100): [{"id": "p1,p1b", "ref": "A", "alts": ["G"], "gt": "1/1", "gs": "0.3", "baf": "1"}], + ("1", 200): [{"id": "p2", "ref": "C", "alts": ["T"], "gt": "0/1", "gs": "0.3", "baf": "0.5"}]} # fmt: skip + 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