SimurgArray/docs/user_guide.md

SimurgArray user guide (v1.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 · Inputs · A complete run · Commands · Presets and reproducibility · QC · Outputs · Accuracy · Limitations · Troubleshooting

Install

Python ≥ 3.12. With uv:

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):

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.0.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):

python3 -m venv simurg-env && simurg-env/bin/pip install simurg-1.0.0-py3-none-any.whl \
    simurg_core-1.0.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):

docker build -t simurg:1.0.0 .
docker run --rm --user "$(id -u):$(id -g)" -v "$PWD":/work simurg:1.0.0 \
    genotype --bpm M.bpm --egt C.egt --idats idats --out run/gtc

conda (recipes in packaging/conda/):

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 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).

# 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 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.

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

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.
• 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 ~5 s, Omni2.5 ~15 s on the Pi).
• Methylation arrays: control-probe QC only.

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.