kemalandClaude Opus 5.5
83654ddaef
- scripts/pgx/wgs_star_truth.py: star alleles from the 1000G 30x phased panel with the DB definitions (wildcard-aware haplotype assignment; checked against GeT-RM on the demo CEPH individuals). - validate_star_getrm.py: --all-genes, --individuals, --coverage (projection onto what the array can see); replicate_concordance.py. - 120 individuals x 3 GSA v1 arrays: array-resolvable 1637/1710 (95.7 %), top 1223/1800 (chip lacks many defining variants), GeT-RM subset 28/35, replicates 99.2 % without the failed array (V2_pgx_star.md, K-029). Co-Authored-By: Claude Opus 5.5 <[email protected]> Claude-Session: https://claude.ai/code/session_01814GKVzZVL3FuF3DeiZXhv
155 lines
10 KiB
Markdown
155 lines
10 KiB
Markdown
# 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).
|
||
* `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
|
||
|
||
* 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 | **13/15** | **15/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 |
|
||
| VKORC1 (rs9923231, c.-1639G>A) | 13/15 | 13/15 |
|
||
| GSTM1 (CN: DEL/DEL vs present) | 9/9 | 9/9 |
|
||
| GSTT1 (CN 0/1/2) | 7/9 | 7/9 |
|
||
| UGT2B17 (CN) | 2/2 | 2/2 |
|
||
| CYP3A4 (PharmVar) | **15/15** | 15/15 |
|
||
| CYP2A6 (PharmVar; *4 deletion from CN) | 7/9 | 8/9 |
|
||
| **total** | **154/168 (91.7 %)** | **158/168 (94.0 %)** |
|
||
|
||
PGX-032 target (≥ 90 % for CYP2D6/CYP2C19/CYP2C9 against truth): **42/45 = 93 %** (top solution), 45/45 within
|
||
the array-equivalent feasible solutions.
|
||
|
||
CYP2D6 structures (PGX-013) come from the three sub-region CNs (`simurg.pgx.star.cyp2d6_structure`): NA12878
|
||
`*3/*68+*4` and NA11832 `*1/*68+*4` are called exactly as GeT-RM reports them; one sample is called `*36+*10/*4`
|
||
(p5 = intron 2 = 3, exon 9 = 2; no GeT-RM truth). The two remaining top-solution misses are *1 vs *2, which this
|
||
array cannot separate (both within the feasible set).
|
||
Details and every discordant call: `V2_pgx_star_getrm16.json`.
|
||
|
||
Not in the total (GeT-RM uses older names or definitions than CPIC/PharmVar today): UGT1A1 (GeT-RM `*60`;
|
||
0/9 top), SLCO1B1 (`*1B/*17/*21` vs CPIC `*37/*15/*20`; 9/15), CYP1A2 (lettered `*1F`; 0/9) and CYP2C8 (PharmVar now
|
||
defines rs2071426 alone as `*19`; the array calls it concordantly with 1000G, GeT-RM reports those carriers as
|
||
`*1`; 8/15).
|
||
|
||
Remaining discordances: 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).
|
||
* 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.
|
||
* `simurg pgx coverage` (PGX-009/035) on GSA v4 PGx with the mask: 438 alleles callable, 68 array-equivalent to
|
||
another allele, 88 not covered, 4 masked. Against the reference's known limitations: CYP4F2 *2 not covered (reported
|
||
as *1) and *4 ≡ *3 — as the reference; UGT1A1 *28 callable (its C>CAT probe matches 1000G) and *37 masked (the
|
||
C>CATAT probe fails); CYP2D6 *40/*58 and NUDT15 *6/*9 are separable here (CPIC core definitions differ from the
|
||
reference's).
|
||
* Detail outputs (PGX-016 … PGX-020, PGX-027/028, PGX-031): `supporting_variants.tsv`, `unresolved_variants.tsv`,
|
||
`variants_not_covered.tsv`, per-allele supporting/unresolved sites for **both** alleles of every candidate,
|
||
`variant_change` vs the top solution (`distant_alternative` above 3 changes), `unexplained_variants`, full HGVS on
|
||
GRCh38 RefSeq accessions.
|
||
* All CPIC gene kinds are called: star alleles, G6PD (named haplotypes; males hemizygous from the VCF ploidy),
|
||
variant-based genes (ABCG2, CACNA1S, CFTR, DPYD, IFNL3, MT-RNR1, RYR1, VKORC1; genotype as `rs… REF/ALT` or
|
||
`Reference`, phenotype diplotype from the two lowest-activity carried variants), HLA-A/B tag variants (PGX-007,
|
||
positions resolved from the manifest's probes, flag `tag_caveat`), and GSTM1/GSTT1/UGT2B17 `*0` genotypes from
|
||
the targeted CN (PGX-006).
|
||
* PharmVar genes without CPIC definitions (CYP1A2, CYP2A6, CYP2A13, CYP2C8, CYP3A4; PharmVar 6.2.29): core alleles
|
||
from the GRCh38 haplotype VCFs, `*1` = reference, prior reference 0.8 (no population frequencies); genotype only
|
||
(no CPIC phenotype). The probe mask was re-derived with these variants (one more record: CYP2A6 19:40843845 C>A,
|
||
concordance 0.47, which had produced false *5 calls).
|
||
* ClinPGx clinical-annotation variants (CC BY-SA 4.0; levels 1A–2B, and level 3 for the PGX-008 parity genes) add
|
||
genotype-only variant-based genes — OPRM1, COMT, DRD2, ANKK1, MTHFR, F5, … (43 genes, 171 of 187 rsIDs resolved on
|
||
GSA v4 PGx from the manifest's probes). The DB now holds 70 genes (16 star-allele, G6PD, 2 HLA tag, 51
|
||
variant-based; spec parity list: 58). `pgx coverage`: 1225 alleles callable, 93 array-equivalent, 437 not
|
||
covered, 5 masked.
|
||
* `--guideline dpwg` (PGX-025): ClinPGx curates DPWG guidelines with the CPIC allele-function and phenotype mappings
|
||
for genes CPIC covers, so those genes reuse the CPIC tables (flag `dpwg_via_cpic_mapping`); DPWG-only genes
|
||
(CYP1A2, CYP3A4, F5) have no public machine-readable mapping and get phenotype `null` with the flag
|
||
`dpwg_mapping_not_available` (K-029).
|
||
|
||
## Independent validation: GSE165845 (GSA v1 C2, 120 1000 Genomes individuals × 3 arrays)
|
||
|
||
A second cohort, independent of the demo used for development (2026-09-29, all on the Pi, $0):
|
||
|
||
* **Arrays:** the 360 GSE165845 GTCs (GSA-MD v1 chips, GSA v1 C2 manifest; K-024) → SNV VCFs; no targeted CN probes
|
||
on this chip, so star alleles are called with `--no-cn`. The chip carries little PGx content: 292 probe-to-variant
|
||
links (GSA v4 PGx: 5 280); `pgx coverage`: 495 alleles callable, 136 array-equivalent, 1 129 not covered — e.g.
|
||
no CYP2D6 *2/*4/*10 probes (rs16947, rs3892097, rs1065852 are absent from the manifest).
|
||
* **Truth:** star alleles derived from the 1000G 30x phased panel with the same DB definitions
|
||
(`scripts/pgx/wgs_star_truth.py`; SNV/indel level — CYP2D6 *5, xN and hybrids invisible; CYP2B6 left out because
|
||
its common core variant rs2279343 is filtered from the panel). The derivation was checked against GeT-RM on the
|
||
demo's CEPH individuals (CYP2D6 14/15, CYP2C9/CYP3A4/CYP3A5 15/15, TPMT 16/16; the rest are nomenclature changes).
|
||
* **Masks:** 2-fold — individuals split in halves; each half is called with the mask derived from the other half
|
||
(A → 0 probes masked, B → 1). After masking the calls agree with 1000G at 99 % or better in every GS bin.
|
||
* **Metrics:** *top* = the top solution equals the WGS diplotype; *array-resolvable* = both diplotypes projected
|
||
onto the alleles this array can tell apart (`pgx coverage`: callable → itself, equivalent_to:X → X, not covered →
|
||
the reference's class) — the fair measure for a chip that lacks a gene's defining variants.
|
||
|
||
| Gene | Top = WGS | Truth within feasible equivalents | Array-resolvable |
|
||
|---|---|---|---|
|
||
| CYP1A2 | 18/120 | 72/120 | 120/120 |
|
||
| CYP2A13 | 97/120 | 119/120 | 117/120 |
|
||
| CYP2A6 | 73/120 | 112/120 | 117/120 |
|
||
| CYP2C19 | 100/120 | 116/120 | 110/120 |
|
||
| CYP2C8 | 114/120 | 120/120 | 115/120 |
|
||
| CYP2C9 | 95/120 | 108/120 | 116/120 |
|
||
| CYP2D6 (SNV level) | 32/120 | 54/120 | 111/120 |
|
||
| CYP3A4 | 112/120 | 117/120 | 117/120 |
|
||
| CYP3A5 | 120/120 | 120/120 | 120/120 |
|
||
| CYP4F2 | 43/120 | 86/120 | 96/120 |
|
||
| NAT2 | 95/120 | 105/120 | 120/120 |
|
||
| NUDT15 | 120/120 | 120/120 | 120/120 |
|
||
| SLCO1B1 | 31/120 | 87/120 | 109/120 |
|
||
| TPMT | 107/120 | 115/120 | 102/103 |
|
||
| UGT1A1 | 66/120 | 120/120 | 47/47 |
|
||
| **total** | **1223/1800 (67.9 %)** | | **1637/1710 (95.7 %)** |
|
||
|
||
* Top-solution concordance is low where the chip lacks the defining variants: the call is then the most probable
|
||
member of an equivalence class the array cannot split (CYP2D6 *1/*2/*10, SLCO1B1 *1/*37, CYP1A2), so the prior
|
||
decides. Projected onto what the array can see, 95.7 % agree (half A 815/854, half B 822/856).
|
||
* CYP4F2 (96/120) is the weakest: the array genotypes at its covered sites match 1000G exactly (rs3093200 180/180),
|
||
so the residue is definitional — *5 and *23 share the only covered variant, and the projection cannot separate
|
||
alleles whose other defining variants are off-array.
|
||
* **GeT-RM** (7 of the 120 individuals, older consensus): 28/35 top, array-resolvable 29/30.
|
||
* **Replicates:** the top genotype of all 70 genes agrees between replicate arrays of one individual in 24 916 of
|
||
25 200 pairs (98.9 %); 87 of the 284 discordant pairs come from NA20362's failed replicate (call rate 0.76, K-024);
|
||
without it 99.2 %.
|
||
|
||
Files: `docs/validation/V2_pgx_star_gse165845.json`; scripts `wgs_star_truth.py`, `validate_star_getrm.py`
|
||
(`--all-genes --individuals --coverage`), `replicate_concordance.py`.
|