SimurgArray/docs/user_guide.md

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 · Inputs · A complete run · Commands · Methylation · CN/LOH · PGx · Cluster training · 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.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):

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

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

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

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

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

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 ~5 s, Omni2.5 ~15 s on the Pi).
• 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.