Release 2.0.0: v2 modules (ArrayTrain, CN/LOH, PGx, methylation QC)
ci / test (3.13) (push) Failing after 30s
ci / test (3.12) (push) Failing after 31s

User guide: ArrayTrain section, v2 commands in the table, v2 accuracy rows,
updated limitations. README status for v2.0.0. Version 2.0.0 (Python package);
the Rust kernels are unchanged and stay simurg_core 1.1.0. uv.lock committed
(uv lock --check clean).

Full suite on the Pi, file by file under the memory guard: 29 files exit 0
(unit: 282 passed), ruff and mypy clean.

Co-Authored-By: Claude Opus 5.5 <[email protected]>
Claude-Session: https://claude.ai/code/session_01814GKVzZVL3FuF3DeiZXhv
This commit is contained in:
kemalandClaude Opus 5.5 committed 2026-09-30 06:24:12 +01:00
1 parent ca770b12d6
commit 6a483cdf5c
6 files changed
+805 -19

No files matched your search

+47 -6
View File
@@ -1,4 +1,4 @@
# SimurgArray user guide (v1.1)
# SimurgArray user guide (v2.0)
SimurgArray turns Illumina Infinium genotyping BeadChip scans (IDAT files) into genotype calls, QC tables and
reports, VCF/BCF and PLINK files. It is a personal, non-commercial **research** tool — not for clinical or
@@ -6,6 +6,8 @@ diagnostic use. Behaviour is specified item by item in `SIMURGARRAY_FUNCTIONAL_S
cited below); the algorithms and their sources are in `SIMURGARRAY_MASTER_PLAN.md`.
Contents: [Install](#install) · [Inputs](#inputs) · [A complete run](#a-complete-run) · [Commands](#commands) ·
[Methylation](#methylation-arrays) · [CN/LOH](#copy-number-loh-and-mosaicism) · [PGx](#pharmacogenomics) ·
[Cluster training](#cluster-training-arraytrain) ·
[Presets and reproducibility](#presets-and-reproducibility) · [QC](#quality-control) · [Outputs](#outputs) ·
[Accuracy](#accuracy) · [Limitations](#limitations) · [Troubleshooting](#troubleshooting)
@@ -33,15 +35,15 @@ Without them the same code runs in pure Python/NumPy; `SIMURG_PURE_PYTHON=1` for
`manylinux_2_28` kernel wheels for x86_64 and aarch64, glibc ≥ 2.28):
```bash
python3 -m venv simurg-env && simurg-env/bin/pip install simurg-1.1.0-py3-none-any.whl \
python3 -m venv simurg-env && simurg-env/bin/pip install simurg-2.0.0-py3-none-any.whl \
simurg_core-1.1.0-cp312-abi3-manylinux_2_28_$(uname -m).whl # the second wheel is optional
```
**Docker** (image ≈ 550 MB, kernels included, runs as a non-root user; data stay in the mounted folder):
```bash
docker build -t simurg:1.1.0 .
docker run --rm --user "$(id -u):$(id -g)" -v "$PWD":/work simurg:1.1.0 \
docker build -t simurg:2.0.0 .
docker run --rm --user "$(id -u):$(id -g)" -v "$PWD":/work simurg:2.0.0 \
genotype --bpm M.bpm --egt C.egt --idats idats --out run/gtc
```
@@ -108,6 +110,10 @@ simurg export finalreport --bpm M.bpm --csv M.csv --egt C.egt --gtcs run/gtc --o
| `simurg qc suggest-thresholds` | data-driven threshold YAML from a cohort's `controls_qc.tsv` |
| `simurg gtc-to-vcf` | harmonised VCF/BCF. Options: `--merge`, `--unsquash-duplicates`, `--merge-snv-indel`, `--auxiliary-loci`, `--info-fields`, `--extra-format`, `--filter-gs`, `--drop-ambiguous`, `--ploidy-policy sex-aware\|diploid`, `--include-loci/--exclude-loci` |
| `simurg export plink\|bedgraph\|finalreport\|locus-summary` | downstream formats |
| `simurg methyl qc` / `methyl convert-manifest` | methylation-array control QC, detection p-values, β/M values (see below) |
| `simurg cn build-model\|call\|annotate` | genome-wide CNV / LOH / mosaic calling from GTCs (see below) |
| `simurg pgx db build\|coverage`, `pgx copy-number train\|call`, `pgx star-allele call`, `pgx phenotype` | pharmacogenomics (see below) |
| `simurg train` | ArrayTrain: an EGT cluster file from a cohort's GTCs (see below) |
| `simurg compare A B` | field-level differential test of GTCs or VCFs (files or folders paired by name): EXACT / NEAR / DISCREPANCY / UNKNOWN with an HTML dashboard; `--fail-on-discrepancy` exits 3 |
Global logging options on every command: `--quiet`, `--log-level`, `--log-json` (JSON Lines on stderr), `--debug`.
@@ -190,6 +196,36 @@ simurg pgx phenotype --calls run/pgx --db pgxdb/pgx_db-….sqlite [--guideline c
* **DPWG**: genes CPIC covers use the CPIC mappings (as ClinPGx does); DPWG-only genes (CYP1A2, CYP3A4, F5) get no
phenotype (K-029).
## Cluster training (ArrayTrain)
`simurg train` builds an EGT cluster file (v3, data block v9) from a cohort's GTCs, so calling does not depend on a
vendor cluster file (plan §12; validation in `docs/validation/V2_arraytrain.md`).
```bash
simurg genotype --bpm GSA.bpm --egt any.egt --idats cohort_idats/ --out cohort_gtc # any EGT: training reads raw X/Y
simurg train --gtcs cohort_gtc --bpm GSA.bpm --out mycohort.egt # ≥ 100 GTCs (--min-samples)
simurg genotype --bpm GSA.bpm --egt mycohort.egt --idats new_idats/ --out run/ # call with the trained file
```
* Training uses each GTC's raw intensities and normalisation transforms, not its genotypes. The cluster file the
cohort was first called with therefore does not matter.
* Samples are chosen by sex: X loci are trained from females, Y from males, and MT is haploid. The sex is the GTC's
`gender` field.
* At most `--max-samples` (2000) GTCs are used, as a seeded random subset.
* **Per locus:**
1. 1-, 2- and 3-cluster θ mixtures are fitted; BIC picks one.
2. Clusters are labelled AA/AB/BB by Hardy–Weinberg fit plus weak position priors, so compressed or asymmetric
loci are handled.
3. Missing clusters are imputed from relations learned across well-clustered loci of the same assay type.
* **Outputs:** `<out>.egt` and `<out>.egt.json`. The JSON gives the counts of 3-, 2- and 1-cluster loci and the
median training score.
* GenTrain and cluster-separation scores in a trained file are SimurgArray's own. GenCall scores are therefore not
comparable with calls made with a vendor cluster file, so re-derive any GenCall-score thresholds (K-030).
* **Accuracy.** GSE165845, trained on 120 arrays and evaluated on 240 other arrays of the same individuals, chr22
against 1000G 30x. The trained file gives 99.44 % concordance at a 99.61 % call rate; the vendor file gives
99.46 % at 99.15 %.
* **Time.** 120 GSA samples take about 70 min on a Raspberry Pi 5, with a 2 GB peak.
## Presets and reproducibility
* `--preset current` (default): normalization 1.2.0, sex rule 2, GTC v5 — the current vendor software behaviour.
@@ -250,13 +286,18 @@ Validated phase by phase (reports in `docs/validation/`):
| Pedigree QC, 2141 Omni2.5 samples | 717/717 parent–offspring pairs; trio Mendel error median 0.031 % |
| 360 GSA arrays (GSE165845) vs NYGC 30x WGS | genotype concordance 0.9934; replicate concordance 0.9966 |
| NA12878 × 3 vs GIAB HG001 v4.2.1 | concordance 0.9985; non-reference F1 0.997 |
| Methylation QC vs sesame (minfiDataEPIC) | β r = 0.99998 |
| CN, 360 GSE165845 arrays vs NYGC 30x SV | PASS precision DEL 44 % / DUP 57 %; rare-DEL sensitivity 30 % |
| PGx, GSA v4 demo vs GeT-RM (16 CEPH) | 154/168 gene calls; GSE165845 array-resolvable 95.7 % |
| ArrayTrain, GSE165845 chr22 vs 1000G 30x | 99.44 % concordance (vendor EGT 99.46 %) at 99.61 % call rate |
A nightly differential run (`scripts/nightly.py`) re-checks these against the oracles.
## Limitations
* Calls need an existing EGT cluster file; cluster training (v2), copy number / LOH / mosaicism and PGx
diplotyping are not in v1.0.
* ArrayTrain was validated on one chip and cohort (GSA v1, GSE165845). CN calling finds rare events, but common
CNVs are normalised away, and mosaic losses below ~20 % are not separable from noise (K-027). PGx is limited
to what the array's probes can resolve: `simurg pgx coverage` lists it. All v2 modules are research-grade.
* The largest cohort tested is 2141 samples (plan goal ≥ 10 000). Memory grows with samples × loci in the cohort
QC steps; on an 8 GB machine keep `--threads` low (≈ 2–3 GB per worker for GSA/Omni-size arrays).
* Speed: with the compiled kernels a GSA sample takes ~1.1 s on a Raspberry Pi 5 core (pure Python ~3.4 s); see