318 lines
22 KiB
Markdown
318 lines
22 KiB
Markdown
<img src="../branding/logo/png/lockup_horizontal_light.png" alt="SimurgArray" width="360">
|
||
|
||
# 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
|
||
diagnostic use. Behaviour is specified item by item in `SIMURGARRAY_FUNCTIONAL_SPEC.md` (IDs such as CLI-021 are
|
||
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)
|
||
|
||
## Install
|
||
|
||
Python ≥ 3.12. With [uv](https://docs.astral.sh/uv/):
|
||
|
||
```bash
|
||
git clone <repository> simurg && cd simurg
|
||
uv venv && uv pip install -e . # add ".[dev]" for tests and linters
|
||
.venv/bin/simurg --version
|
||
```
|
||
|
||
**Optional compiled kernels** (≈ 3–4× faster genotyping, 4× faster manifest loading; results are byte-identical):
|
||
|
||
```bash
|
||
curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh -s -- -y --profile minimal # Rust (needs a C linker, e.g. gcc)
|
||
uv pip install ./rust/simurg-core
|
||
.venv/bin/python -c "from simurg._core import backend; print(backend())" # rust 1.1.0
|
||
```
|
||
|
||
Without them the same code runs in pure Python/NumPy; `SIMURG_PURE_PYTHON=1` forces that path.
|
||
|
||
**From release files** (`bash scripts/build_dist.sh` writes them to `dist/`: a pure wheel, an sdist and
|
||
`manylinux_2_28` kernel wheels for x86_64 and aarch64, glibc ≥ 2.28):
|
||
|
||
```bash
|
||
python3 -m venv simurg-env && simurg-env/bin/pip install simurgarray-2.0.0-py3-none-any.whl \
|
||
simurgarray_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: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
|
||
```
|
||
|
||
**conda** (recipes in `packaging/conda/`):
|
||
|
||
```bash
|
||
conda build -c conda-forge --output-folder conda-bld packaging/conda/simurg-core packaging/conda/simurg
|
||
conda create -n simurgarray -c ./conda-bld -c conda-forge simurgarray simurgarray-core
|
||
```
|
||
|
||
All four installs give identical results on x86_64 and aarch64 (`docs/validation/P10.md`).
|
||
Only pure-Python dependencies are needed (numpy, pyarrow, typer, rich, pyyaml, openpyxl). VCF/BCF, tabix/CSI
|
||
indexes and PLINK files are written by SimurgArray itself; bcftools/plink2 are not required (they are used only
|
||
as test oracles).
|
||
|
||
## Inputs
|
||
|
||
| File | Needed by | Notes |
|
||
|---|---|---|
|
||
| **IDAT pairs** (`<barcode>_<RxxCxx>_Grn.idat` + `_Red.idat`, optionally `.gz`) | `genotype`, `ingest`, `qc controls` | `--idats DIR` is searched recursively; the sample ID is `<barcode>_<position>`. Encrypted (v1) IDATs are refused. Duplicate pairs: `--on-duplicate-idat error\|first\|newest` |
|
||
| **BPM manifest** | all commands | defines loci, bead addresses, normalization bins, positions (its genome build decides the VCF build) |
|
||
| **EGT cluster file** | `genotype`, `qc loci`, optional elsewhere | must describe the manifest's loci (`E_EGT_MANIFEST_MISMATCH` otherwise); a cluster file of a sibling manifest version works if it covers the loci (a warning is printed) |
|
||
| **CSV manifest** | `gtc-to-vcf` (indels), `export finalreport` (strands); optional for `genotype` (cross-check) | same product/version as the BPM |
|
||
| **Reference FASTA** | `gtc-to-vcf` | plain or bgzip; contig names must match the manifest (`1…22, X, Y, MT`) or use `--contig-alias auto\|FILE`. An mmap cache is built next to it on first use (~3 GB for GRCh38) |
|
||
| **Sample sheet** (CSV) | optional everywhere | GenomeStudio style (`[Header]`, `[Manifests]`, `[Data]`) or a plain table. Key columns `SentrixBarcode_A`, `SentrixPosition_A`; recognised: `Sample_ID`, `Sample_Name`, `Sample_Plate`, `Sample_Well`, `Sample_Group`, `Sex`/`Gender`, `Parent1`, `Parent2`, `Replicate`. Other columns are carried into the QC table. Example: `examples/sample_sheet.csv` |
|
||
|
||
`simurg inspect FILE…` prints the header of any of these (add `--json` for scripts).
|
||
|
||
## A complete run
|
||
|
||
`examples/gsa_demo.sh` runs everything below on the public 36-sample GSA v4 PGx demo
|
||
(`bash scripts/fetch_data.sh p0 && bash scripts/fetch_data.sh p6` first).
|
||
|
||
```bash
|
||
# 1. IDAT → GTC v5 per sample + sample_summary.tsv + control-probe QC + contamination + relatedness
|
||
simurg genotype --bpm M.bpm --egt C.egt --csv M.csv --idats idats/ --sheet sheet.csv \
|
||
--out run/gtc --qc-profile genotyping --threads 0
|
||
|
||
# 2. self-contained HTML QC report (+ xlsx/csv QC table)
|
||
simurg qc report --data run/gtc --out run/qc_report --output-format both
|
||
|
||
# 3. plus-strand VCF.gz + .tbi per sample, or one cohort BCF + .csi (--merge)
|
||
simurg gtc-to-vcf --bpm M.bpm --csv M.csv --fasta GRCh38.fa --gtcs run/gtc --out run/vcf
|
||
simurg gtc-to-vcf --bpm M.bpm --csv M.csv --fasta GRCh38.fa --gtcs run/gtc --out run/bcf --merge
|
||
|
||
# 4. exports
|
||
simurg export plink --vcfs run/vcf --gtcs run/gtc --out run/plink --format both # .bed/.bim/.fam + .pgen/.pvar/.psam
|
||
simurg export bedgraph --bpm M.bpm --gtcs run/gtc --out run/tracks --chrom-style ucsc
|
||
simurg export finalreport --bpm M.bpm --csv M.csv --egt C.egt --gtcs run/gtc --out run/report # --matrix for loci × samples
|
||
```
|
||
|
||
`--threads 0` uses all cores. Every command refuses to overwrite existing outputs unless `--force` is given.
|
||
|
||
## Commands
|
||
|
||
| Command | Purpose |
|
||
|---|---|
|
||
| `simurg inspect FILES… [--json]` | headers / tables of contents of IDAT, BPM, CSV, EGT, GTC, sample sheets (CLI-020) |
|
||
| `simurg genotype` | IDAT → normalization → GenCall-equivalent calls → GTC; sample summary, controls, cohort QC. Key options: `--cutoff` (GenCall score no-call threshold, 0.15), `--preset`, `--norm-version`, `--sex-version`, `--gtc-version 3\|5`, `--autocall-date`, `--no-controls`, `--no-cohort-qc`, QC threshold options |
|
||
| `simurg ingest` / `simurg normalize` | the first two stages on their own: raw per-sample intensities (Parquet), then fitted transforms and normalized X/Y/θ/R. `ingest --resume` keeps outputs whose inputs are unchanged |
|
||
| `simurg qc controls` | control-probe metrics from IDATs without genotyping (`--array-type genotyping\|methylation`) |
|
||
| `simurg qc samples\|contamination\|relatedness\|loci` | sample QC, BAF-regression contamination, KING-robust kinship (+ Mendel errors from `Parent1/Parent2`), per-locus summary — from existing GTCs |
|
||
| `simurg qc report` | HTML report (views, linked selection, thresholds editor) and QC table from one or more dataset folders |
|
||
| `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`.
|
||
Run `simurg COMMAND --help` for the full list.
|
||
|
||
## Methylation arrays
|
||
|
||
`simurg methyl qc` runs control-probe QC and a signal summary for Infinium methylation arrays (EPIC v1/v2, MSA,
|
||
custom designs) from IDATs and the methylation CSV manifest (spec MET-001 … MET-018). It is a QC tool: there is no
|
||
differential-methylation analysis.
|
||
|
||
```bash
|
||
simurg methyl qc --manifest infinium-methylationepic-v-1-0-b5-manifest-file.csv --idats idats/ \
|
||
--sheet SampleSheet.csv --out run/methyl --profile auto
|
||
```
|
||
|
||
* 21 control metrics (staining, extension, hybridization, target removal I/II, bisulfite conversion I/II with their
|
||
background ratios, specificity I/II, non-polymorphic, restoration); thresholds from the `methylation_epic` /
|
||
`methylation_msa` profiles (`auto` picks by manifest name), `--config` or `--threshold KEY=VALUE`.
|
||
* Per sample: detection p-values (pOOBAH: out-of-band Infinium I signal and negative controls), noob background
|
||
correction, dye-bias correction (`--dye-bias linear|nonlinear`), β = M/(M+U+α) (`--beta-offset`, default 0) and
|
||
M-values.
|
||
* `pass_qc` = every enabled, evaluable metric passes; `--min-prop-passing X` also requires that fraction of probes
|
||
to pass detection (`--pval-threshold`, default 0.05).
|
||
* Outputs: `methyl_qc_summary.tsv/.xlsx/.html` (detection histogram and a β PCA coloured by `Sample_Group`),
|
||
`samples/<id>/<id>_cgs.tsv.gz` (probe_id, beta, m_value, detection_p), `<id>_ctrl.tsv.gz`, `pcs.tsv.gz`,
|
||
`probes.tsv`/`controls.tsv` (the manifest in the open-source toolkits' layout; also written by
|
||
`simurg methyl convert-manifest`), and `controls_raw.tsv`/`controls_qc.tsv`, so `simurg qc report` accepts the folder.
|
||
* Sample sheets may use `Sentrix_ID`/`Sentrix_Position` (minfi/sesame convention) as key columns.
|
||
* `simurg qc controls --array-type methylation --csv MANIFEST.csv` computes only the control tables.
|
||
|
||
## Copy number, LOH and mosaicism
|
||
|
||
Genome-wide CNV / LOH / mosaic calling from GTCs (spec CYTO-001 … CYTO-042; validation in
|
||
`docs/validation/V2_cn.md`).
|
||
|
||
```bash
|
||
# once per manifest + batch: GC covariates, population BAF, and a batch reference from ≥ 24 normal GTCs
|
||
simurg cn build-model --bpm GSA.bpm --fasta GRCh38.fa --egt GSA.egt --reference-gtcs normals/ --out gsa.simurg-cyto
|
||
simurg cn call --bpm GSA.bpm --model gsa.simurg-cyto --gtcs run/ --fasta GRCh38.fa --out run/cn
|
||
simurg cn annotate --vcfs run/cn --annotation-bundle simurg-cyto-annot-GRCh38-….zip # ISCN, cytobands, genes
|
||
```
|
||
|
||
* Per sample `<id>.cyto.cnv.vcf.gz` (+ `.tbi`) with DEL/DUP/LOH/mosaic calls, `<id>.cn_segments.parquet`, and
|
||
`cn_sample_summary.tsv` / `cn_call_manifest.json` for the batch. `scripts/build_annotation_bundle.sh` builds the
|
||
annotation bundle from UCSC tables.
|
||
* Mosaic fractions are estimated from the noise-normalised BAF split (CYTO-041); losses without heterozygotes are
|
||
promoted to constitutional deletions. Per-probe LRR noise from the reference is on by default
|
||
(`--no-probe-lrr-sd` turns it off).
|
||
* Array CN calling finds rare events; CNVs common in the reference batch are normalised away (K-027).
|
||
|
||
## Pharmacogenomics
|
||
|
||
Targeted copy number, star alleles and phenotypes (spec PCN-001 … PCN-017, PGX-001 … PGX-035; validation in
|
||
`docs/validation/V2_pgx_cn.md` and `V2_pgx_star.md`). GRCh38 only.
|
||
|
||
```bash
|
||
bash scripts/pgx/fetch_sources.sh # CPIC (CC0), PharmVar, ClinPGx → data/pgx/sources/
|
||
simurg pgx db build --sources data/pgx/sources --fasta GRCh38.fa --manifests GSA-PGx.bpm --csv GSA-PGx.csv \
|
||
[--mask GSA-PGx.mask.tsv] --out pgxdb/ # → pgx_db-1-<sha12>.sqlite
|
||
simurg pgx coverage --bpm GSA-PGx.bpm --db pgxdb/pgx_db-….sqlite # which alleles this array can call
|
||
simurg pgx copy-number train --gtcs cohort/ --bpm GSA-PGx.bpm --csv GSA-PGx.csv --fasta GRCh38.fa --out cnmodel/
|
||
simurg pgx copy-number call --gtcs run/ --model cnmodel/GSA-PGx.simurg-pgxcn --fasta GRCh38.fa --out run/vcf
|
||
simurg gtc-to-vcf --gtcs run/ --bpm GSA-PGx.bpm --csv GSA-PGx.csv --fasta GRCh38.fa --out run/vcf
|
||
simurg pgx star-allele call --vcfs run/vcf --db pgxdb/pgx_db-….sqlite --out run/pgx
|
||
simurg pgx phenotype --calls run/pgx --db pgxdb/pgx_db-….sqlite [--guideline cpic|dpwg]
|
||
```
|
||
|
||
* **Targeted CN** for CYP2D6 (5′ region, intron 2, exon 9), CYP2A6, GSTM1, GSTT1 (alt contig), UGT2B17, SULT1A1 and
|
||
CYP2E1: one record per region in `<id>.pgx_cn.vcf.gz`. Training needs ≥ 96 samples (`--min-samples` lowers it
|
||
for research use) and the platform from `--csv` or `--platform`.
|
||
* **Star alleles**: CPIC star-allele genes, PharmVar genes CPIC does not define (CYP1A2, CYP2A6, CYP2A13, CYP2C8,
|
||
CYP3A4), G6PD, HLA-A/B tag variants, genotype-only variant-based genes (CPIC and ClinPGx), and `*0` genotypes
|
||
for GSTM1/GSTT1/UGT2B17. CYP2D6 structures (`*5`, `xN`, `*68+*4`, `*36+*10`, `*13`) come from the sub-region CNs.
|
||
* **Probe mask** (PGX-010): `scripts/pgx/derive_masks.py` compares array calls with 1000 Genomes 30x for the batch's
|
||
1000G individuals and masks probes below 95 % concordance.
|
||
* **Outputs**: `diplotype_summary.tsv`, `supporting_variants.tsv`, `unresolved_variants.tsv`,
|
||
`variants_not_covered.tsv`, `pgx_calls.json`; after `phenotype`, `<id>.pgx.json` (`simurg.pgx/1`) and
|
||
`diplotype_annotations.tsv`. Array-equivalent alleles, no-call impact and alternatives are listed per gene.
|
||
* **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.
|
||
* `--preset autocall-1.6.2`: normalization 1.1.2, sex rule 1, GTC v3 with legacy sample statistics — reproduces the
|
||
public 1000 Genomes Omni2.5 GTCs byte for byte (2138 of 2141; the 3 exceptions are in `docs/known_unknowns.md`, K-009).
|
||
* Individual settings override the preset (`--norm-version`, `--sex-version`, `--gtc-version`).
|
||
* The default numerics are bit-exact with the reference implementations (including their float32 and median quirks);
|
||
where a "clean-math" choice exists it is an explicit option, e.g. `qc samples --median exact`.
|
||
* Outputs are deterministic for the same inputs and options. `--autocall-date` fixes the one time-stamp written
|
||
into GTCs. Every run writes `run_manifest.json` (version, git revision, argv, effective options and where each
|
||
came from, environment, seed) and a log under `logs/`.
|
||
|
||
## Quality control
|
||
|
||
`genotype` computes per-sample call rate (all / autosomal), LogR deviation, GC10/GC50, intensity percentiles,
|
||
estimated sex and heterozygosity; control-probe metrics (staining, extension, hybridisation, target removal,
|
||
stringency, non-specific binding, non-polymorphic, restoration); BAF-regression contamination (≥ 10 samples) and
|
||
KING-robust relatedness (duplicates, first/second/third degree, unexpected relatives, Mendel errors when the sheet
|
||
has parents).
|
||
|
||
Each sample gets two statuses: `functional_status` (call rate, LogR deviation, …) and `control_status` (control
|
||
probes), plus a `flags` list (`low_call_rate`, `high_lrr_sd`, `het_outlier`, `duplicate_of:…`, `sex_mismatch`, …).
|
||
Thresholds come from a **profile** (`--qc-profile genotyping|genotyping_example|genotyping_ex|methylation_epic|methylation_msa|auto`,
|
||
files in `src/simurg/qc/profiles/`), a YAML file (`--qc-config`), or single overrides
|
||
(`--threshold-call-rate 0.97`, `--threshold KEY=VALUE`). Without a profile, flags are reported but no sample fails.
|
||
`simurg qc suggest-thresholds` proposes cohort-specific values.
|
||
|
||
## Outputs
|
||
|
||
`simurg genotype --out run/gtc`:
|
||
|
||
| Path | Content |
|
||
|---|---|
|
||
| `samples/<id>/<id>.gtc` | GTC v5 (or v3): calls, scores, raw and normalized intensities, BAF, LRR, transforms, sample statistics — readable by BeadArrayFiles, Picard, gtc2vcf |
|
||
| `sample_summary.tsv` / `.json` | one row per sample: status, call rates, LogR dev, sex (estimated/reported/check), GC10/50, intensity percentiles, het rate, contamination, control status, flags, sheet columns |
|
||
| `controls_raw.tsv`, `controls_qc.tsv` | control-probe intensities and derived metrics |
|
||
| `relatedness.tsv` | related / duplicate pairs: kinship, concordance, degree, expected relation |
|
||
| `mendel.tsv` | trio Mendelian error counts (when the sheet has `Parent1/Parent2`) |
|
||
| `run_manifest.json`, `logs/` | provenance |
|
||
|
||
`simurg gtc-to-vcf --out run/vcf`: `<id>.snv.vcf.gz` + `.tbi` per sample, or `cohort.snv.bcf` + `.csi` with
|
||
`--merge`. FORMAT `GT:GS:BAF:LRR` (GS = GenCall score); alleles on the plus strand of the reference; multi-probe
|
||
designs of one site merged (`--unsquash-duplicates` keeps one record per probe); indels resolved from the CSV
|
||
`SourceSeq`. `audit/loci.parquet` records the strand decision of every locus, `audit/merges.parquet` every merge.
|
||
|
||
Exports: PLINK 1 (`.bed/.bim/.fam`) and PLINK 2 (`.pgen/.pvar/.psam`); gzipped BAF/LRR bedGraph per sample
|
||
(`--smoothed-lrr K` adds a moving-median track); a long final report (`final_report.tsv.gz`: Top/Forward/Plus
|
||
alleles, GC and GT scores, LRR, BAF, position) or a genotype matrix. File schemas are OUT-xxx items of the spec.
|
||
|
||
## Accuracy
|
||
|
||
Validated phase by phase (reports in `docs/validation/`):
|
||
|
||
| Check | Result |
|
||
|---|---|
|
||
| Normalization and calls vs 2141 vendor GTCs (1000 Genomes Omni2.5, `--preset autocall-1.6.2`) | bit-identical transforms; GTCs byte-identical in 2138/2141 |
|
||
| GSA v4 demo vs gtc2vcf / GTCtoVCF / plink2 | 0 unexplained differences in GTC fields, VCF records and PLINK files |
|
||
| 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
|
||
|
||
* 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
|
||
`docs/validation/P9.md` for server numbers and thread scaling. Loading a manifest/cluster file is a fixed cost
|
||
per run (GSA ~4 s, Omni2.5 ~18 s on the Pi with the kernels).
|
||
* Methylation arrays: control-probe QC, detection and β values only (no differential analysis).
|
||
|
||
## Troubleshooting
|
||
|
||
Errors and warnings carry a stable code (`E_…`, `W_…`) — all are listed in `docs/messages.md` with their spec item.
|
||
Add `--debug` for a traceback or `--log-json` for machine-readable logs. Exit codes (CLI-014): 0 success,
|
||
1 internal error, 2 usage or pre-flight validation error (nothing processed), 3 some samples failed (or
|
||
`compare --fail-on-discrepancy` found discrepancies), 4 no sample succeeded, 130/143 interrupted, 141 the reader of
|
||
stdout went away (e.g. `simurg inspect … | head`; nothing is reported). Findings about vendor behaviour that SimurgArray reproduces
|
||
(or deliberately does not) are in `docs/known_unknowns.md`.
|