From 7b13720a4f374711a86b3e4b31a7109887984434 Mon Sep 17 00:00:00 2001 From: Kemal Yaylali Date: Wed, 30 Sep 2026 11:48:17 +0100 Subject: [PATCH] README: badges, full install options, input files, quick start, usage for every command, outputs, performance, development --- README.md | 387 ++++++++++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 349 insertions(+), 38 deletions(-) diff --git a/README.md b/README.md index 64d1a4b..9370f05 100644 --- a/README.md +++ b/README.md @@ -1,55 +1,366 @@ -# SimurgArray -

SimurgArray: Thirty birds, one picture.

-Personal research toolkit for Illumina Infinium genotyping BeadChip data -(IDAT → normalization → genotype calling → BAF/LRR → QC → VCF/PLINK; CN/LOH/mosaicism, PGx, cluster training, -methylation control QC). +

+ 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 +

-- **SIMURGARRAY_MASTER_PLAN.md** — architecture, algorithms, roadmap -- **SIMURGARRAY_FUNCTIONAL_SPEC.md** — 416 itemised requirements (CLI, inputs, outputs, QC, CN, PGx, methylation, edge cases) -- `docs/validation/` — per-phase validation reports · `docs/known_unknowns.md` — findings log -- `branding/` — logo, icons, templates and the brand guide · `blog/` — five posts about the project +

+ 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 +

-## Status +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. -**v2.0.0** (2026-09-30) — see `CHANGELOG.md`. v1: IDAT → normalization → GenCall-equivalent genotype calling → -GTC v5 → sample/control QC with an HTML report → plus-strand VCF/BCF → PLINK 1/2, bedGraph, final report; plus a -differential-testing harness. v2 adds cluster training (ArrayTrain), genome-wide CN/LOH/mosaic calling, -pharmacogenomics (targeted CN, star alleles, phenotypes) and methylation control QC. Validation (`docs/validation/`): +> **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, plink2 | zero unexplained differences | -| 1000 Genomes pedigrees (Omni2.5) | 717/717 parent–offspring pairs; Mendel error median 0.031 % | +| 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 EGT, GSE165845 chr22 vs 1000G 30x | 99.44 % concordance (vendor EGT 99.46 %), call rate 99.61 % | -| CN, 360 GSE165845 arrays vs NYGC 30x SV | PASS precision DEL 44 % / DUP 57 % | -| PGx, GSA v4 demo vs GeT-RM | 154/168 gene calls | +| 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 | -Commands: `simurg inspect`, `ingest`, `normalize`, `genotype`, `qc controls|samples|contamination|relatedness|loci|report|suggest-thresholds`, -`gtc-to-vcf`, `export plink|bedgraph|finalreport|locus-summary`, `compare`, `train`, `cn build-model|call|annotate`, -`pgx db build|coverage|copy-number|star-allele|phenotype`, `methyl qc|convert-manifest`. **User guide: -`docs/user_guide.md`**; worked example: `examples/gsa_demo.sh`. Optional Rust kernels (byte-identical outputs, -~80 GSA samples/min per server core, `docs/validation/P9.md`); wheels, Docker image and conda recipes -(`docs/validation/P10.md`). +## 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 -uv venv && uv pip install -e ".[dev]" -uv run pytest tests/unit # no downloads needed -bash scripts/bootstrap_refs.sh # reference open-source repos → refs/ (oracles) -bash scripts/fetch_data.sh p0 # ~1.4 GB of public test data → data/ -bash scripts/make_oracle_gtcs.sh # builds gtc2vcf, writes oracle GTCs -bash scripts/make_oracle_qc.sh # locus-summary oracle (BeadArrayFiles, separate venv) -bash scripts/fetch_data.sh p6 # GRCh38 FASTA (~0.9 GB download, 3.1 GB unpacked) -bash scripts/make_oracle_vcf.sh # gtc2vcf oracle VCF for the P6 tests -uv run pytest tests/oracle # real-data acceptance tests (add -m "not slow" to skip the long ones) -uv run simurg inspect data/p0/gsa_v4_pgx/idats_Demo_36/idats_Demo_36/207895980019/207895980019_R03C01_Grn.idat -bash examples/gsa_demo.sh results/gsa_demo 6 # end-to-end run on 6 demo samples +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 ``` -Vendor data files are downloaded for personal research use only and are never committed. See `THIRD_PARTY_NOTICES.md`. +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.