From c3e0681b6e005333c92c453b6c02d234e9e1b57c Mon Sep 17 00:00:00 2001 From: Kemal Yaylali Date: Tue, 29 Sep 2026 13:16:19 +0100 Subject: [PATCH] PGx CN: pgx copy-number train|call (CLI-031/032, PCN-001..017, IN-050), region BED v1, demo identities (K-028) - simurg/pgx/cn.py: per-region CN mixture, mu_c = delta + gamma*log2((c+1)/3), free lower CN-0 component, Hardy-Weinberg-constrained weights over haplotype copy numbers, EM from a grid of label starts. - pipeline/pgx_cn.py: model file (.simurg-pgxcn), batch rules (PCN-004), per-sample CN VCF + tbi, batch and sample summaries. - The GSA v4 PGx demo arrays are 30 CEPH/1000G individuals (scripts/pgx/demo36_identity.tsv); against NYGC 30x SV RD_CN 224/240 region calls agree (GSTM1, UGT2B17 30/30), in-sample (docs/validation/V2_pgx_cn.md). --- THIRD_PARTY_NOTICES.md | 2 +- docs/known_unknowns.md | 1 + docs/messages.md | 9 + docs/validation/V2_pgx_cn.md | 48 ++ docs/validation/V2_pgx_cn_demo36.json | 89 +++ scripts/pgx/build_cn_regions.py | 85 +++ scripts/pgx/demo36_identity.tsv | 39 ++ scripts/pgx/validate_cn_truth.py | 72 +++ src/simurg/cli/main.py | 72 +++ src/simurg/core/errors.py | 20 + src/simurg/pgx/cn.py | 204 +++++++ .../pgx/data/pgx_cn_regions.GRCh38.v1.bed | 11 + src/simurg/pipeline/pgx_cn.py | 510 ++++++++++++++++++ tests/unit/test_pgx_cn.py | 145 +++++ 14 files changed, 1306 insertions(+), 1 deletion(-) create mode 100644 docs/validation/V2_pgx_cn.md create mode 100644 docs/validation/V2_pgx_cn_demo36.json create mode 100644 scripts/pgx/build_cn_regions.py create mode 100644 scripts/pgx/demo36_identity.tsv create mode 100644 scripts/pgx/validate_cn_truth.py create mode 100644 src/simurg/pgx/cn.py create mode 100644 src/simurg/pgx/data/pgx_cn_regions.GRCh38.v1.bed create mode 100644 src/simurg/pipeline/pgx_cn.py create mode 100644 tests/unit/test_pgx_cn.py diff --git a/THIRD_PARTY_NOTICES.md b/THIRD_PARTY_NOTICES.md index 89fe351..b9647c8 100644 --- a/THIRD_PARTY_NOTICES.md +++ b/THIRD_PARTY_NOTICES.md @@ -13,7 +13,7 @@ used only as external subprocess oracles and are not copied or linked. | [Illumina GTCtoVCF](https://github.com/Illumina/GTCtoVCF) — Copyright Illumina, Inc. | Apache-2.0 | Indel source-sequence handling ported to `simurg/harmonize/sites.py` (design-strand orientation, left shift, IUPAC flank matching, deletion/insertion context score); multi-probe consistency idea in `simurg/harmonize/combine.py`; test oracle (run under Python 3 with two local compatibility patches, never shipped) | | [Picard `picard.arrays`](https://github.com/broadinstitute/picard) | MIT | Second format reference | | [PyO3](https://github.com/PyO3/pyo3), [rust-numpy](https://github.com/PyO3/rust-numpy), [ndarray](https://github.com/rust-ndarray/ndarray) and their Rust dependencies | MIT / Apache-2.0 / BSD-2-Clause | Compiled into the optional `simurg_core` extension (`rust/simurg-core`); list with `cargo tree` | -| [UCSC Genome Browser](https://genome.ucsc.edu/) `cytoBand` tables (hg38, hg19) and `ncbiRefSeq` (RefSeq transcripts, NCBI) | Free for all uses (UCSC data terms; RefSeq is public domain) | `simurg/cn/data/cytoBand_GRCh3{7,8}.txt.gz` and the derived centromere/acrocentric mask BEDs are shipped; `scripts/build_annotation_bundle.sh` downloads cytoBand + ncbiRefSeq for gene lists | +| [UCSC Genome Browser](https://genome.ucsc.edu/) `cytoBand` tables (hg38, hg19) and `ncbiRefSeq` (RefSeq transcripts, NCBI) | Free for all uses (UCSC data terms; RefSeq is public domain) | `simurg/cn/data/cytoBand_GRCh3{7,8}.txt.gz` and the derived centromere/acrocentric mask BEDs are shipped; `scripts/build_annotation_bundle.sh` downloads cytoBand + ncbiRefSeq for gene lists; `simurg/pgx/data/pgx_cn_regions.GRCh38.v1.bed` is derived from ncbiRefSeq exon coordinates (`scripts/pgx/build_cn_regions.py`) | | [illuminaio](https://github.com/HenrikBengtsson/illuminaio) | GPL-2 | Documentation of IDAT field codes only; oracle via subprocess | ## BSD-2-Clause (BeadArrayFiles, verbatim from its README) diff --git a/docs/known_unknowns.md b/docs/known_unknowns.md index 626a2dc..a2ccb06 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-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. | | K-025 | 2026-09-28 | Scaling limits (P9) | With the Rust kernels a Graviton3 core genotypes 78–80 GSA samples/min. At > 8 workers a 125 MB/s gp3 volume caps the run at ~620 samples/min (GTC ≈ 13 MB each); on a RAM disk efficiency is 0.77 at 16 and 0.48 at 32 workers for 360 samples, explained by a fixed start-up cost (≈ 2.5 s + 0.2 s per worker). Linear scaling for large batches is inferred from that model, not measured. | `docs/validation/P9.md`; measure a ≥ 5000-sample batch when one is at hand. | diff --git a/docs/messages.md b/docs/messages.md index 804104f..fcd114c 100644 --- a/docs/messages.md +++ b/docs/messages.md @@ -36,11 +36,15 @@ Generated from `src/simurg/core/errors.py` — do not edit by hand. | `E_METH_NO_CONTROLS` | error | Missing control probes in manifest | MET-002 | | `E_METH_NO_IDATS` | error | No IDATs found | MET-004 | | `E_METH_PREP_FAILED` | error | Prep failed (sample) | MET-010 | +| `E_MODEL_CLUSTER_MISMATCH` | error | Input GTCs were called with another cluster file than the model | IN-050 | | `E_MODEL_MANIFEST_MISMATCH` | error | Model was built for a different manifest probe set | IN-051 | | `E_MODEL_VERSION` | error | Model schema is too old; rebuild it | IN-051 | | `E_NO_INPUT` | error | Neither an input folder nor a sample sheet was given | CLI-006 | | `E_NO_INPUT_FOLDER` | error | Sample sheet has no RootFolder and no folder was given | CLI-006 | | `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_TRAIN_TOO_SMALL` | error | Too few (passing) samples to train a PGx CN model | PCN-005 | | `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 | @@ -92,7 +96,12 @@ Generated from `src/simurg/core/errors.py` — do not edit by hand. | `W_METH_MISSING_CONTROL_TYPES` | warning | Control types missing from the manifest; metric set to NA | MET-009 | | `W_METH_MISSING_SAMPLE` | warning | Missing sample: IDATs absent; the sample fails, the run continues | MET-015 | | `W_METH_NO_FEATURES` | warning | No features, skipping PCA plot | MET-013 | +| `W_MODEL_CLUSTER_MISMATCH` | warning | Input GTCs were called with another cluster file than the model | IN-050 | | `W_NO_BGZIP` | warning | Uncompressed VCF written; downstream commands need bgzip + index | OUT-010 | +| `W_PGX_BATCH_BELOW_RECOMMENDED` | warning | Fewer passing samples than recommended | PCN-004 | +| `W_PGX_BATCH_COMPOSITION` | warning | Batch spans several lots/batches or > 14 days of scans | PCN-006 | +| `W_PGX_NO_ALT_CONTIG` | warning | FASTA lacks a contig carrying a PGx region; REF written as N | IN-062 | +| `W_PGX_REGION_FEW_PROBES` | warning | PGx CN region has fewer than 10 unmasked probes | PCN-008 | | `W_QC_LABEL_SANITIZED` | warning | Dataset label was sanitized | QC-042 | | `W_QC_NO_BAF` | warning | GTC has no B-allele frequencies; contamination not estimated | QC-064 | | `W_QC_PEDIGREE` | warning | Parent1/Parent2/Replicate value does not name a sample in the set | IN-018 | diff --git a/docs/validation/V2_pgx_cn.md b/docs/validation/V2_pgx_cn.md new file mode 100644 index 0000000..8b22b01 --- /dev/null +++ b/docs/validation/V2_pgx_cn.md @@ -0,0 +1,48 @@ +# V2 — Targeted PGx copy number (PCN-001 … PCN-017) + +Status: **implemented, interim validation** (2026-09-29). `simurg pgx copy-number train|call` (CLI-031/032), model +file `.simurg-pgxcn` (IN-050), region BED v1 built from UCSC ncbiRefSeq (`scripts/pgx/build_cn_regions.py`, +packaged as `simurg/pgx/data/pgx_cn_regions.GRCh38.v1.bed`), per-sample `*.pgx_cn.vcf.gz` + `.tbi`, +`pgx_cn_batch_summary.tsv`, `pgx_cn_sample_summary.tsv`. Unit tests: `tests/unit/test_pgx_cn.py`. + +## Data + +* **Arrays:** the GSA v4 PGx demo (`idats_Demo_36`, 36 arrays), genotyped by SimurgArray. Its CNV probes cover every + region (51–357 probes; GSTT1: 160 probes on `chr22_KI270879v1_alt` in the GRCh38 manifest). +* **Identities:** the demo arrays are 30 CEPH / 1000 Genomes individuals (6 replicate arrays). Identified by + genotype concordance with the 3202-sample 1000G 30x panel on 4647 chr22 SNPs: best match 0.992–0.996, next best + ≤ 0.91 (relatives) — `scripts/pgx/demo36_identity.tsv`. +* **Truth:** NYGC 30x SV callset (1KGP_3202 GATK-SV/svtools freeze V3), RD_CN per individual + (`scripts/pgx/validate_cn_truth.py`). +* The GSE165845 GSA-MD v1 arrays (360 arrays, CN/LOH validation) have 0–9 probes per PGx region: not usable here. + +## Result (30 individuals, model trained on the same 36 arrays with `--min-samples 30`) + +| Region | Concordance with NYGC RD_CN | Mismatches (truth → call) | +|---|---|---| +| GSTM1 (whole gene) | **30/30** | – | +| UGT2B17 (whole gene) | **30/30** | – | +| SULT1A1 | 29/30 | 2→1 | +| CYP2A6.exon.1 | 29/30 | 3→2 | +| CYP2A6.intron.7 | 28/30 | 3→2, 1→2 | +| CYP2D6.p5 | 28/30 | 4→3, 2→1 | +| CYP2D6.intron.2 | 25/30 | 3→2 ×3, 4→3, 2→1 | +| CYP2D6.exon.9 | 25/30 | 3→2 ×3, 4→2, 2→1 | +| **all** | **224/240 (93.3 %)** | | + +GSTT1 (no SV truth on the alt contig): CN 0/1/2 = 5/16/9, in Hardy–Weinberg proportions (null allele 0.43; CEU +literature ≈ 0.4–0.5). Details: `V2_pgx_cn_demo36.json`. + +* CYP2D6 truth is one gene-level record; the exon-9 / intron-2 "misses" at truth CN 3 are expected for tandem + structures such as `*68+*4` or `*36+*10`, where not every sub-region is duplicated. Star-allele truth (GeT-RM) + is the right yardstick and follows with the star-allele caller. +* PCN-017 target (≥ 95 % per region) is met for GSTM1, UGT2B17, SULT1A1, CYP2A6.exon.1; not (yet) shown for the + CYP2D6 sub-regions. + +## Caveats (K-028) + +* Trained and evaluated on the same 30 individuals (the spec requires ≥ 96 training samples; the priors were + chosen while looking at this truth). This is an optimistic, in-sample figure until an independent cohort is + available. +* Deviations from PCN-009: reference from the training cohort (not the calling batch); no GC-correction step; CN + model = compressed-intensity template + HWE-constrained weights (`simurg/pgx/cn.py`). diff --git a/docs/validation/V2_pgx_cn_demo36.json b/docs/validation/V2_pgx_cn_demo36.json new file mode 100644 index 0000000..168ff2b --- /dev/null +++ b/docs/validation/V2_pgx_cn_demo36.json @@ -0,0 +1,89 @@ +{ + "individuals": 30, + "regions": { + "GSTM1": { + "record": "HGSV_9691", + "concordant": 30, + "n": 30, + "rate": 1.0, + "truth_to_call_mismatches": {} + }, + "UGT2B17": { + "record": "HGSV_63272", + "concordant": 30, + "n": 30, + "rate": 1.0, + "truth_to_call_mismatches": {} + }, + "SULT1A1": { + "record": "HGSV_207805", + "concordant": 29, + "n": 30, + "rate": 0.9667, + "truth_to_call_mismatches": { + "2->1": 1 + } + }, + "CYP2A6.intron.7": { + "record": "HGSV_232737", + "concordant": 28, + "n": 30, + "rate": 0.9333, + "truth_to_call_mismatches": { + "3->2": 1, + "1->2": 1 + } + }, + "CYP2A6.exon.1": { + "record": "HGSV_232737", + "concordant": 29, + "n": 30, + "rate": 0.9667, + "truth_to_call_mismatches": { + "3->2": 1 + } + }, + "CYP2D6.exon.9": { + "record": "HGSV_247317", + "concordant": 25, + "n": 30, + "rate": 0.8333, + "truth_to_call_mismatches": { + "4->2": 1, + "3->2": 3, + "2->1": 1 + } + }, + "CYP2D6.intron.2": { + "record": "HGSV_247317", + "concordant": 25, + "n": 30, + "rate": 0.8333, + "truth_to_call_mismatches": { + "4->3": 1, + "3->2": 3, + "2->1": 1 + } + }, + "CYP2D6.p5": { + "record": "HGSV_247317", + "concordant": 28, + "n": 30, + "rate": 0.9333, + "truth_to_call_mismatches": { + "4->3": 1, + "2->1": 1 + } + } + }, + "overall": { + "concordant": 224, + "n": 240, + "rate": 0.9333 + }, + "gstt1_cn_distribution": { + "0": 5, + "1": 16, + "2": 9 + } +} \ No newline at end of file diff --git a/scripts/pgx/build_cn_regions.py b/scripts/pgx/build_cn_regions.py new file mode 100644 index 0000000..64f3cf0 --- /dev/null +++ b/scripts/pgx/build_cn_regions.py @@ -0,0 +1,85 @@ +#!/usr/bin/env python3 +"""Build the PGx CN target-region BED (spec PCN-001) from UCSC ncbiRefSeq (GRCh38). + + python scripts/pgx/build_cn_regions.py ncbiRefSeq.hg38.txt.gz OUT.bed [--version N] + +Regions (0-based half-open BED: chrom start0 end region gene), from the MANE/first-listed transcript: + CYP2D6.p5 2 kb 5' flank of CYP2D6 (upstream of exon 1; the gene is on the minus strand) + CYP2D6.intron.2 CYP2D6 intron 2 + CYP2D6.exon.9 CYP2D6 exon 9 + CYP2A6.exon.1 CYP2A6 exon 1 (± 500 bp) + CYP2A6.intron.7 CYP2A6 intron 7 + GSTM1, GSTT1 (on chr22_KI270879v1_alt, IN-062), UGT2B17, SULT1A1, CYP2E1: whole gene. +Exons and introns are numbered in transcript order (strand-aware). [DB: NCBI RefSeq via UCSC ncbiRefSeq] +""" + +from __future__ import annotations + +import argparse +import gzip +from pathlib import Path + +TRANSCRIPT = {"CYP2D6": "NM_000106", "CYP2A6": "NM_000762", "GSTM1": "NM_000561", "GSTT1": "NM_000853", + "UGT2B17": "NM_001077", "SULT1A1": "NM_001055", "CYP2E1": "NM_000773"} # fmt: skip +FLANK = 2000 +PAD = 500 + + +def main() -> int: + ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("refseq", type=Path) + ap.add_argument("out", type=Path) + ap.add_argument("--version", type=int, default=1) + a = ap.parse_args() + tx: dict[str, tuple[str, str, list[int], list[int]]] = {} + with gzip.open(a.refseq, "rt") as fh: + for line in fh: + f = line.rstrip("\n").split("\t") + name, chrom, strand, gene = f[1], f[2], f[3], f[12] + if gene not in TRANSCRIPT or not name.startswith(TRANSCRIPT[gene] + "."): + continue + if "_" in chrom and not (gene == "GSTT1" and chrom == "chr22_KI270879v1_alt"): + continue # primary assembly, except GSTT1 which is only on the alt contig + starts = [int(x) for x in f[9].rstrip(",").split(",")] + ends = [int(x) for x in f[10].rstrip(",").split(",")] + tx.setdefault(gene, (chrom, strand, starts, ends)) + missing = set(TRANSCRIPT) - set(tx) + if missing: + raise SystemExit(f"transcripts not found: {sorted(missing)}") + rows: list[tuple[str, int, int, str, str]] = [] + + def exon(gene: str, k: int) -> tuple[int, int]: + chrom, 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] + a_, b_ = exon(gene, k), exon(gene, k + 1) + return (a_[1], b_[0]) if strand == "+" else (b_[1], a_[0]) + + for gene in ("CYP2D6", "CYP2A6"): + chrom, strand, s, e = tx[gene] + if gene == "CYP2D6": + p5 = (e[-1], e[-1] + FLANK) if strand == "-" else (s[0] - FLANK, s[0]) + rows += [(chrom, *p5, "CYP2D6.p5", gene), (chrom, *intron(gene, 2), "CYP2D6.intron.2", gene), + (chrom, *exon(gene, 9), "CYP2D6.exon.9", gene)] # fmt: skip + else: + x1 = exon(gene, 1) + rows += [(chrom, x1[0] - PAD, x1[1] + PAD, "CYP2A6.exon.1", gene), + (chrom, *intron(gene, 7), "CYP2A6.intron.7", gene)] # fmt: skip + for gene in ("GSTM1", "GSTT1", "UGT2B17", "SULT1A1", "CYP2E1"): + chrom, strand, s, e = tx[gene] + rows.append((chrom, s[0], e[-1], gene, gene)) + order = {f"chr{i}": i for i in range(1, 23)} + rows.sort(key=lambda r: (order.get(r[0].split("_")[0], 99), r[0], r[1])) + with a.out.open("w") as fh: + fh.write(f"# simurg PGx CN regions v{a.version}, GRCh38, from UCSC ncbiRefSeq (PCN-001)\n") + for r in rows: + fh.write("\t".join(map(str, r)) + "\n") + print(a.out.read_text()) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/scripts/pgx/demo36_identity.tsv b/scripts/pgx/demo36_identity.tsv new file mode 100644 index 0000000..1dadfc9 --- /dev/null +++ b/scripts/pgx/demo36_identity.tsv @@ -0,0 +1,39 @@ +# GSA v4 PGx demo (idats_Demo_36) sample identities: best match among the 3202 1000 Genomes 30x individuals +# by genotype concordance on 4647 chr22 SNPs with 1000G AF 0.05-0.95 (next best <= 0.91: relatives) +sample_id individual concordance second_best second_concordance +207895980019_R03C01 NA12875 0.9920 NA12865 0.7038 +207895980019_R03C02 NA12156 0.9933 NA10831 0.7021 +207895980019_R04C01 NA12156 0.9931 NA10831 0.7027 +207895980019_R05C01 NA12044 0.9931 NA10857 0.6969 +207895980019_R06C01 NA11992 0.9931 NA10860 0.6848 +207895980019_R07C01 NA06993 0.9940 NA06991 0.7209 +207895980019_R07C03 NA07019 0.9953 NA07056 0.7051 +207895980019_R08C01 NA12239 0.9938 NA10847 0.7022 +207895980019_R09C01 NA12716 0.9944 NA12707 0.6896 +207895980019_R09C03 NA06991 0.9931 NA06993 0.7215 +207895980019_R10C01 NA12878 0.9933 NA12892 0.7022 +207895980019_R11C01 NA12760 0.9948 NA12752 0.7171 +207895980019_R11C02 NA06985 0.9948 NA06991 0.6996 +207895980019_R12C01 NA11993 0.9961 NA10860 0.7168 +207895980019_R12C02 NA12003 0.9935 NA10838 0.6998 +207895980019_R13C03 NA06994 0.9948 NA07029 0.7114 +207895980019_R14C03 NA11993 0.9959 NA10860 0.7163 +207895980019_R15C02 NA12003 0.9931 NA10838 0.6993 +207895980035_R03C03 NA07029 0.9952 NA07000 0.7237 +207895980035_R04C03 NA12892 0.9957 NA12878 0.7005 +207895980035_R05C03 NA12248 0.9950 NA10835 0.7143 +207895980035_R06C03 NA10846 0.9920 NA12145 0.7056 +207895980035_R07C01 NA07000 0.9957 NA07029 0.7231 +207895980035_R07C03 NA12155 0.9935 NA10831 0.6903 +207895980035_R08C01 NA11832 0.9942 NA10855 0.7252 +207895980035_R08C03 NA06993 0.9940 NA06991 0.7205 +207895980035_R09C03 NA10838 0.9948 NA12004 0.7007 +207895980035_R10C01 NA12813 0.9937 NA07045 0.9103 +207895980035_R10C03 NA12249 0.9938 NA10835 0.6858 +207895980035_R11C03 NA12812 0.9942 NA12801 0.7061 +207895980035_R12C03 NA11881 0.9935 NA10859 0.7138 +207895980035_R13C01 NA11995 0.9953 NA10861 0.6857 +207895980035_R13C03 NA12248 0.9950 NA10835 0.7145 +207895980035_R14C01 NA12891 0.9942 NA12878 0.6998 +207895980035_R15C03 NA12043 0.9946 NA10857 0.7162 +207895980042_R15C02 NA11832 0.9942 NA10855 0.7253 diff --git a/scripts/pgx/validate_cn_truth.py b/scripts/pgx/validate_cn_truth.py new file mode 100644 index 0000000..474049b --- /dev/null +++ b/scripts/pgx/validate_cn_truth.py @@ -0,0 +1,72 @@ +#!/usr/bin/env python3 +"""Validate `simurg pgx copy-number call` against the NYGC 30x SV callset (1KGP_3202 freeze V3) — PCN-017. + + python scripts/pgx/validate_cn_truth.py --calls CALL_DIR --identity scripts/pgx/demo36_identity.tsv \ + --sv 1KGP_3202.sv_freeze_V3.vcf.gz --bcftools BCFTOOLS --out summary.json + +Truth per region = RD_CN of one SV record per region (below, chosen as the record spanning the whole region; for +GSTM1 the read-depth DUP/CNV record, whose RD_CN gives 0/1/2 in Hardy–Weinberg proportions — the DEL record there +reports CN 1 for 27/30 individuals). CYP2D6 sub-regions share one gene-level record, so exon-9 conversions (*36) and +2D7-2D6 hybrids (*13, *68) legitimately differ from it. GSTT1 lies on an alt contig with no SV calls: only its CN +distribution is reported. Replicate arrays of one individual are counted once (the first array). +""" + +from __future__ import annotations + +import argparse +import collections +import csv +import json +import subprocess +from pathlib import Path + +TRUTH_RECORD = {"GSTM1": "HGSV_9691", "UGT2B17": "HGSV_63272", "SULT1A1": "HGSV_207805", + "CYP2A6.intron.7": "HGSV_232737", "CYP2A6.exon.1": "HGSV_232737", "CYP2D6.exon.9": "HGSV_247317", + "CYP2D6.intron.2": "HGSV_247317", "CYP2D6.p5": "HGSV_247317"} # fmt: skip + + +def main() -> int: + ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("--calls", required=True, type=Path) + ap.add_argument("--identity", required=True, type=Path) + ap.add_argument("--sv", required=True) + ap.add_argument("--bcftools", default="bcftools") + 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"] + calls: dict[str, dict[str, str]] = {} + with (a.calls / "pgx_cn_sample_summary.tsv").open() as fh: + for r in csv.DictReader(fh, delimiter="\t"): + ind = ident.get(r["sample_id"]) + if ind and ind not in calls: + calls[ind] = {k[3:]: v for k, v in r.items() if k.startswith("cn_")} + ids = sorted(set(TRUTH_RECORD.values())) + q = subprocess.run([a.bcftools, "query", "-H", "-s", ",".join(sorted(calls)), "--force-samples", "-i", + "ID=@/dev/stdin", "-f", "%ID[\t%RD_CN]\n", a.sv], input="\n".join(ids), capture_output=True, + text=True, check=True).stdout.splitlines() # fmt: skip + samples = [h.split("]")[1].split(":")[0] for h in q[0].split("\t")[1:]] + truth = {ln.split("\t")[0]: dict(zip(samples, ln.split("\t")[1:], strict=True)) for ln in q[1:]} + res: dict[str, dict[str, object]] = {} + tot = [0, 0] + for region, rid in TRUTH_RECORD.items(): + pairs = [(truth[rid][ind], c[region]) for ind, c in calls.items() if ind in truth.get(rid, {})] + pairs = [(t, c) for t, c in pairs if t not in (".", "")] + ok = sum(t == c for t, c in pairs) + tot[0] += ok + 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 + 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 + a.out.write_text(json.dumps(summary, indent=1)) + print(json.dumps(summary, indent=1)) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/src/simurg/cli/main.py b/src/simurg/cli/main.py index 38f9d6b..2051591 100644 --- a/src/simurg/cli/main.py +++ b/src/simurg/cli/main.py @@ -803,6 +803,78 @@ def cn_annotate_cmd( run_guarded(rep, lambda: run_cn_annotate(opts, rep)) +# -------------------------------------------------------------------------------------------------------- +pgx_app = typer.Typer(help="Pharmacogenomics: targeted copy number, star alleles, phenotypes (v2).") +app.add_typer(pgx_app, name="pgx") +pgx_cn_app = typer.Typer(help="Targeted PGx copy number (PCN-001).") +pgx_app.add_typer(pgx_cn_app, name="copy-number") + + +@pgx_cn_app.command("train") +def pgx_cn_train_cmd( + gtcs: Annotated[Path, typer.Option("--gtcs", help="Folder searched recursively for training GTCs (required).")], + bpm: Annotated[Path, typer.Option("--bpm", help="BPM manifest (required).")], + fasta: Annotated[Path, typer.Option("--fasta", help="GRCh38 FASTA (required).")], + out: Annotated[Path, typer.Option("--out", help="Output folder.")] = Path("."), + mask: Annotated[Path | None, typer.Option("--mask", help="Probe mask TSV (IN-052).")] = None, + regions: Annotated[Path | None, typer.Option("--regions", help="Region BED (default: built-in v1).")] = None, + summary: Annotated[Path | None, typer.Option("--summary", help="sample_summary.tsv (assay control).")] = None, + platform: Annotated[str, typer.Option("--platform", help="auto|LCG|EX")] = "auto", + min_samples: Annotated[int, typer.Option("--min-samples", min=1)] = 96, + min_passing: Annotated[int, typer.Option("--min-passing", min=1)] = 90, + max_samples: Annotated[int, typer.Option("--max-samples", min=1)] = 3000, + train_lrr_sd_max: Annotated[float, typer.Option("--train-lrr-sd-max")] = 0.20, + train_call_rate_min: Annotated[float, typer.Option("--train-call-rate-min")] = 0.99, + train_assay_control_min: Annotated[float, typer.Option("--train-assay-control-min")] = 1.0, + force: ForceOpt = False, + log_json: LogJson = False, + debug: Debug = False, + quiet: Quiet = False, +) -> None: + """Train a PGx CN model (CLI-032, PCN-005).""" + from simurg.pipeline.pgx_cn import PgxTrainOptions, run_train + + if platform not in ("auto", "LCG", "EX"): + raise typer.BadParameter("--platform auto|LCG|EX") + rep = _reporter("pgx copy-number train", out, log_json, debug, quiet) + o = PgxTrainOptions(gtcs=gtcs, bpm=bpm, fasta=fasta, out=out, mask=mask, regions=regions, summary=summary, + platform=platform, min_samples=min_samples, min_passing=min_passing, max_samples=max_samples, + lrr_sd_max=train_lrr_sd_max, call_rate_min=train_call_rate_min, + assay_control_min=train_assay_control_min, force=force) # fmt: skip + run_guarded(rep, lambda: run_train(o, rep)) + + +@pgx_cn_app.command("call") +def pgx_cn_call_cmd( + gtcs: Annotated[Path, typer.Option("--gtcs", help="Folder searched recursively for GTC files (required).")], + model: Annotated[Path, typer.Option("--model", help="PGx CN model (.simurg-pgxcn, required).")], + fasta: Annotated[Path, typer.Option("--fasta", help="GRCh38 FASTA with alt contigs (required).")], + out: Annotated[Path, typer.Option("--out", help="Output folder (required).")], + min_batch: Annotated[int, typer.Option("--min-batch", min=1)] = 24, + min_passing: Annotated[int, typer.Option("--min-passing", min=1)] = 22, + recommended_passing: Annotated[int, typer.Option("--recommended-passing", min=1)] = 96, + lrr_sd_max: Annotated[float, typer.Option("--lrr-sd-max")] = 0.20, + cluster_mismatch: Annotated[str, typer.Option("--cluster-mismatch", help="error|warn")] = "error", + no_bgzip: Annotated[bool, typer.Option("--no-bgzip")] = False, + force: ForceOpt = False, + log_json: LogJson = False, + debug: Debug = False, + quiet: Quiet = False, +) -> None: + """Call targeted PGx copy number per region (CLI-031, PCN-009).""" + from simurg.pipeline.pgx_cn import PgxCallOptions, run_call + + if cluster_mismatch not in ("error", "warn"): + raise typer.BadParameter("--cluster-mismatch error|warn") + if min_passing > min_batch: + raise typer.BadParameter("--min-passing must be <= --min-batch") + rep = _reporter("pgx copy-number call", out, log_json, debug, quiet) + o = PgxCallOptions(gtcs=gtcs, model=model, fasta=fasta, out=out, min_batch=min_batch, min_passing=min_passing, + recommended_passing=recommended_passing, lrr_sd_max=lrr_sd_max, + cluster_mismatch=cluster_mismatch, bgzip=not no_bgzip, force=force) # fmt: skip + run_guarded(rep, lambda: run_call(o, rep)) + + # -------------------------------------------------------------------------------------------------------- methyl_app = typer.Typer(help="Methylation arrays: control QC, detection, beta values (MET-001).") app.add_typer(methyl_app, name="methyl") diff --git a/src/simurg/core/errors.py b/src/simurg/core/errors.py index 404c278..9869a54 100644 --- a/src/simurg/core/errors.py +++ b/src/simurg/core/errors.py @@ -151,6 +151,26 @@ CATALOGUE: dict[str, MessageDef] = dict( _m("E_CN_SIDECAR_MISSING", "error", "cn_call_manifest.json not found in the cn call output", "CYTO-032"), _m("E_CN_VCF_UNINDEXED", "error", "CNV VCFs must be bgzipped and indexed for annotate", "CYTO-032"), _m("E_CN_FAILED", "error", "CN calling failed for the sample", "CYTO-040"), + # --- targeted PGx copy number (PCN-001 .. PCN-017) --------------------------------------------------- + _m("E_PGX_BATCH_TOO_SMALL", "error", "Too few (passing) samples in the batch for PGx CN calling", "PCN-004"), + _m("W_PGX_BATCH_BELOW_RECOMMENDED", "warning", "Fewer passing samples than recommended", "PCN-004"), + _m("E_PGX_TRAIN_TOO_SMALL", "error", "Too few (passing) samples to train a PGx CN model", "PCN-005"), + _m( + "E_MODEL_CLUSTER_MISMATCH", + "error", + "Input GTCs were called with another cluster file than the model", + "IN-050", + ), + _m( + "W_MODEL_CLUSTER_MISMATCH", + "warning", + "Input GTCs were called with another cluster file than the model", + "IN-050", + ), + _m("W_PGX_REGION_FEW_PROBES", "warning", "PGx CN region has fewer than 10 unmasked probes", "PCN-008"), + _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"), # --- methylation control QC (MET-001 .. MET-018) ------------------------------------------------- _m("E_METH_NO_ASSAY", "error", "Probe section not found", "MET-002"), _m("E_METH_NO_CONTROLS", "error", "Missing control probes in manifest", "MET-002"), diff --git a/src/simurg/pgx/cn.py b/src/simurg/pgx/cn.py new file mode 100644 index 0000000..a9fc740 --- /dev/null +++ b/src/simurg/pgx/cn.py @@ -0,0 +1,204 @@ +"""Targeted PGx copy number — spec PCN-001, PCN-009. Pure functions, no I/O. + +Region value x = median over the region's unmasked probes of (sample-centred LRR − the probe's training reference), +minus the sample's scale (median of the same quantity over non-target autosomal probes). + +CN mixture [OURS; array intensity compression as in cnvPartition/PennCNV, LIT: Wang et al. 2007 Genome Res 17:1665]: +component means µ_c = δ_r + γ_r · log2((c + 1) / 3) for CN c = 0 … 5 — a probe signal proportional to (c + β) with +a background β = 1 copy-equivalent, compressed by γ_r. δ_r (the region's CN-2 location), γ_r, σ_r and the weights +π_r are per region; γ_r has a Gaussian prior N(0.25, 0.08²). CN 0 may sit lower than the template (its own shift +ε₀ ≤ 0 and SD): GSTT1 probes drop by ~2 at CN 0, GSTM1 probes by ~0.3. The weights are constrained to Hardy–Weinberg +proportions of haplotype copy numbers (0 … 3 copies per haplotype, Dirichlet prior favouring one copy). Fitted by +EM from a grid of starts (each observed cluster labelled CN 0 … 4, CN 0 on the template or on the lowest values). +HWE and the template's spacing ratios identify the labels in multi-cluster regions (GSTM1 0/1/2 rather than 1/2/3, +which would imply unobserved CN-0 samples). On GSA v4 PGx CNV probes γ_r ≈ 0.19–0.56 (K-028). +""" + +from __future__ import annotations + +import math +from dataclasses import dataclass + +import numpy as np +import numpy.typing as npt + +F64 = npt.NDArray[np.float64] +CN_MAX = 5 # "5 or more" +CNS = np.arange(CN_MAX + 1) +TEMPLATE = np.log2((CNS + 1) / 3.0) +HAP_ALPHA = np.array([1.2, 3.0, 1.2, 1.05]) # Dirichlet prior on haplotype copy frequencies f_0..f_3 (f_1 favoured) +SIGMA_FLOOR = 0.012 +GAMMA_PRIOR = (0.25, 0.08) # mean, SD +GAMMA_RANGE = (0.15, 0.80) +PI_FLOOR = 0.002 +CN0_PRIOR_N = 2.0 + + +@dataclass(frozen=True) +class Region: + chrom: str + start0: int + end: int + name: str + gene: str + + +def parse_regions(text: str) -> list[Region]: + """PCN-001 BED: ``chrom start0 end region gene``; ``#`` comments.""" + out = [] + for ln in text.splitlines(): + if not ln.strip() or ln.startswith("#"): + continue + c, s, e, r, g = ln.split("\t")[:5] + out.append(Region(c, int(s), int(e), r, g)) + return out + + +@dataclass +class RegionFit: + delta: float + gamma: float + sigma: float + weights: list[float] + n: int + cn0_shift: float = 0.0 # ≤ 0: extra drop of the CN-0 component below the template (probe-specific background) + cn0_sigma: float = 0.0 # SD of the CN-0 component (≥ sigma; 0 = sigma) + + def means(self) -> F64: + m = self.delta + self.gamma * TEMPLATE + m[0] += self.cn0_shift + return m + + def sds(self) -> F64: + s = np.full(len(CNS), self.sigma) + s[0] = max(self.sigma, self.cn0_sigma) + return s + + +def hwe_weights(f: F64) -> F64: + """CN weights under Hardy–Weinberg from haplotype copy-number frequencies f_k (k = 0 … 3 copies): + π_c = Σ_{i+j=c} f_i f_j, CN ≥ 5 pooled into 5 [LIT: HWE; OURS: use as the CN mixture prior].""" + pi = np.zeros(len(CNS)) + for i in range(len(f)): + for j in range(len(f)): + pi[min(i + j, CN_MAX)] += f[i] * f[j] + return pi + + +def _hap_update(r: F64, f: F64) -> F64: + """M-step for f: split each sample's CN responsibilities over the haplotype pairs (i, j) with i + j = c.""" + counts = HAP_ALPHA - 1.0 + counts = counts.astype(np.float64).copy() + pi = hwe_weights(f) + for i in range(len(f)): + for j in range(len(f)): + c = min(i + j, CN_MAX) + share = f[i] * f[j] / max(pi[c], 1e-300) + w = float(r[:, c].sum()) * share + counts[i] += w + counts[j] += w + counts = np.maximum(counts, 1e-9) + out: F64 = counts / counts.sum() + return out + + +def _norm_pdf(x: F64, mu: F64, sd: float | F64) -> F64: + sd_ = np.broadcast_to(np.asarray(sd, dtype=np.float64), mu.shape)[None, :] + return np.exp(-0.5 * ((x[:, None] - mu[None, :]) / sd_) ** 2) / (sd_ * math.sqrt(2 * math.pi)) + + +def _em_region(x: F64, gamma: float, delta: float, iters: int, eps0: float = 0.0) -> tuple[RegionFit, float]: + """EM for one region; returns the fit and its penalised log-likelihood (HWE weights from haplotype copy + frequencies with a Dirichlet prior, Gaussian prior on γ). + CN 0 has its own shift ε₀ ≤ 0 below the template and its own SD, both shrunk towards the template with a prior + weight of ``CN0_PRIOR_N`` samples (on some probe sets CN 0 is far lower than the compressed template, GSTT1).""" + n = len(x) + f = np.array([0.1, 0.8, 0.07, 0.03]) + pi = hwe_weights(f) + sd = max(SIGMA_FLOOR, float(np.std(x)) if n > 1 else 0.05) + sd0 = sd + if eps0 < 0: # a start with a far CN-0 cluster: the other components get the spread of the rest + rest = x[x > delta + gamma * TEMPLATE[0] + eps0 / 2] + sd = max(SIGMA_FLOOR, float(np.std(rest)) if len(rest) > 1 else 0.05) + sd0 = sd + g0, gs = GAMMA_PRIOR + ll = -math.inf + for _ in range(iters): + mu = delta + gamma * TEMPLATE + mu[0] += eps0 + sds = np.full(len(CNS), sd) + sds[0] = sd0 + dens = pi[None, :] * _norm_pdf(x, mu, sds) + 1e-300 + tot = dens.sum(axis=1, keepdims=True) + ll = float(np.log(tot).sum() + ((HAP_ALPHA - 1) * np.log(np.maximum(f, 1e-12))).sum() + - 0.5 * ((gamma - g0) / gs) ** 2) # fmt: skip + r = dens / tot + f = _hap_update(r, f) + pi = hwe_weights(f) + # δ and γ: weighted least squares of x on (1, t_c) over CN ≥ 1, γ shrunk towards its prior + rr = r.copy() + rr[:, 0] = 0.0 + w = rr.sum(axis=0) + sw, st, stt = float(w.sum()), float((w * TEMPLATE).sum()), float((w * TEMPLATE**2).sum()) + sx, stx = float((rr * x[:, None]).sum()), float((rr * x[:, None] * TEMPLATE[None, :]).sum()) + lam = sd**2 / gs**2 + det = sw * (stt + lam) - st * st + if det > 0: + delta = (sx * (stt + lam) - st * (stx + lam * g0)) / det + gamma = float(np.clip((sw * (stx + lam * g0) - st * sx) / det, *GAMMA_RANGE)) + mu = delta + gamma * TEMPLATE + n1 = float(w.sum()) + if n1 > 0: + sd = max(SIGMA_FLOOR, math.sqrt(float((rr * (x[:, None] - mu[None, :]) ** 2).sum() / n1))) + r0 = r[:, 0] + n0 = float(r0.sum()) + eps0 = min(0.0, float((r0 * (x - mu[0])).sum()) / (n0 + CN0_PRIOR_N)) + v0 = float((r0 * (x - mu[0] - eps0) ** 2).sum()) + sd0 = max(sd, math.sqrt((v0 + CN0_PRIOR_N * sd**2) / (n0 + CN0_PRIOR_N))) + return RegionFit(float(delta), float(gamma), sd, pi.tolist(), n, float(eps0), float(sd0)), ll + + +def fit_region(x: F64, iters: int = 200) -> RegionFit: + """PCN-009 step 4 for one region: the best of EM runs started with each observed cluster labelled CN 0 … 4, each + with CN 0 on the template or on the lowest values.""" + g0 = GAMMA_PRIOR[0] + cands = sorted({round(float(q) - g0 * float(t), 4) for q in np.quantile(x, [0.1, 0.3, 0.5, 0.7, 0.9]) + for t in TEMPLATE[:5]}) # fmt: skip + lo = float(np.quantile(x, 0.02)) + best: tuple[float, float, float] | None = None + for d0 in cands: + for e0 in sorted({0.0, round(min(0.0, lo - (d0 + g0 * float(TEMPLATE[0]))), 4)}): + _, ll = _em_region(x, g0, d0, 80, e0) + if best is None or ll > best[0] + 1e-9: + best = (ll, d0, e0) + assert best is not None + return _em_region(x, g0, best[1], iters, best[2])[0] + + +def fit_mixture(values: dict[str, F64], iters: int = 200) -> dict[str, RegionFit]: + """Per-region fits for every region with ≥ 3 finite values.""" + out = {} + for name, x in values.items(): + x = x[np.isfinite(x)] + if len(x) >= 3: + out[name] = fit_region(x, iters) + return out + + +def posterior(x: float, fit: RegionFit) -> F64: + """P(CN = c | x) for c = 0 … 5, with the training weights (floored) as prior.""" + pi = np.maximum(np.asarray(fit.weights), PI_FLOOR) + pi = pi / pi.sum() + d = pi * _norm_pdf(np.array([x]), fit.means(), fit.sds())[0] + s = float(d.sum()) + out: F64 = d / s if s > 0 else pi + return out + + +def call_cn(x: float, fit: RegionFit) -> tuple[int, int, F64]: + """PCN-009 steps 5–6: (CN, QUAL, posterior).""" + p = posterior(x, fit) + cn = int(np.argmax(p)) + pmax = float(p[cn]) + qual = 60 if pmax >= 1.0 else int(min(60, max(0, round(-10 * math.log10(max(1e-300, 1 - pmax)))))) + return cn, qual, p diff --git a/src/simurg/pgx/data/pgx_cn_regions.GRCh38.v1.bed b/src/simurg/pgx/data/pgx_cn_regions.GRCh38.v1.bed new file mode 100644 index 0000000..9919272 --- /dev/null +++ b/src/simurg/pgx/data/pgx_cn_regions.GRCh38.v1.bed @@ -0,0 +1,11 @@ +# simurg PGx CN regions v1, GRCh38, from UCSC ncbiRefSeq (PCN-001) +chr1 109687816 109693745 GSTM1 GSTM1 +chr4 68537172 68576322 UGT2B17 UGT2B17 +chr10 133527362 133539123 CYP2E1 CYP2E1 +chr16 28605257 28610044 SULT1A1 SULT1A1 +chr19 40844772 40845293 CYP2A6.intron.7 CYP2A6 +chr19 40849746 40850947 CYP2A6.exon.1 CYP2A6 +chr22 42126498 42126752 CYP2D6.exon.9 CYP2D6 +chr22 42129185 42129737 CYP2D6.intron.2 CYP2D6 +chr22 42130810 42132810 CYP2D6.p5 CYP2D6 +chr22_KI270879v1_alt 270313 278446 GSTT1 GSTT1 diff --git a/src/simurg/pipeline/pgx_cn.py b/src/simurg/pipeline/pgx_cn.py new file mode 100644 index 0000000..0ca78e1 --- /dev/null +++ b/src/simurg/pipeline/pgx_cn.py @@ -0,0 +1,510 @@ +"""``simurg pgx copy-number train|call`` — targeted PGx copy number (spec PCN-001 … PCN-017, CLI-031/032, IN-050). + +I/O and orchestration only; the numerics live in ``simurg.pgx.cn``. + +Deviations (K-028): the per-probe reference is the training cohort's median, stored in the model (not the calling +batch's median, which is not CN 2 in regions with common deletions such as GSTM1); no GC correction step (the region +median over many probes plus the sample scale already removes most of it); GSTT1 needs alt-contig probe coordinates, +which the GSA v4 PGx GRCh38 manifest carries (160 probes on `chr22_KI270879v1_alt`); manifests without them give +`NoCall` for GSTT1. +""" + +from __future__ import annotations + +import hashlib +import io +import json +import math +import zipfile +from dataclasses import dataclass, field +from datetime import UTC, datetime +from importlib import resources +from pathlib import Path +from typing import Any + +import numpy as np +import numpy.typing as npt +import pyarrow as pa +import pyarrow.parquet as pq + +from simurg import __version__ +from simurg.core.errors import ExitCode, SimurgError +from simurg.core.log import Reporter +from simurg.formats.bpm import read_bpm +from simurg.pgx.cn import CN_MAX, Region, RegionFit, call_cn, fit_mixture, parse_regions + +MODEL_SCHEMA = 1 +REGIONS_FILE = "pgx_cn_regions.GRCh38.v1.bed" +N_SCALE_PROBES = 20_000 +MIN_REGION_PROBES = 10 # PCN-008 +MIN_FINITE_FRAC = 0.5 # PCN-002: fewer finite region values → NoCall [OURS] +F32 = npt.NDArray[np.float32] + + +def default_regions() -> tuple[str, str]: + """(file name, text) of the packaged region BED (PCN-001).""" + return REGIONS_FILE, resources.files("simurg.pgx.data").joinpath(REGIONS_FILE).read_text(encoding="utf-8") + + +def _core(chrom: str) -> str: + c = chrom[3:] if chrom.lower().startswith("chr") else chrom + return "X" if c == "XY" else c + + +@dataclass +class PgxCnModel: + header: dict[str, Any] + regions: list[Region] + regions_text: str + probe_idx: npt.NDArray[np.int64] # manifest indices of all probes used (region probes, then scale probes) + probe_region: list[str] # region name, "" for scale probes + probe_ref: F32 # training reference (median sample-centred LRR) + fits: dict[str, RegionFit] = field(default_factory=dict) + + def region_sel(self, name: str) -> npt.NDArray[np.int64]: + return np.flatnonzero(np.asarray(self.probe_region) == name) + + +def region_probes( + regions: list[Region], chroms: list[str], pos: npt.NDArray[np.int64], masked: set[int] +) -> dict[str, npt.NDArray[np.int64]]: + """Manifest indices inside each region (1-based position in (start0, end]), mask removed (PCN-008).""" + core = np.array([_core(c) for c in chroms]) + out = {} + for r in regions: + sel = np.flatnonzero((core == _core(r.chrom)) & (pos > r.start0) & (pos <= r.end)) + out[r.name] = np.array([i for i in sel.tolist() if i not in masked], dtype=np.int64) + return out + + +def scale_probes( + chroms: list[str], pos: npt.NDArray[np.int64], exclude: npt.NDArray[np.int64], n: int = N_SCALE_PROBES +) -> npt.NDArray[np.int64]: + """Seeded random autosomal non-target probes for the per-sample scale (PCN-009 step 2).""" + core = np.array([_core(c) for c in chroms]) + ok = np.isin(core, [str(i) for i in range(1, 23)]) & (pos > 0) + ok[exclude] = False + cand = np.flatnonzero(ok) + rng = np.random.default_rng(20260929) + return np.sort(rng.choice(cand, min(n, len(cand)), replace=False)).astype(np.int64) + + +def _centred(lrr: npt.NDArray[np.float32], auto_idx: npt.NDArray[np.int64]) -> npt.NDArray[np.float64]: + x = lrr.astype(np.float64) + x[~np.isfinite(x)] = np.nan + med = float(np.nanmedian(x[auto_idx])) if np.isfinite(x[auto_idx]).any() else 0.0 + return x - med + + +def region_values(d: npt.NDArray[np.float64], model: PgxCnModel) -> dict[str, tuple[float, int]]: + """Per region: (x, n finite probes). ``d`` = centred LRR at ``model.probe_idx`` minus the reference.""" + scale_sel = model.region_sel("") + sc = float(np.nanmedian(d[scale_sel])) if np.isfinite(d[scale_sel]).any() else 0.0 + out = {} + for r in model.regions: + sel = model.region_sel(r.name) + v = d[sel] + nf = int(np.isfinite(v).sum()) + out[r.name] = (float(np.nanmedian(v)) - sc if nf else math.nan, nf) + return out + + +# ------------------------------------------------------------------------------------------------ model file +def write_model(path: Path, m: PgxCnModel) -> None: + table = pa.table({"index": m.probe_idx, "region": m.probe_region, "ref": m.probe_ref, + "name": m.header.pop("_names")}) # fmt: skip + buf = io.BytesIO() + pq.write_table(table, buf, compression="zstd") + mixture = {name: {"delta": f.delta, "gamma": f.gamma, "sigma": f.sigma, "weights": f.weights, "n": f.n, + "cn0_shift": f.cn0_shift, "cn0_sigma": f.cn0_sigma} for name, f in m.fits.items()} # fmt: skip + path.parent.mkdir(parents=True, exist_ok=True) + tmp = path.with_suffix(path.suffix + ".part") + with zipfile.ZipFile(tmp, "w", compression=zipfile.ZIP_STORED) as z: + z.writestr("header.json", json.dumps(m.header, indent=1)) + z.writestr("probes.parquet", buf.getvalue()) + z.writestr("regions.bed", m.regions_text) + z.writestr("mixture.json", json.dumps(mixture, indent=1)) + tmp.replace(path) + + +def read_model(path: Path) -> PgxCnModel: + """IN-050.""" + with zipfile.ZipFile(path) as z: + header = json.loads(z.read("header.json")) + table = pq.read_table(io.BytesIO(z.read("probes.parquet"))) + regions_text = z.read("regions.bed").decode() + mixture = json.loads(z.read("mixture.json")) + if int(header.get("schema_version", 0)) < MODEL_SCHEMA: + raise SimurgError("E_MODEL_VERSION", f"{path}: PGx CN model schema {header.get('schema_version')} is too old", + hint="retrain it with simurg pgx copy-number train") # fmt: skip + fits = { + k: RegionFit( + v["delta"], + v["gamma"], + v["sigma"], + list(v["weights"]), + int(v["n"]), + v.get("cn0_shift", 0.0), + v.get("cn0_sigma", 0.0), + ) + for k, v in mixture.items() + } + return PgxCnModel(header, parse_regions(regions_text), regions_text, + table.column("index").to_numpy().astype(np.int64), table.column("region").to_pylist(), + table.column("ref").to_numpy().astype(np.float32), fits) # fmt: skip + + +# ------------------------------------------------------------------------------------------------ train +@dataclass +class PgxTrainOptions: + gtcs: Path + bpm: Path + fasta: Path + out: Path + mask: Path | None = None + regions: Path | None = None + summary: Path | None = None + platform: str = "auto" + min_samples: int = 96 + min_passing: int = 90 + max_samples: int = 3000 + lrr_sd_max: float = 0.20 + call_rate_min: float = 0.99 + assay_control_min: float = 1.0 + force: bool = False + argv: list[str] = field(default_factory=list) + + +def read_mask_tsv(path: Path | None) -> tuple[set[str], str]: + """IN-052: ``probe_namereason``; returns (names, sha256).""" + if path is None: + return set(), "" + data = path.read_bytes() + names = set() + for ln in data.decode("utf-8").splitlines(): + if ln.strip() and not ln.startswith("#"): + names.add(ln.split("\t")[0].strip()) + return names, hashlib.sha256(data).hexdigest() + + +def _summary(path: Path | None) -> dict[str, dict[str, str]]: + if path is None or not path.exists(): + return {} + lines = path.read_text(encoding="utf-8").splitlines() + head = lines[0].split("\t") + return {r[0]: dict(zip(head, r, strict=False)) for r in (ln.split("\t") for ln in lines[1:] if ln)} + + +def _assay_control(row: dict[str, str]) -> float: + for k, v in row.items(): + if k.endswith("_norm_r") and k.startswith("TGA_Ctrl"): + try: + return float(v) + except ValueError: + return math.nan + return math.nan + + +def _probe_set_sha(names: list[str], chroms: list[str], pos: npt.NDArray[np.int64]) -> str: + text = "\n".join(f"{n}\t{c}\t{p}" for n, c, p in zip(names, chroms, pos.tolist(), strict=True)) + return hashlib.sha256(text.encode()).hexdigest() + + +def run_train(o: PgxTrainOptions, rep: Reporter) -> int: + """``simurg pgx copy-number train`` (CLI-032, PCN-005).""" + from simurg.core.atomic import check_outputs_absent + from simurg.formats.discovery import discover_gtcs + from simurg.formats.fasta import Genome + from simurg.formats.gtc import read_gtc + from simurg.pipeline.qc import write_tsv + + man = read_bpm(o.bpm) + out_model = o.out / f"{o.bpm.stem}.simurg-pgxcn" + check_outputs_absent([out_model, o.out / "train_exclusions.tsv"], force=o.force) + genome = Genome(o.fasta) + if genome.build not in ("GRCh38", "unknown", ""): + raise SimurgError("E_BUILD_MISMATCH", f"PGx CN supports GRCh38 only (FASTA is {genome.build})") # PCN-015 + reg_name, reg_text = (o.regions.name, o.regions.read_text()) if o.regions else default_regions() + regions = parse_regions(reg_text) + mask_path = o.mask or ( + o.bpm.with_name(o.bpm.stem + ".mask.tsv") if (o.bpm.with_name(o.bpm.stem + ".mask.tsv")).exists() else None + ) + mask_names, mask_sha = read_mask_tsv(mask_path) + masked = {i for i, n in enumerate(man.names) if n in mask_names} + pos = man.map_infos.astype(np.int64) + rp = region_probes(regions, list(man.chroms), pos, masked) + for r in regions: + if len(rp[r.name]) < MIN_REGION_PROBES: + rep.warn("W_PGX_REGION_FEW_PROBES", f"region {r.name}: {len(rp[r.name])} unmasked probe(s) on " + f"{man.manifest_name}; it will be NoCall") # fmt: skip + reg_idx = np.concatenate([rp[r.name] for r in regions]) if regions else np.zeros(0, np.int64) + sc_idx = scale_probes(list(man.chroms), pos, reg_idx) + probe_idx = np.concatenate([reg_idx, sc_idx]) + probe_region = [r.name for r in regions for _ in range(len(rp[r.name]))] + [""] * len(sc_idx) + gtcs = discover_gtcs(o.gtcs) + keys = sorted(gtcs) + if len(keys) > o.max_samples: + rng = np.random.default_rng(42) + keys = sorted(rng.choice(keys, o.max_samples, replace=False).tolist()) + rep.info(f"training on a seeded random subset of {o.max_samples} of {len(gtcs)} samples (PCN-005)") + summ = _summary(o.summary or (o.gtcs / "sample_summary.tsv")) + rows = np.full((len(keys), len(probe_idx)), np.nan, dtype=np.float32) + exclusions: list[dict[str, str]] = [] + passing = np.zeros(len(keys), dtype=np.bool_) + cluster_files = set() + for i, k in enumerate(keys): + g = read_gtc(gtcs[k]) + if len(g.logr_ratios) != man.num_loci: + raise SimurgError("E_MODEL_MANIFEST_MISMATCH", f"{gtcs[k]} has {len(g.logr_ratios)} loci, the manifest " + f"{man.num_loci}") # fmt: skip + cluster_files.add(g.cluster_file) + rows[i] = _centred(g.logr_ratios, sc_idx)[probe_idx] + reasons = [] + if not g.logr_dev <= o.lrr_sd_max: + reasons.append(f"logr_dev {g.logr_dev:.3f} > {o.lrr_sd_max}") + if not g.call_rate >= o.call_rate_min: + reasons.append(f"call_rate {g.call_rate:.4f} < {o.call_rate_min}") + ac = _assay_control(summ.get(k, {})) + if math.isfinite(ac) and ac < o.assay_control_min: + reasons.append(f"assay_control_norm_r {ac:.3f} < {o.assay_control_min}") + if reasons: + exclusions.append({"sample_id": k, "reason": "; ".join(reasons)}) + else: + passing[i] = True + n_pass = int(passing.sum()) + if len(keys) < o.min_samples or n_pass < o.min_passing: + raise SimurgError("E_PGX_TRAIN_TOO_SMALL", f"{n_pass} of {len(keys)} sample(s) passed QC; training requires " + f"at least {o.min_samples} samples of which {o.min_passing} pass", + hint="lower --min-samples/--min-passing only for research use") # fmt: skip + ref = np.nanmedian(rows[passing], axis=0) + d = rows - ref + sc_sel = np.flatnonzero(np.asarray(probe_region) == "") + sc = np.nanmedian(d[:, sc_sel], axis=1) + mad = float(np.nanmedian(np.abs(sc[passing] - np.nanmedian(sc[passing])))) * 1.4826 or 1e-9 + z = np.abs(sc - np.nanmedian(sc[passing])) / mad + use = passing & (z <= 5) # PCN-009: outlier samples out of the reference + for i in np.flatnonzero(passing & ~use).tolist(): + exclusions.append({"sample_id": keys[i], "reason": f"genome-wide NR outlier (|z| = {z[i]:.1f} > 5)"}) + if (~use & passing).any(): + ref = np.nanmedian(rows[use], axis=0) + d = rows - ref + model = PgxCnModel({}, regions, reg_text, probe_idx, probe_region, ref.astype(np.float32)) + values = {} + for r in regions: + if len(rp[r.name]) < MIN_REGION_PROBES: + continue + values[r.name] = np.array([region_values(d[i], model)[r.name][0] for i in np.flatnonzero(use)]) + model.fits = fit_mixture(values) + csv_fmt = "unknown" if o.platform == "auto" else o.platform + model.header = { + "schema": "simurg.pgx_cn_model", "schema_version": MODEL_SCHEMA, "software_version": __version__, + "manifest_name": man.manifest_name, "manifest_revision": getattr(man, "version", None), + "probe_set_sha256": _probe_set_sha(man.names, list(man.chroms), pos), + "n_probes": man.num_loci, + "cluster_file_name": sorted(cluster_files)[0] if len(cluster_files) == 1 else sorted(cluster_files), + "genome_build": "GRCh38", "fasta_contigs_sha256": genome.sha256, "platform": csv_fmt, + "mask_sha256": mask_sha, "regions_bed_version": reg_name, "n_training_samples": int(use.sum()), + "training_exclusions": len(exclusions), "created": datetime.now(UTC).strftime("%Y-%m-%dT%H:%M:%SZ"), + "regions_few_probes": [r.name for r in regions if len(rp[r.name]) < MIN_REGION_PROBES], + "model": "mu_c = delta_r + gamma_r * log2((c + 1) / 3), c = 0..5 (simurg.pgx.cn)", + "_names": [man.names[i] for i in probe_idx.tolist()], + } # fmt: skip + o.out.mkdir(parents=True, exist_ok=True) + write_model(out_model, model) + write_tsv(o.out / "train_exclusions.tsv", ["sample_id", "reason"], exclusions) + for name, f in model.fits.items(): + rep.info(f"{name}: CN2 at {f.delta:+.3f}, gamma {f.gamma:.3f}, sd {f.sigma:.3f}, " + f"weights {' '.join(f'{w:.2f}' for w in f.weights)}") # fmt: skip + rep.info(f"PGx CN model from {int(use.sum())} sample(s) written to {out_model}") + return ExitCode.OK + + +# ------------------------------------------------------------------------------------------------ call +@dataclass +class PgxCallOptions: + gtcs: Path + model: Path + fasta: Path + out: Path + min_batch: int = 24 + min_passing: int = 22 + recommended_passing: int = 96 + lrr_sd_max: float = 0.20 + cluster_mismatch: str = "error" + bgzip: bool = True + force: bool = False + argv: list[str] = field(default_factory=list) + + +def _vcf_header(sample: str, contigs: list[tuple[str, int]], meta: dict[str, str]) -> str: + h = ["##fileformat=VCFv4.3", f"##source=SimurgArray {__version__} pgx-cn", f"##reference=file://{meta['fasta']}", + "##simurg_genome_build=GRCh38", f"##simurg_manifest={meta['manifest']}", + f"##simurg_pgx_cn_model={meta['model']}", f"##simurg_batch={meta['batch']}", + f"##OverallPloidy={meta['ploidy']}", "##GCCorrect=false"] # fmt: skip + h += [f"##contig=" for c, n in contigs] + h += ['##ALT=', + '##FILTER=', + '##FILTER=', + '##FILTER=', + '##FILTER=', + '##INFO=', + '##INFO=', + '##INFO=', + '##INFO=', + '##INFO=', + '##INFO=', + '##FORMAT=', + '##FORMAT=', + "\t".join(["#CHROM", "POS", "ID", "REF", "ALT", "QUAL", "FILTER", "INFO", "FORMAT", sample])] # fmt: skip + return "\n".join(h) + "\n" + + +def call_sample_regions(d: npt.NDArray[np.float64], model: PgxCnModel, passing: bool) -> list[dict[str, Any]]: + """PCN-002/009: one result per region.""" + vals = region_values(d, model) + out = [] + for r in model.regions: + x, nf = vals[r.name] + n_all = len(model.region_sel(r.name)) + fit = model.fits.get(r.name) + filters: list[str] = [] + if fit is None or n_all < MIN_REGION_PROBES or nf < MIN_FINITE_FRAC * n_all or not math.isfinite(x): + cn, qual, nr = None, 0, math.nan + filters.append("NoCall") + else: + cn, qual, _ = call_cn(x, fit) + # NR ≈ CN / 2: the continuous CN that the region's mixture means imply for x [OURS] + nr = max(0.0, 3.0 * 2.0 ** ((x - fit.delta) / fit.gamma) - 1.0) / 2.0 if cn else 0.0 + if qual < 7: + filters.append("Q7") + if not passing: + filters.append("SampleQuality") + out.append( + {"region": r, "cn": cn, "qual": qual, "nr": nr, "x": x, "probes": nf, "filters": filters or ["PASS"]} + ) + return out + + +def run_call(o: PgxCallOptions, rep: Reporter) -> int: + """``simurg pgx copy-number call`` (CLI-031).""" + from simurg.core.atomic import atomic_write_text, check_outputs_absent + from simurg.formats.discovery import discover_gtcs + from simurg.formats.fasta import ContigMap, Genome + from simurg.formats.gtc import read_gtc + from simurg.pipeline.cn import write_vcf + from simurg.pipeline.qc import write_tsv + + model = read_model(o.model) + genome = Genome(o.fasta) + if genome.build not in ("GRCh38", "unknown", ""): + raise SimurgError("E_BUILD_MISMATCH", f"PGx CN supports GRCh38 only (FASTA is {genome.build})") + cmap = ContigMap([e.name for e in genome.entries], "auto") + gtcs = discover_gtcs(o.gtcs) + ext = ".pgx_cn.vcf.gz" if o.bgzip else ".pgx_cn.vcf" + check_outputs_absent([o.out / f"{k}{ext}" for k in gtcs] + [o.out / "pgx_cn_batch_summary.tsv"], force=o.force) + n_probes = int(model.header["n_probes"]) + mcl = model.header.get("cluster_file_name") + # pass 1: batch QC from GTC headers (PCN-004) + heads = {} + for k, p in gtcs.items(): + g = read_gtc(p) + heads[k] = (g.logr_dev, g.sample_name or k, g.cluster_file, len(g.logr_ratios)) + n_pass = sum(1 for v in heads.values() if v[0] < o.lrr_sd_max) + if len(heads) < o.min_batch or n_pass < o.min_passing: + raise SimurgError("E_PGX_BATCH_TOO_SMALL", f"{n_pass} sample(s) passed QC. Requires at least " + f"{max(o.min_passing, 1)} samples to proceed.") # fmt: skip + if n_pass < o.recommended_passing: + rep.warn("W_PGX_BATCH_BELOW_RECOMMENDED", f"{n_pass} passing sample(s); {o.recommended_passing} recommended") + for k, (_, _, cf, nl) in heads.items(): + if nl != n_probes: + raise SimurgError("E_MODEL_MANIFEST_MISMATCH", f"{k} has {nl} loci, the model's manifest {n_probes}") + if isinstance(mcl, str) and cf and cf != mcl: + msg = f"Cluster file ({cf}) differs from the CN model's cluster file ({mcl})" + if o.cluster_mismatch == "error": + raise SimurgError("E_MODEL_CLUSTER_MISMATCH", msg, hint="retrain the model (PCN-005)") + rep.warn("W_MODEL_CLUSTER_MISMATCH", msg) + break + contigs: list[tuple[str, int]] = [] + contig_name: dict[str, str | None] = {} + for r in model.regions: + if r.chrom in contig_name: + continue + c = cmap.get(r.chrom) + contig_name[r.chrom] = c + if c is None: + rep.warn("W_PGX_NO_ALT_CONTIG", f"{r.chrom} is not in {o.fasta.name}; REF written as N (IN-062)") + contigs.append((r.chrom, genome.contig_length(c) if c else r.end)) + order = {c: i for i, (c, _) in enumerate(contigs)} + model_sha = hashlib.sha256(o.model.read_bytes()).hexdigest() + mh = model.header + meta = { + "fasta": str(o.fasta.resolve()), + "manifest": f"{mh['manifest_name']};sha256={mh['probe_set_sha256']}", + "model": f"{o.model.name};sha256={model_sha};cluster_file={mcl};regions={mh['regions_bed_version']}", + "batch": f"n_samples={len(heads)};n_passing={n_pass};lrr_sd_max={o.lrr_sd_max}", + } + centre_idx = model.probe_idx[model.region_sel("")] # the scale probes (autosomal), as in training + o.out.mkdir(parents=True, exist_ok=True) + srows, per_region = [], {r.name: {"n_gain": 0, "n_loss": 0, "n_neutral": 0, "n_nocall": 0, "n_filtered": 0} + for r in model.regions} # fmt: skip + failed = 0 + for k, p in gtcs.items(): + try: + g = read_gtc(p) + d = _centred(g.logr_ratios, centre_idx)[model.probe_idx] - model.probe_ref + passing = heads[k][0] < o.lrr_sd_max + res = call_sample_regions(d, model, passing) + except SimurgError: + raise + except Exception as e: # PCN: one sample never stops the batch + failed += 1 + rep.error("E_PGX_CN_FAILED", f"{k}: {e}", sample_id=k) + continue + called = [(x["region"], x["cn"]) for x in res if x["cn"] is not None] + tot = sum(r.end - r.start0 for r, _ in called) + ploidy = sum((r.end - r.start0) * min(cn, CN_MAX) for r, cn in called) / tot if tot else math.nan + meta_s = {**meta, "ploidy": f"{ploidy:.3f}" if math.isfinite(ploidy) else "."} + recs = [] + for x in sorted(res, key=lambda y: (order[y["region"].chrom], y["region"].start0)): + r = x["region"] + c = contig_name.get(r.chrom) + ref = (genome.fetch(c, r.start0 - 1, r.start0) or "N").upper() if c and r.start0 >= 1 else "N" + s, e_ = r.start0 + 1, r.end + cn = "." if x["cn"] is None else str(x["cn"]) + nr = "." if not math.isfinite(x["nr"]) else f"{x['nr']:.6f}" + line = "\t".join([r.chrom, str(s - 1), f"CNV:{r.name}:{r.chrom}:{s}-{e_}", ref, "", str(x["qual"]), + ";".join(x["filters"]), f"END={e_};SVTYPE=CNV;CNVLEN={e_ - s + 1};PROBE={x['probes']};" + f"REGION={r.name};GENE={r.gene}", "CN:NR", f"{cn}:{nr}"]) # fmt: skip + recs.append((r.chrom, s - 2, e_, line)) + pr = per_region[r.name] + if x["cn"] is None: + pr["n_nocall"] += 1 + elif x["cn"] > 2: + pr["n_gain"] += 1 + elif x["cn"] < 2: + pr["n_loss"] += 1 + else: + pr["n_neutral"] += 1 + pr["n_filtered"] += x["filters"] != ["PASS"] + write_vcf(o.out / f"{k}{ext}", _vcf_header(heads[k][1], contigs, meta_s), recs, o.bgzip) + srows.append({"sample_id": k, "sample_name": heads[k][1], "sample_folder": str(Path(p).parent), + "logr_dev": f"{heads[k][0]:.6f}", "passed_batch_qc": str(passing).lower(), + "used_in_reference": "false", + **{f"cn_{x['region'].name}": ("." if x["cn"] is None else x["cn"]) for x in res}, + "flags": ""}) # fmt: skip + cols = ["region", "gene", "n_gain", "n_loss", "n_neutral", "n_nocall", "n_filtered"] + rows: list[dict[str, Any]] = [ + {"region": r.name, "gene": r.gene, **per_region[r.name]} for r in sorted(model.regions, key=lambda r: r.name) + ] + totals = {c: sum(per_region[r.name][c] for r in model.regions) for c in cols[2:]} + rows.append({"region": "ALL", "gene": "", **totals}) + write_tsv(o.out / "pgx_cn_batch_summary.tsv", cols, rows) + scols = ["sample_id", "sample_name", "sample_folder", "logr_dev", "passed_batch_qc", "used_in_reference"] + scols += [f"cn_{r.name}" for r in model.regions] + ["flags"] + write_tsv(o.out / "pgx_cn_sample_summary.tsv", scols, srows) + run = {"software": f"SimurgArray {__version__}", **meta, "samples": len(gtcs), "failed": failed} + atomic_write_text(o.out / "pgx_cn_run.json", json.dumps(run, indent=1)) + rep.info(f"PGx CN: {len(srows)}/{len(gtcs)} sample(s) called") + if failed == 0: + return ExitCode.OK + return ExitCode.PARTIAL if srows else ExitCode.ALL_FAILED diff --git a/tests/unit/test_pgx_cn.py b/tests/unit/test_pgx_cn.py new file mode 100644 index 0000000..471327e --- /dev/null +++ b/tests/unit/test_pgx_cn.py @@ -0,0 +1,145 @@ +"""Targeted PGx copy number (spec PCN-001 … PCN-017, CLI-031/032, IN-050).""" + +from __future__ import annotations + +import gzip +from pathlib import Path + +import numpy as np +from typer.testing import CliRunner + +from simurg.cli.main import app +from simurg.formats.gtc import Gtc, write_gtc +from simurg.pgx.cn import TEMPLATE, call_cn, fit_region, hwe_weights, parse_regions +from simurg.pipeline.pgx_cn import default_regions +from simurg.synth.files import SynthLocus, encode_bpm + + +def _cluster_values(rng: np.random.Generator, cns: np.ndarray, delta: float, gamma: float, sd: float, + cn0_drop: float = 0.0) -> np.ndarray: # fmt: skip + x = delta + gamma * TEMPLATE[cns] + rng.normal(0, sd, len(cns)) + return np.where(cns == 0, x + cn0_drop, x) + + +def _hwe_cns(rng: np.random.Generator, n: int, f: list[float]) -> np.ndarray: + return rng.choice(len(f), n, p=f) + rng.choice(len(f), n, p=f) + + +def test_hwe_weights() -> None: # PCN-009 [OURS] + w = hwe_weights(np.array([0.3, 0.6, 0.1, 0.0])) + assert abs(w.sum() - 1) < 1e-12 + assert np.allclose(w[:4], [0.09, 0.36, 0.06 + 0.36, 0.12]) + + +def test_labels_identified() -> None: # PCN-009: GSTM1-like, GSTT1-like and one-cluster regions + rng = np.random.default_rng(1) + cns = _hwe_cns(rng, 60, [0.7, 0.3]) # common deletion: mostly CN 0 and 1 + x = _cluster_values(rng, cns, 0.1, 0.22, 0.02) + f = fit_region(x) + assert np.mean([call_cn(v, f)[0] for v in x] == cns) > 0.95 + cns = _hwe_cns(rng, 60, [0.45, 0.55]) + x = _cluster_values(rng, cns, 0.3, 0.55, 0.02, cn0_drop=-1.2) # CN 0 far below the template + f = fit_region(x) + assert f.cn0_shift < -0.5 + assert np.mean([call_cn(v, f)[0] for v in x] == cns) > 0.95 + cns = np.full(40, 2) + cns[:2] = 1 + x = _cluster_values(rng, cns, -0.05, 0.25, 0.015) + f = fit_region(x) + assert [call_cn(v, f)[0] for v in x] == cns.tolist() + + +def test_packaged_regions() -> None: # PCN-001 + name, text = default_regions() + regs = parse_regions(text) + assert name.endswith(".v1.bed") and len(regs) == 10 + assert {r.gene for r in regs} == {"CYP2D6", "CYP2A6", "GSTM1", "GSTT1", "UGT2B17", "SULT1A1", "CYP2E1"} + assert next(r for r in regs if r.gene == "GSTT1").chrom == "chr22_KI270879v1_alt" + + +def _dataset(tmp: Path, n_samples: int = 40) -> tuple[Path, Path, Path, Path]: + rng = np.random.default_rng(5) + loci = [] + for c in ("1", "2"): + for k in range(1500): + loci.append(SynthLocus(f"cnv{c}_{k}", c, 2000 * (k + 1), "[A/G]", len(loci) + 1, 0, 0)) + bpm = tmp / "m.bpm" + bpm.write_bytes(encode_bpm(loci)) + fa = tmp / "g.fa" + body = {c: "".join(rng.choice(list("ACGT"), 1500 * 2000 + 10_000)) for c in ("1", "2")} + fa.write_text("".join(f">{c}\n" + "\n".join(s[i : i + 60] for i in range(0, len(s), 60)) + "\n" + for c, s in body.items())) # fmt: skip + regions = tmp / "regions.bed" + regions.write_text("# test regions\n1\t100000\t160000\tGENEA\tGENEA\n2\t300000\t340000\tGENEB\tGENEB\n" + "3\t1000\t5000\tGENEC\tGENEC\n") # fmt: skip + n = len(loci) + pos = np.array([lc.pos for lc in loci]) + chrom = np.array([lc.chrom for lc in loci]) + in_a = (chrom == "1") & (pos > 100_000) & (pos <= 160_000) + in_b = (chrom == "2") & (pos > 300_000) & (pos <= 340_000) + cn_a = _hwe_cns(rng, n_samples, [0.5, 0.5]) + gtcs = tmp / "gtcs" + truth = {} + for i in range(n_samples): + sid = f"2000000000{i:02d}_R01C01" + lrr = rng.normal(0, 0.08, n) + lrr[in_a] += 0.24 * TEMPLATE[cn_a[i]] + (-0.2 if cn_a[i] == 0 else 0.0) + cn_b = 3 if i == 0 else 2 + lrr[in_b] += 0.25 * TEMPLATE[cn_b] + truth[sid] = (int(cn_a[i]), cn_b) + g = Gtc(num_snps=n, genotypes=np.ones(n, np.uint8), genotype_scores=np.full(n, 0.9, np.float32), + b_allele_freqs=np.zeros(n, np.float32), logr_ratios=lrr.astype(np.float32), + raw_x=np.full(n, 1000, np.uint16), raw_y=np.full(n, 1000, np.uint16), + base_calls=np.zeros((n, 2), np.uint8), gender=b"F", call_rate=0.995, logr_dev=0.12, + cluster_file="c.egt") # fmt: skip + (gtcs / sid).mkdir(parents=True) + write_gtc(g, gtcs / sid / f"{sid}.gtc") + (tmp / "truth.txt").write_text("\n".join(f"{k}\t{a}\t{b}" for k, (a, b) in truth.items())) + return bpm, fa, gtcs, regions + + +def test_pgx_cn_train_and_call(tmp_path: Path) -> None: # CLI-031/032, PCN-002/004/011/012/013, IN-050 + bpm, fa, gtcs, regions = _dataset(tmp_path) + run = CliRunner() + r = run.invoke(app, ["pgx", "copy-number", "train", "--gtcs", str(gtcs), "--bpm", str(bpm), "--fasta", str(fa), + "--out", str(tmp_path / "model"), "--regions", str(regions), "--quiet"]) # fmt: skip + assert r.exit_code == 2 and "E_PGX_TRAIN_TOO_SMALL" in r.output # 40 < 96 (PCN-005) + r = run.invoke(app, ["pgx", "copy-number", "train", "--gtcs", str(gtcs), "--bpm", str(bpm), "--fasta", str(fa), + "--out", str(tmp_path / "model"), "--regions", str(regions), "--min-samples", "30", + "--min-passing", "30", "--quiet"]) # fmt: skip + assert r.exit_code == 0, r.output + model = tmp_path / "model" / "m.simurg-pgxcn" + assert model.exists() and (tmp_path / "model" / "train_exclusions.tsv").exists() + out = tmp_path / "call" + r = run.invoke(app, ["pgx", "copy-number", "call", "--gtcs", str(gtcs), "--model", str(model), "--fasta", str(fa), + "--out", str(out), "--min-batch", "50", "--quiet"]) # fmt: skip + assert r.exit_code == 2 and "E_PGX_BATCH_TOO_SMALL" in r.output + r = run.invoke(app, ["pgx", "copy-number", "call", "--gtcs", str(gtcs), "--model", str(model), "--fasta", str(fa), + "--out", str(out), "--quiet"]) # fmt: skip + assert r.exit_code == 0, r.output + truth = {ln.split("\t")[0]: (int(ln.split("\t")[1]), int(ln.split("\t")[2])) + for ln in (tmp_path / "truth.txt").read_text().splitlines()} # fmt: skip + ok_a = ok_b = 0 + for sid, (a, b) in truth.items(): + with gzip.open(out / f"{sid}.pgx_cn.vcf.gz", "rt") as fh: + text = fh.read() + assert (out / f"{sid}.pgx_cn.vcf.gz.tbi").exists() + recs = [ln.split("\t") for ln in text.splitlines() if not ln.startswith("#")] + assert [x[7].split("REGION=")[1].split(";")[0] for x in recs] == ["GENEA", "GENEB", "GENEC"] # PCN-002 + for x in recs: # PCN-011: POS = S − 1, ID coordinates = POS + 1 … END + end = int(x[7].split("END=")[1].split(";")[0]) + s_, e_ = x[2].split(":")[-1].split("-") + assert int(s_) == int(x[1]) + 1 and int(e_) == end and x[8] == "CN:NR" + genec = recs[2] + assert genec[6] == "NoCall" and genec[9].startswith(".:") + ok_a += int(recs[0][9].split(":")[0]) == a + ok_b += int(recs[1][9].split(":")[0]) == b + assert ok_a >= 0.95 * len(truth) and ok_b >= 0.9 * len(truth), (ok_a, ok_b) + rows = [ln.split("\t") for ln in (out / "pgx_cn_batch_summary.tsv").read_text().splitlines()] + assert rows[0] == ["region", "gene", "n_gain", "n_loss", "n_neutral", "n_nocall", "n_filtered"] + assert [r_[0] for r_ in rows[1:]] == ["GENEA", "GENEB", "GENEC", "ALL"] + for r_ in rows[1:-1]: # PCN-012: the four counts add up to N + assert sum(int(v) for v in r_[2:6]) == len(truth) + head = (out / "pgx_cn_sample_summary.tsv").read_text().splitlines()[0].split("\t") + assert head[:6] == ["sample_id", "sample_name", "sample_folder", "logr_dev", "passed_batch_qc", + "used_in_reference"] and head[-1] == "flags" # fmt: skip