SimurgArray: Thirty birds, one picture.

CI Version 2.0.0 Python 3.12 and 3.13 MIT licence Linux x86_64 and aarch64 Optional Rust kernels Docker ready Runs on a Raspberry Pi 5 Research use only

Open analysis of Illumina Infinium genotyping arrays, from IDAT to genotypes, QC, VCF, copy number and pharmacogenomics.
Install · Quick start · Usage · User guide · Validation · Changelog

SimurgArray reads the raw scanner files of Infinium BeadChips and takes them all the way to analysis ready outputs. Its genotype calling reproduces the vendor's GTC files byte for byte on public data, every algorithm names its source, and every place where it differs from a reference tool is written down in `docs/known_unknowns.md`. It runs on anything from a Raspberry Pi to a cloud server. > **Research use only.** SimurgArray is not a medical device and must not be used for diagnosis or treatment > decisions, including its pharmacogenomics output. ## Contents * [Features](#features) * [Validation](#validation) * [Requirements](#requirements) * [Install](#install) * [Input files](#input-files) * [Quick start](#quick-start) * [Usage](#usage) * [Outputs at a glance](#outputs-at-a-glance) * [Performance and memory](#performance-and-memory) * [Development](#development) * [Repository layout](#repository-layout) * [Documentation](#documentation) * [Licence, data and trademarks](#licence-data-and-trademarks) ## Features | | Area | What you get | |---|---|---| | 🧬 | **Genotyping** | IDAT → normalization → GenCall equivalent calls, scores, BAF and LRR → GTC v3/v5, bit exact with the vendor software (`--preset current` or `autocall-1.6.2`) | | 📊 | **Quality control** | 21 control probe metrics, sample QC, BAF regression contamination, KING kinship and Mendel errors, sex checks, per locus summaries, and a self contained interactive HTML report | | 📄 | **Outputs** | Plus strand VCF/BCF with tabix/CSI indexes, PLINK 1 and 2, bedGraph tracks, final reports and genotype matrices, all written natively | | 🪶 | **ArrayTrain** | Train your own EGT cluster file from a cohort's GTCs instead of borrowing one | | 🧩 | **Copy number** | Genome wide CNV, LOH and mosaicism calling with an HMM, ISCN 2020 annotation, cytobands and genes | | 💊 | **Pharmacogenomics** | Targeted copy number (CYP2D6 and others), star allele diplotypes for 70 genes, CPIC and DPWG phenotypes, from open databases | | 🧪 | **Methylation QC** | EPIC and MSA control metrics, pOOBAH detection, noob, dye bias correction, β and M values | | 🔍 | **Differential testing** | `simurg compare` checks two runs field by field (EXACT / NEAR / DISCREPANCY / UNKNOWN) with an HTML dashboard | | ⚡ | **Speed** | Optional Rust kernels with identical outputs: about 80 GSA samples a minute per server core, about one sample a second on a Pi core | ## Validation Each phase has its own report in `docs/validation/`. | Check | Result | |---|---| | 2141 vendor GTCs, 1000 Genomes Omni2.5 (`--preset autocall-1.6.2`) | normalization bit identical; GTCs byte identical in 2138/2141 | | GSA demo vs gtc2vcf, GTCtoVCF and plink2 | zero unexplained differences | | 1000 Genomes pedigrees (Omni2.5) | 717/717 parent and offspring pairs; median Mendel error 0.031 % | | 360 GSA arrays (GSE165845) vs NYGC 30x WGS / GIAB HG001 | concordance 0.9934 / 0.9985; replicates 0.9966 | | ArrayTrain cluster file, GSE165845 chr22 vs 1000G 30x | 99.44 % concordance (vendor file 99.46 %) at a 99.61 % call rate | | Copy number, 360 GSE165845 arrays vs NYGC 30x SV | PASS precision DEL 44 % / DUP 57 %; rare deletions found 30 % | | Pharmacogenomics, GSA v4 demo vs GeT RM | 154/168 gene calls | | Methylation QC vs sesame | β r = 0.99998 | ## Requirements * **Python 3.12 or 3.13** on 64 bit Linux (x86_64 or aarch64). Tested on Ubuntu 24.04 and Raspberry Pi OS. * **Memory:** about 2 to 3 GB per worker process for GSA or Omni sized arrays. 8 GB is enough for a Pi with two workers. * **Optional:** Rust and a C compiler for the compiled kernels; Docker or conda for those installs. * Only pure Python dependencies are needed (numpy, pyarrow, typer, rich, pyyaml, openpyxl). bcftools and plink2 are **not** required; SimurgArray writes VCF, BCF, indexes and PLINK files itself. They are used only as test oracles. ## Install ### From source (recommended) With [uv](https://docs.astral.sh/uv/): ```bash git clone https://gitea.yaylali.uk/kemal/SimurgArray.git cd SimurgArray uv venv && uv pip install -e . # add ".[dev]" for tests, ruff and mypy source .venv/bin/activate simurg --version ``` With plain pip: ```bash python3 -m venv .venv && .venv/bin/pip install -e . ``` ### Optional compiled kernels Three to four times faster genotyping and four times 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 gcc or clang) uv pip install ./rust/simurg-core python -c "from simurg._core import backend; print(backend())" # prints: rust 1.1.0 ``` Without the kernels the same code runs in pure Python and NumPy. `SIMURG_PURE_PYTHON=1` forces that path. ### From release files `bash scripts/build_dist.sh` writes a pure wheel, an sdist and `manylinux_2_28` kernel wheels for x86_64 and aarch64 (glibc 2.28 or newer) to `dist/`: ```bash 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 kernel wheel is optional ``` ### Docker The image is about 550 MB, includes the kernels and runs as a non root user. Your 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 ```bash conda build -c conda-forge --output-folder conda-bld packaging/conda/simurg-core packaging/conda/simurg conda create -n simurg -c ./conda-bld -c conda-forge simurg simurg-core ``` All installs give identical results on x86_64 and aarch64 (`docs/validation/P10.md`). ## Input files | File | Needed by | Notes | |---|---|---| | **IDAT pairs** (`__Grn.idat` and `_Red.idat`, optionally `.gz`) | `genotype`, `qc controls`, `methyl qc` | searched recursively; the sample ID is `_` | | **BPM manifest** | almost every command | loci, bead addresses, normalization bins, positions; its genome build sets the VCF build | | **EGT cluster file** | `genotype`, `qc loci` | from the vendor, or trained with `simurg train` | | **CSV manifest** | `gtc-to-vcf` (indels), `export finalreport`, PGx | same product and version as the BPM | | **Reference FASTA** | `gtc-to-vcf`, `cn`, `pgx` | GRCh38 or GRCh37, plain or bgzip; `bash scripts/fetch_data.sh p6` downloads GRCh38 | | **Sample sheet** (optional) | most commands | `Sample_ID`, `Sample_Name`, `Sample_Group`, `Parent1`/`Parent2`, sex and custom columns | Manifests and cluster files come from the array vendor's support pages for your product. They are licensed by the vendor and are **never** committed to this repository; `.gitignore` blocks IDAT, BPM, EGT and GTC files. ## Quick start A complete run on the public GSA v4 PGx demo (36 samples): ```bash bash scripts/fetch_data.sh p0 # demo IDATs, manifest and cluster file (~1.4 GB) bash scripts/fetch_data.sh p6 # GRCh38 FASTA for VCF output (~0.9 GB download) bash examples/gsa_demo.sh results/gsa_demo 6 # end to end on 6 samples; drop the 6 for all 36 ``` The demo writes GTCs, a sample summary, an HTML QC report, per sample VCFs and a cohort BCF, PLINK files, bedGraph tracks and a final report to `results/gsa_demo/`. On a Raspberry Pi 5 all 36 samples take about ten minutes. ## Usage Every command prints its options with `simurg COMMAND --help`. Global options on every command: `--quiet`, `--log-level`, `--log-json` (JSON Lines on stderr) and `--debug`. ### Look inside files ```bash simurg inspect idats/204112350001_R01C01_Grn.idat GSA.bpm GSA.egt simurg inspect run/gtc/samples/S1/S1.gtc --json ``` ### Genotype: IDAT → GTC ```bash simurg genotype --bpm GSA.bpm --egt GSA.egt --csv GSA.csv \ --idats idats/ --sheet SampleSheet.csv --out run/gtc --threads 0 ``` * `--threads 0` uses every core. On an 8 GB machine keep it at 2 or 3. * `--preset current` (default) behaves like the current vendor software; `--preset autocall-1.6.2` reproduces the legacy AutoCall GTCs. `--cutoff` sets the GenCall no call threshold (0.15). * Writes one GTC per sample, `sample_summary.tsv`, control probe QC (`controls_qc.tsv`), relatedness, and `run_manifest.json` with every version, option and input checksum. ### Quality control ```bash simurg qc report --data run/gtc --out run/qc_report --output-format both # HTML report + xlsx/csv table simurg qc samples --bpm GSA.bpm --gtcs run/gtc --out run/qc simurg qc contamination --bpm GSA.bpm --gtcs run/gtc --out run/qc simurg qc relatedness --bpm GSA.bpm --gtcs run/gtc --sheet SampleSheet.csv --out run/qc simurg qc loci --bpm GSA.bpm --egt GSA.egt --gtcs run/gtc --out run/qc simurg qc controls --array-type genotyping --idats idats/ --out run/controls # no genotyping needed simurg qc suggest-thresholds --data run/gtc --out thresholds.yaml ``` The HTML report is a single file you can share: linked views, a thresholds editor, chip and plate heatmaps, trend plots and a PCA card. Use `--no-metadata` or `--drop-metadata COL` before sharing reports with sample sheet columns. ### VCF and exports ```bash simurg gtc-to-vcf --bpm GSA.bpm --csv GSA.csv --fasta GRCh38.fa --gtcs run/gtc --out run/vcf # per sample simurg gtc-to-vcf --bpm GSA.bpm --csv GSA.csv --fasta GRCh38.fa --gtcs run/gtc --out run/bcf --merge # one cohort BCF simurg export plink --vcfs run/vcf --gtcs run/gtc --out run/plink --format both # .bed/.bim/.fam + .pgen/.pvar/.psam simurg export bedgraph --bpm GSA.bpm --gtcs run/gtc --out run/tracks --chrom-style ucsc simurg export finalreport --bpm GSA.bpm --csv GSA.csv --egt GSA.egt --gtcs run/gtc --out run/report # --matrix for loci × samples ``` Useful `gtc-to-vcf` options: `--info-fields`, `--extra-format`, `--merge-snv-indel`, `--auxiliary-loci`, `--filter-gs`, `--drop-ambiguous`, `--ploidy-policy sex-aware|diploid`, `--include-loci` and `--exclude-loci`. ### Train your own cluster file (ArrayTrain) ```bash simurg train --gtcs cohort_gtc --bpm GSA.bpm --out mycohort.egt # at least 100 GTCs simurg genotype --bpm GSA.bpm --egt mycohort.egt --idats new_idats/ --out run2/gtc ``` Training reads raw intensities, not genotypes, so it does not matter which cluster file the cohort was first called with. Samples are chosen by sex for X and Y loci. See `docs/validation/V2_arraytrain.md`. ### Copy number, LOH and mosaicism ```bash 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/gtc --fasta GRCh38.fa --out run/cn simurg cn annotate --vcfs run/cn ``` Build the model once per manifest and batch from at least 24 normal samples. Each sample gets a CNV VCF, a segments file and, after `annotate`, ISCN strings with cytobands and genes. ### Pharmacogenomics (GRCh38) ```bash bash scripts/pgx/fetch_sources.sh # CPIC, PharmVar, ClinPGx → data/pgx/sources simurg pgx db build --sources data/pgx/sources --fasta GRCh38.fa --manifests GSA_PGx.bpm --csv GSA_PGx.csv --out pgxdb/ 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/gtc --model cnmodel/GSA_PGx.simurg-pgxcn --fasta GRCh38.fa --out run/vcf simurg gtc-to-vcf --gtcs run/gtc --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 # or dpwg ``` `pgx db build` names the database after its content hash (`pgx_db-1-.sqlite`); use that file name. Copy number training needs at least 96 samples. Outputs per sample: diplotypes, supporting and missing variants, alternative diplotypes the array cannot tell apart, and a JSON report with phenotype and evidence. ### Methylation arrays ```bash simurg methyl qc --manifest EPIC_manifest.csv --idats idats/ --sheet SampleSheet.csv --out run/methyl --profile auto simurg methyl convert-manifest --manifest EPIC_manifest.csv --out manifest_tables/ ``` ### Compare two runs ```bash simurg compare runA/gtc runB/gtc --out compare_report --fail-on-discrepancy ``` Every field of every GTC or VCF pair is compared. `compare_report/dashboard.html` shows the result and `discrepancies.tsv.gz` lists every difference. ### Exit codes `0` success · `1` internal error · `2` usage or validation error (nothing processed) · `3` some samples failed, or `compare --fail-on-discrepancy` found differences · `4` no sample succeeded · `130`/`143` interrupted. Every error and warning has a stable code (`E_…`, `W_…`) listed in `docs/messages.md`. ## Outputs at a glance | Step | Main files | |---|---| | `genotype` | `samples//.gtc`, `sample_summary.tsv`, `controls_qc.tsv`, `relatedness.tsv`, `run_manifest.json` | | `qc samples` | `sample_qc.tsv`, `contamination.tsv`, `relatedness.tsv` (+ `mendel.tsv` with parents in the sheet) | | `qc report` | `qc_report_.html`, `qc_table_.xlsx` | | `gtc-to-vcf` | `.snv.vcf.gz` + `.tbi`, or `cohort.snv.bcf` + `.csi` | | `export plink` | `.bed/.bim/.fam`, `.pgen/.pvar/.psam` | | `cn call` | `.cyto.cnv.vcf.gz`, `.cn_segments.parquet`, `cn_sample_summary.tsv` | | `pgx phenotype` | `.pgx.json`, `diplotype_summary.tsv`, `diplotype_annotations.tsv` | | `train` | `.egt`, `.egt.json` | Every HTML report carries the software version and a checksum of the data it shows. ## Performance and memory * With the kernels a GSA sample takes about 1.1 s on one Raspberry Pi 5 core (3.4 s pure Python) and a server core handles about 80 samples a minute (`docs/validation/P9.md`). * Loading a manifest is a fixed cost per command: about 5 s for GSA and 15 s for Omni2.5 on a Pi. * Cohort QC steps grow with samples × loci. On small machines use `--threads 2`, and `--no-cohort-qc` on very large `genotype` runs, then run the `qc` commands separately. * `simurg train` needs about 2 GB for 120 GSA samples and takes about 70 minutes on a Pi. ## Development ```bash uv venv && uv pip install -e ".[dev]" uv run ruff check . && uv run ruff format --check . uv run mypy uv run pytest tests/unit # no downloads needed ``` Oracle tests compare against reference tools on real public data: ```bash bash scripts/bootstrap_refs.sh # reference open source repositories → refs/ bash scripts/fetch_data.sh p0 # public test data → data/ bash scripts/make_oracle_gtcs.sh # builds gtc2vcf and writes oracle GTCs bash scripts/make_oracle_qc.sh # locus summary oracle (BeadArrayFiles) bash scripts/fetch_data.sh p6 && bash scripts/make_oracle_vcf.sh uv run pytest tests/oracle # add -m "not slow" to skip the long ones ``` CI runs lint, type checks, a secret scan and the unit tests with and without the Rust kernels on Python 3.12 and 3.13 for every push. `scripts/nightly.py` reruns the differential tests against the oracles. ## Repository layout ``` src/simurg/ the package: formats, signal, genotype, qc, harmonize, cn, pgx, train, pipeline, cli rust/simurg-core/ optional compiled kernels (PyO3) tests/ unit and oracle tests scripts/ data fetchers, oracle builders, validation and cloud helpers, brand and blog image scripts docs/ user guide, messages, validation reports, known unknowns examples/ end to end demo branding/ logo, icons, templates and the brand guide blog/ five posts about the project, with featured images packaging/ conda recipes (plus Dockerfile at the root) ``` ## Documentation * [`docs/user_guide.md`](docs/user_guide.md): every command in detail, presets, QC, outputs, accuracy, limits, troubleshooting * [`SIMURGARRAY_FUNCTIONAL_SPEC.md`](SIMURGARRAY_FUNCTIONAL_SPEC.md): 416 numbered requirements (CLI, inputs, outputs, QC, CN, PGx, methylation, edge cases) * [`SIMURGARRAY_MASTER_PLAN.md`](SIMURGARRAY_MASTER_PLAN.md): architecture, algorithms with their sources, roadmap * [`docs/validation/`](docs/validation/): one report per phase and module * [`docs/known_unknowns.md`](docs/known_unknowns.md): every difference from reference tools, with dates and explanations * [`docs/messages.md`](docs/messages.md): all error and warning codes * [`branding/README.md`](branding/README.md) and [`blog/`](blog/) ## Licence, data and trademarks SimurgArray is released under the [MIT licence](LICENSE). Code adapted from MIT, BSD and Apache licensed projects is credited in [`THIRD_PARTY_NOTICES.md`](THIRD_PARTY_NOTICES.md); GPL tools are only run as separate programs for comparison. Brand artwork is CC BY 4.0. Vendor manifests, cluster files and scan data are downloaded by each user for their own research and are never part of this repository. PGx sources keep their own terms (CPIC CC0, PharmVar terms of use, ClinPGx CC BY SA 4.0); SimurgArray ships the database builder, not the database. Illumina, Infinium, GenCall, GenTrain and GenomeStudio are trademarks of Illumina, Inc. SimurgArray is an independent project, not affiliated with or endorsed by Illumina.