# SimurgArray — Master Plan & Technical Specification > **SimurgArray** is an open, reproducible, personal-research toolkit for analysing Illumina Infinium > genotyping BeadChip data: IDAT → normalization → genotype calling → BAF/LRR → QC → VCF/PLINK, with later > modules for copy number, LOH/mosaicism, pharmacogenomics and cluster-file training. > > This document is **self-contained**: a developer on a fresh machine should be able to > build the whole project from it without any other notes. Version 1.0 · 2026-09-27. > > *Name:* the Simurg (Anka / Zümrüdüanka) is the mythical bird of Turkish–Persian legend. CLI and Python package: `simurg`. --- ## Table of contents 0. How to use this document 1. Scope, goals, non-goals 2. Ground rules (licensing & provenance) 3. Target platform & environment setup 4. Architecture & repository layout 5. Canonical data model 6. File-format specifications (IDAT, BPM, CSV manifest, EGT, GTC, sample sheet) 7. Core algorithms (exact math) 8. Strand, alleles & VCF export 9. QC module 10. Copy number, LOH & mosaicism module 11. Pharmacogenomics module 12. Cluster-file training module (ArrayTrain) 13. SimurgArray output specification 14. CLI specification 15. Provenance & reproducibility 16. Testing & validation strategy 17. Datasets (with URLs) and fetch plan 18. Phased roadmap with acceptance criteria 19. Performance plan 20. Known unknowns & planned experiments 21. Risk register 22. References - Appendix A — numeric helper definitions (percentile, trimmean, robustfit, fuzzy memberships) --- ## 0. How to use this document - **Detailed functional specification:** `SIMURGARRAY_FUNCTIONAL_SPEC.md` is the authoritative spec for externally visible behaviour: CLI options, validation, file schemas, QC semantics, messages, edge cases and acceptance tests. Each of its items has an ID and is traced to a § and phase of this plan. This plan remains authoritative for algorithms and numerics. - Work **phase by phase** (§18). Do not start a phase before the previous phase's acceptance criteria pass. - Every algorithm in §7 is written as exact math with its **evidence source** in brackets: **[DOC]** Illumina public documentation / tech notes · **[LIT]** papers & expired patents · **[OSS]** open-source implementation (file:line where useful) · **[INF]** our inference · **[HYP]** hypothesis to test. - Where the spec says "bit-exact mode" it means: reproduce Illumina's published GTC outputs (see §16 oracles) including their numeric quirks. "Clean mode" means the mathematically tidy version. Both modes are kept; the default is bit-exact mode so results are comparable with files produced by Illumina software. --- ## 1. Scope, goals, non-goals ### 1.1 Goals (v1.0) 1. Read IDAT (v3 binary), BPM, CSV manifest, EGT, GTC (v3–v5), sample sheets. 2. Assemble per-locus raw X/Y from Red/Green IDATs. 3. Fit and apply the per-bin affine intensity normalization. 4. Call genotypes from an existing cluster file (EGT) with a GenCall-equivalent score. 5. Compute θ, R, BAF, LRR, LogR deviation, call rates, GC10/GC50, intensity percentiles, sex estimate. 6. Control-probe QC metrics; sample- and locus-level QC tables; HTML QC report. 7. Write GTC v5 (interoperable with Illumina's BeadArrayFiles, Picard, gtc2vcf), VCF/BCF (with GT:GS:BAF:LRR), PLINK 1.9/2 files, bedGraph tracks, TSV/JSON summaries. 8. Allele/strand harmonization with a full audit trail. 9. Deterministic, provenance-tracked, resumable batch pipeline; scalable to ≥10 000 samples. ### 1.2 Goals (v2.0) 10. ArrayTrain: our own cluster training → EGT files (so SimurgArray is not dependent on vendor EGTs). 11. Genome-wide CNV/LOH/mosaic calling from LRR/BAF (HMM), with ISCN-like annotation. 12. Targeted PGx copy number (e.g. CYP2D6 sub-regions) and star-allele diplotype calling from array VCFs. 13. Contamination (BAF regression) and relatedness/sample-identity QC. 14. Methylation-array control QC (optional). ### 1.3 Non-goals - Encrypted (DES/XML, "version 1") IDATs — detect and refuse with a clear error. - Clinical/diagnostic use. SimurgArray is a research tool. - Any use of vendor binaries' internals (see §2). --- ## 2. Ground rules (licensing & provenance) SimurgArray is a **personal, non-commercial research** project. 1. **Allowed sources:** public documentation and tech notes, papers, expired patents (US7035740 and US7467117, both expired 2024), public file-format specs, open-source code under its licence, public datasets, and our own experiments. 2. **Open-source reuse is permitted and encouraged**, with attribution (keep a `THIRD_PARTY_NOTICES.md`): - **MIT / BSD / Apache-2.0 / Unlicense** code may be ported or vendored: gtc2vcf/idat2gtc/BAFregress (MIT), mocha (MIT), Illumina BeadArrayFiles (BSD-2), Illumina GTCtoVCF (Apache-2.0), Picard `picard.arrays` (MIT), PennCNV (Unlicense), PyPGx (MIT), bcftools (MIT). - **GPL** tools (plink2, illuminaio, verifyIDintensity, GenotypeHarmonizer, akt) are used only as **external subprocess oracles**; never copied or linked into SimurgArray so SimurgArray can stay MIT-licensed. - **Do not use** code without a licence or with non-commercial/PolyForm terms (zCall, argyle, Aldy, Cyrius, Stargazer, QuantiSNP binaries) — ideas from their *papers* are fine. 3. **Never** decompile, disassemble or inspect vendor software binaries, and never bypass licence checks. 4. **Vendor data files** (Illumina BPM/CSV/EGT, demo IDATs) are downloaded by the user for personal research use; they are **not committed** to the repository. The repo only contains a fetch script with checksums. 5. SimurgArray documentation and code refer to *Illumina file formats* (IDAT, BPM, EGT, GTC) as interoperability facts. It does not use or imply any vendor product brand. --- ## 3. Target platform & environment setup ### 3.1 Why move off the Raspberry Pi The code itself would run on ARM, but the Pi (4 cores, 8 GB RAM, aarch64) is too small for the validation data (1000 Genomes Omni set is 107 GB compressed; cohort-level steps need tens of GB of RAM) and too slow for iteration. Some optional reference tools only ship x86-64 binaries. ### 3.2 Recommended machine | | Minimum | Recommended | |---|---|---| | CPU | x86-64, 8 cores | 16–32 cores (AVX2/AVX-512) | | RAM | 32 GB | 64–128 GB | | Disk | 500 GB SSD | 1–2 TB NVMe (+ HDD for raw archives) | | OS | Ubuntu 22.04/24.04 LTS | Ubuntu 24.04 LTS | | GPU | not needed | optional (only for ArrayTrain at ≥10k samples) | A cloud VM works (e.g. 16 vCPU / 64 GB / 1 TB SSD). macOS (Apple Silicon) also works for development but use Linux for large validation runs. ### 3.3 Toolchain ```bash # system sudo apt update && sudo apt install -y build-essential git curl wget pkg-config cmake \ zlib1g-dev libbz2-dev liblzma-dev libcurl4-openssl-dev libssl-dev libdeflate-dev \ tabix samtools bcftools aria2 pigz jq # Python via uv (fast, reproducible) curl -LsSf https://astral.sh/uv/install.sh | sh uv python install 3.12 # Rust (Phase 9) curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh -s -- -y # Reference/oracle tools (subprocess only) # gtc2vcf plugins for bcftools: build from https://github.com/freeseek/gtc2vcf (instructions in its README) # plink2: https://www.cog-genomics.org/plink/2.0/ # R + Bioconductor illuminaio (optional oracle): install.packages("BiocManager"); BiocManager::install("illuminaio") # Java 17 + Picard jar (optional oracle) ``` ### 3.4 Python dependencies (justify each) | Package | Why | |---|---| | numpy | all numerics (float32/float64 control) | | scipy | stats, special functions, optimisation, HMM helpers | | pyarrow | canonical columnar model, Parquet/Arrow IPC, memory-mapping | | polars | fast tabular QC/summary operations | | pysam or cyvcf2 | VCF/BCF writing/reading via htslib | | typer + rich | CLI and progress reporting | | pydantic | config/schema validation | | pyyaml | pipeline config | | structlog | JSON structured logging | | hypothesis, pytest, pytest-benchmark | tests | | scikit-learn (v2 only) | GMM baselines for ArrayTrain experiments | | maturin + PyO3 (Phase 9) | Rust extension | ### 3.5 Moving the project The canonical repository is the Gitea repo **https://gitea.yaylali.uk/kemal/SimurgArray**. ```bash # on the new machine git clone https://gitea.yaylali.uk/kemal/SimurgArray.git && cd SimurgArray bash scripts/bootstrap_refs.sh && bash scripts/fetch_data.sh p0 ``` --- ## 4. Architecture & repository layout ### 4.1 Principles - Scientific algorithms are **pure functions** over NumPy/Arrow arrays: no file I/O, no logging side-effects, no global state. Each returns `(result, AuditRecord)`. - Layers: `formats` (I/O) → `core` (data model) → `signal` / `genotype` / `cn` / `qc` / `pgx` (science) → `pipeline` (orchestration) → `cli` / `api` (interfaces). - Per-sample work is embarrassingly parallel; only cluster training, PGx-CN batch models, locus summaries and relatedness are cohort-wide. - Every numeric step has a declared dtype policy (float32 where the GTC stores float32; float64 internally when "clean mode"). ### 4.2 Layout ``` simurg/ ├── pyproject.toml # uv/hatch; package "simurg"; entry point `simurg` ├── SIMURGARRAY_MASTER_PLAN.md # this file ├── THIRD_PARTY_NOTICES.md ├── scripts/ │ ├── bootstrap_refs.sh # clones reference repos into refs/ (git-ignored) │ └── fetch_data.sh # downloads datasets into data/ with sha256 checks (git-ignored) ├── src/simurg/ │ ├── core/ # schemas.py (Arrow), types.py, audit.py, provenance.py, errors.py │ ├── formats/ │ │ ├── binutil.py # little-endian readers, .NET 7-bit varint strings │ │ ├── idat.py bpm.py csv_manifest.py egt.py gtc.py samplesheet.py │ │ ├── vcf.py plink.py bedgraph.py finalreport.py │ ├── signal/ │ │ ├── raw_xy.py # IDAT→locus X/Y assembly (§7.1) │ │ ├── normalization/ # fit.py (§7.2), apply.py (§7.3), numerics.py (Appendix A) │ │ └── polar.py # θ, R (§7.4) │ ├── genotype/ │ │ ├── gencall.py # EGT-based caller + score (§7.5) │ │ ├── baf_lrr.py # §7.6 │ │ ├── sex.py # §7.8 │ │ └── train/ # ArrayTrain (§12) │ ├── qc/ # controls.py, sample.py, locus.py, contamination.py, relatedness.py, report/ │ ├── harmonize/ # strand.py, alleles.py, indels.py (§8) │ ├── cn/ # gc_wave.py, hmm.py, segment.py, mosaic.py, annotate.py (§10) │ ├── pgx/ # db.py, cn_model.py, star_allele.py, phenotype.py (§11) │ ├── pipeline/ # dag.py, runner.py, cache.py, config.py │ ├── cli/ # main.py (typer) │ └── synth/ # synthetic data generator (§16.4) ├── rust/simurg-core/ # Phase 9: idat parse, normalization, calling (PyO3) ├── tests/ │ ├── unit/ property/ integration/ regression/ oracle/ perf/ │ └── fixtures/ # ONLY redistributable data (Bioconductor Artistic-2.0 files, synthetic) ├── docs/ # user guide, algorithm notes, validation reports ├── refs/ # (git-ignored) cloned third-party repos └── data/ # (git-ignored) downloaded datasets ``` --- ## 5. Canonical data model Arrow schemas (Parquet on disk). Names are SimurgArray's own. **`loci`** (one row per manifest locus, BPM order) `locus_idx:int32, name:str, ilmn_id:str, chrom:str, pos:int64, snp:str ("[A/G]"), ilmn_strand:str (TOP/BOT/PLUS/MINUS), source_strand:str, ref_strand:str(+/-/null), address_a:int32, address_b:int32, assay_type:uint8 (0 InfII, 1 InfI-AT red, 2 InfI-CG green), norm_id:uint8, norm_bin:uint16 (after §6.2 rule), exp_clusters:int8, intensity_only:bool, ploidy:str, genome_build:str, source_seq:str, top_genomic_seq:str, is_control:bool` **`clusters`** (EGT, joined to loci by name) `locus_idx, n_aa, n_ab, n_bb:int32, theta_mean_aa/ab/bb, theta_dev_aa/ab/bb, r_mean_aa/ab/bb, r_dev_aa/ab/bb:float32, intensity_threshold:float32, cluster_sep:float32, gentrain_score:float32, orig_score:float32, edited:bool` **`sample_intensity`** (one per sample) `raw_x, raw_y:uint16[n_loci]` plus control arrays; metadata: barcode, position, sample_id, scanner info, idat hashes. **`sample_calls`** (one per sample) `genotype:uint8 (0 NC,1 AA,2 AB,3 BB), score:float32, norm_x, norm_y, theta, r, baf, lrr:float32` + `transforms` table (one row per norm_bin: version, offset_x, offset_y, scale_x, scale_y, shear, theta, nn12, rr12, taa, tbb). **`sample_summary`**: call_rate, autosomal_call_rate, logr_dev, gc10, gc50, sex_est, n_calls, n_nocalls, n_intensity_only, p05/p50/p95 x & y, contamination, qc flags. **`audit`**: every transformation record `{step, locus_idx|null, sample, input, output, rule, reason, evidence_tag}`. --- ## 6. File-format specifications All binary formats are **little-endian**. Strings in BPM/EGT/GTC are **.NET BinaryWriter strings**: a 7-bit varint (LEB128) length prefix followed by UTF-8 bytes (the GTC PDF implies a single length byte — only true for < 128 bytes). [OSS: BeadArrayFiles `BeadArrayUtility.read_string`] ### 6.1 IDAT (binary, version 3) [OSS: illuminaio, BeadArrayFiles, gtc2vcf] ``` char[4] magic = "IDAT" int64 version # must be 3; version 1 = encrypted XML → refuse int32 n_fields n_fields × { uint16 code; int64 offset } # table of contents ``` | code | meaning | type at offset | |---|---|---| | 1000 | number of addresses N | int32 | | 102 | address IDs ("IlluminaID") | int32[N] | | 103 | SD | uint16[N] | | 104 | Mean (scanner trimmed mean) | uint16[N] | | 107 | NBeads | uint8[N] | | 200 | MidBlock | int32 count + array (**element type disputed: int32 vs uint8 — check against files**) | | 300 | RunInfo | int32 n, then n × 5 strings (time, block type, pars, code, version) | | 400 | RedGreen | 4 bytes | | 401 | Manifest | string | | 402 | Barcode | string | | 403 | ChipType | string | | 404 | SentrixPosition (R01C01) | string | | 405 | Beadset | string | | 406 | SampleID | string | | 407 | Description | string | | 408 | Plate | string | | 409 | Well | string | | 410 | SampleCount | int32 | | 510 | unknown | string | - One file per colour: `__Grn.idat`, `..._Red.idat` (optionally `.gz`). - **Join Grn and Red by address ID, not by position** (build address→index hash). - Validation: magic, version, TOC offsets within file, N consistent across arrays, both colours share barcode/position. - Implementation: memory-map; decode arrays with `numpy.frombuffer` at TOC offsets (zero-copy). ### 6.2 BPM manifest (binary) [OSS: BeadArrayFiles `BeadPoolManifest.py`; gtc2vcf] ``` char[3] "BPM"; uint8 1 int32 version # if version & 0x1000: version ^= 0x1000 ; supported 3..5 string manifest_name string control_config # (version > 1) CSV-ish lines "addr[:addr..],Category,Color,Name" int32 num_loci int32[num_loci] indexes # 1-based string[num_loci] names uint8[num_loci] norm_ids # < 100 num_loci × LocusEntry [int32 m; m × string header lines] # trailer (optional; gtc2vcf reads it) ``` **LocusEntry** (int32 entry version ∈ {4,6,7,8}): strings IlmnID, Name, 3× unused, **int32 index**, unused string, IlmnStrand, SNP, Chr, Ploidy, Species, MapInfo (string), TopGenomicSeq (v4 only), CustomerStrand; int32 AddressA, int32 AddressB; strings AlleleA_ProbeSeq, AlleleB_ProbeSeq (v4 only), GenomeBuild, Source, SourceVersion, SourceStrand, SourceSeq (v4 only); **v≥6:** uint8 unused, int8 Exp_Clusters, int8 Intensity_Only, uint8 AssayType (0 InfII, 1 InfI A/T, 2 InfI G/C; rule: AddressB==0 ⇔ AssayType==0); **v≥7:** 4 float (frac A,C,T,G); **v8:** string RefStrand. For v<6 infer AssayType from the bases around `[x/y]` in SourceSeq (gtc2vcf `get_assay_type`). **Normalization bin rule (critical):** ``` norm_id' = (norm_id + 100 * assay_type) mod 256 # mimics vendor byte wrap (affects Omni5) norm_bin = rank of norm_id' among sorted distinct norm_id' values transform for locus i = GTC.transforms[norm_bin[i]] ``` ### 6.3 CSV manifest [DOC: Illumina manifest column descriptions] Sections `[Heading]`, `[Assay]`, `[Controls]`. Assay columns: IlmnID, Name, IlmnStrand, SNP, AddressA_ID, AlleleA_ProbeSeq, AddressB_ID, AlleleB_ProbeSeq, GenomeBuild, Chr, MapInfo, Ploidy, Species, Source, SourceVersion, SourceStrand, SourceSeq, TopGenomicSeq, BeadSetID, Exp_Clusters, Intensity_Only, RefStrand. Exp_Clusters: 1 non-polymorphic, 2 mtDNA/Y, 3 diploid. Intensity_Only=1: always no-call, used for CN. The CSV is required for indel REF/ALT resolution (SourceSeq/TopGenomicSeq). ### 6.4 EGT cluster file [OSS: BeadArrayFiles `ClusterFile.py`, gtc2vcf `egt_init`, Picard `InfiniumEGTFile`] ``` int32 version (=3) string gencall_version, cluster_version, call_version, normalization_version, date_created uint8 is_wgt (=1) string manifest_name int32 data_block_version # 5..9 string opa int32 num_records num_records × ClusterRecord { int32 N_AA, N_AB, N_BB float Rdev_AA, Rdev_AB, Rdev_BB float Rmean_AA, Rmean_AB, Rmean_BB float Tdev_AA, Tdev_AB, Tdev_BB float Tmean_AA, Tmean_AB, Tmean_BB dbv ≥ 7 (gtc2vcf) or == 9 (BeadArrayFiles): float intensity_threshold + 14 unused floats # check which! elif dbv ≥ 6: 15 unused floats ; elif dbv < 5: 2 unused floats } num_records × { float cluster_separation, float total_score (=GenTrain score), float original_score, uint8 edited } num_records × string genotype (unused) num_records × string locus_name # join to BPM by NAME num_records × int32 address num_records × int32[3] cluster counts (assert == N_AA, N_AB, N_BB) if dbv == 9: num_records × float (unknown) ``` Derived: `R̄ = (N_AA·Rmean_AA + N_AB·Rmean_AB + N_BB·Rmean_BB)/(N_AA+N_AB+N_BB)` = expected R for intensity-only loci. ### 6.5 GTC genotype file v3–v5 [DOC: BeadArrayFiles `docs/GTC_File_Format_v5.pdf`; OSS: `GenotypeCalls.py`] ``` char[3] "gtc"; uint8 version (3,4,5) int32 M; M × { int16 id; uint32 offset } # for ids 1,2,3 the "offset" IS the value arrays = int32 count + elements ``` | ID | Name | Type | |---|---|---| | 1 | NumSNPs | int32 in TOC | | 2 | Ploidy | int32 in TOC | | 3 | PloidyType (1 diploid, 2 autopoly, 3 allopoly) | int32 in TOC | | 10/11/12 | sample name / plate / well | string | | 100 | cluster file name | string | | 101 | manifest name | string | | 200 | imaging date | string | | 201 | autocall date | string | | 300 | autocall version | string | | 400 | normalization transforms | XForm[] (52 bytes each) | | 500/501 | raw control X / Y | uint16[] | | 1000/1001 | raw X / raw Y | uint16[] | | 1002 | genotypes | uint8[] (0 NC,1 AA,2 AB,3 BB,4 NULL,5 A,6 B, 7.. polyploid) | | 1003 | base calls (TOP strand, '-' = NC) | char[2][] | | 1004 | GenCall scores | float[] | | 1005 | scanner data (name, pmt_green int32, pmt_red int32, version, user) | struct | | 1006 | call rate | float | | 1007 | gender 'M'/'F'/'U' | char | | 1008 | LogR dev | float | | 1009 | p10 GC | float | | 1010 | DX flag | int32 | | 1011 | extended sample data: p50GC float, NumCalls, NumNoCalls, NumIntensityOnly int32 | struct | | 1012 | B allele freq | float[] (v≥4) | | 1013 | LogR ratio | float[] (v≥4) | | 1014/1015 | intensity percentiles X / Y (P05,P50,P95) | 3×uint16 (v5) | | 1016 | Sentrix ID | string (v5) | - **XForm (52 bytes):** int32 version; float offset_x, offset_y, scale_x, scale_y, shear, theta; 6 reserved floats (gtc2vcf stores cvx, cvy, nn12, rr12, taa, tbb there). - GTC does **not** store normalized X/Y, θ, R — they are recomputed (§7.3–7.4). - Write-complete checks: v3 NumIntensityOnly readable; v4 LogR readable; v5 Sentrix ID readable. - Open question: GTC PDF calls raw X "green"; open implementations use **X = red** for Infinium II. Verify (§20). ### 6.6 Sample sheet GenomeStudio-compatible CSV; optional `[Header]` / `[Manifests]` / `[Data]` sections. Required columns `SentrixBarcode_A, SentrixPosition_A`; optional `Sample_ID, Sample_Name, Sample_Plate, Sample_Well, Gender, Parent1, Parent2, Replicate`. Case-insensitive headers. Sample name precedence: sheet `Sample_Name` → IDAT SampleID → `barcode_position`. Without a sheet, discover IDAT pairs recursively in a folder. ### 6.7 GenomeStudio Final Report (read-only, for comparison) Tab-delimited, `[Header]` then `[Data]`; columns such as SNP Name, Sample ID, Allele1/2 - Top/Forward/AB/Plus, GC Score, Theta, R, X, Y, X Raw, Y Raw, B Allele Freq, Log R Ratio. Used only as an oracle (§16). --- ## 7. Core algorithms (exact math) Notation: per locus i raw intensities (x, y) as float; P_p(v) = Illumina percentile (Appendix A); trimmeanK(v) = mean of v within [P_{K/2}, P_{100−K/2}]. ### 7.1 Raw X/Y assembly [OSS: gtc2vcf `fill_array` idat2gtc.c:3711-3745; LIT: Steemers 2006] - Infinium II (assay 0): **X = Red.Mean[addrA], Y = Grn.Mean[addrA]**. - Infinium I A/T (assay 1): X = Red.Mean[addrA], Y = Red.Mean[addrB]. - Infinium I C/G (assay 2): X = Grn.Mean[addrA], Y = Grn.Mean[addrB]. - If any used address has **NBeads < 2** or is missing → (0, 0) ("missing"). - Controls: control_x = Red, control_y = Grn for each control address in BPM control config. - Chemistry reason: Infinium II extends the SNP base itself (A/T → red dye, C/G → green), so [A/T] and [C/G] SNPs need Infinium I (two probes, same colour). ### 7.2 Normalization fit (per sample, per bin) [LIT: US7035740, Peiffer 2006; DOC: GenTrain 3.0 tech note; OSS: idat2gtc.c:2352-2945] Two versions exist: **1.1.2** (older, "GenTrain 2") and **1.2.0** (current, "GenTrain 3"). Default = 1.2.0. Constants: `SAMPLING = 400`, `ROBUST_THRESHOLD = 192`. A bin with < 10 loci is an error. In bit-exact mode, missing loci (0,0) **stay in** the bin as zeros (vendor quirk). **(a) Small Infinium-I bin, v1.2.0 only** (norm_id' ≥ 100 and n < 192) — 2-parameter single-channel model: 1. Mirror each point so that x ≥ y. 2. c = P_2(y); x ← x − c, y ← y − c. 3. r = |x|+|y|, t = (2/π)·atan2(y, x). 4. AA = {t < 0.1}. ȳ = trimmean20(y_AA) − 50. 5. r̄ = trimmean20(r with NaNs removed and array padded back with zeros — quirk) − 2ȳ. 6. offset_x = offset_y = c + ȳ; θ = 0; shear = 0; scale_x = scale_y = r̄. **(b) Otherwise — six-parameter fit:** 1. **Outlier removal** (only if n ≥ 400): for each of x, y and t = y/(x+y): low = min(5th smallest, P_1), high = max(5th largest, P_99); drop points with any value ≤ low or ≥ high. 2. **Offset** (if n < 192: offset = 0): - shift x, y by their minima (mx, my); - x-sweep: 400 points linspace(0, max x) on the x-axis; for each, nearest data point (Euclidean) → unique set = AA candidates; y-sweep analogously → BB candidates. (For ≥2000 points the vendor uses a binned approximate nearest-neighbour search `findClosestSitesToPointsAlongAxis` — reproduce it in bit-exact mode; see gtc2vcf `nearest_neighbor.c`.) - robust fits (bisquare, Appendix A): y = b_a + m_a·x on AA; x = b_b + m_b·y on BB; - intersection o_x = (b_b + m_b·b_a)/(1 − m_b·m_a), o_y = b_a + m_a·o_x; - offset = (o_x + mx, o_y + my); subtract. 3. **nn12 / rr12** diagnostics: nn12 = fraction with x < 0.2·P_95(x) and y < 0.2·P_95(y); rr12 = mean(x+y of those) / mean(x+y of all). 4. **Rotation** (if n < 192: θ = 0): new x-sweep over [min x, max x], union with AA; robust fit through origin y = m·x (`robustfit0`, see quirk in Appendix A); θ = atan(m); rotate x′ = cosθ·x + sinθ·y, y′ = −sinθ·x + cosθ·y. 5. **Shear** (if n < 192: shear = 0): y-sweep over [0, max y], union with BB; fit x = m·y through origin; shear = m; x ← x − shear·y. 6. **Scale**: - n ≥ 192: AA = {x > 5y and x > 0.1·P_99(x)}, BB = {y > 5x and y > 0.1·P_99(y)}; scale_x = trimmean50(x_AA), scale_y = trimmean50(y_BB). - n < 192: clamp negatives to 0; angles in degrees; AA = {θ° < 10}, BB = {θ° > 80} (v1.1.2 also counts zero points as AA; v1.2.0 does not use missing loci); scale = trimmean20. 7. Store diagnostics: taa = θ·180/π, tbb = (atan(shear) − θ)·180/π. Differences 1.1.2 → 1.2.0: new `madsigma`; scale step ignores missing loci in bins < 192; new single-channel Infinium-I model for small bins. ### 7.3 Apply the transform [DOC/OSS: BeadArrayFiles `GenotypeCalls.py:668-698`] ``` if x_raw <= 0 and y_raw <= 0: (xn, yn) = (NaN, NaN) tx = x_raw − offset_x ; ty = y_raw − offset_y tx2 = cosθ·tx + sinθ·ty ; ty2 = −sinθ·tx + cosθ·ty tx3 = tx2 − shear·ty2 xn = max(tx3, 0) / scale_x ; yn = max(ty2, 0) / scale_y ``` Bit-exact mode: cosθ, sinθ computed in float32 (gtc2vcf does this). ### 7.4 Polar coordinates [DOC: GenTrain 2 tech note] θ = (2/π)·atan2(yn, xn) ∈ [0,1]; R = xn + yn. (0,0) → NaN. Bit-exact: R = |xn|+|yn|; both negative → R = FLT_MIN·FLT_EPSILON. ### 7.5 Genotype calling from an EGT (GenCall-equivalent) [DOC: GenCall tech note; LIT: Kermani patents; OSS: idat2gtc.c:3022-3405] Inputs per locus: cluster params (μθ_k, σθ_k, μR_k, σR_k, N_k) for k ∈ {AA, AB, BB}, GenTrain score G, intensity threshold τ; sample point (t, r). 1. Intensity-only locus → NC (counted separately). If r < τ → NC, score NaN. 2. Floor: σR_k ← max(σR_k, 0.1). 3. **Dispersion adjustment ("flair"):** z1 = ½(μθ_AB − μθ_AA)/(σθ_AB + σθ_AA); z2 = ½(μθ_BB − μθ_AB)/(σθ_BB + σθ_AB); mz = min(z1, z2); η = max(mz/3, 1)·min(mz/3, 1); σθ_k ← max(η·σθ_k, 0.02); m = max(0.2, median of original σθ); σR_k ← max(σR_k, max(μR_k/7, m)). 4. **Likelihoods:** d_k = 100·|t − μθ_k| + |r − μR_k|; choose the two nearest clusters; for those, use half-flattened Gaussian pdfs (density held at peak value on the "outward" side, depending on cluster ordering in θ and R — six cases, see idat2gtc `modilik` :3077-3159); third cluster L = 0; if no case applies use plain Gaussians for all: L_k = φ*(t; μθ_k, σθ_k)·φ*(r; μR_k, σR_k). 5. **Priors:** P_k = (N_k/N + 1)/4. **Posteriors:** AP_k = L_k·P_k / Σ_j L_j·P_j. **Call** = argmax_k AP_k. 6. **Sub-scores** (g = called cluster): - s_ap = smf(log10(AP_max / AP_2nd); 0, 2) - s_r = smf(r; 0, 0.1) · smf(r; 0, μR_g/10) · smf(r; μR_g − 6σR_g, μR_g − 2σR_g) · zmf(r; μR_g + 6σR_g, μR_g + 12σR_g) - s_θ: AA → zmf(t; μθ+2σθ, μθ+6σθ); AB → zmf(|t−μθ|/σθ; 2, 6); BB → smf(t; μθ−6σθ, μθ−2σθ) - s_mis = 0.7 if only one cluster present and it is not AB; 0.95 if two clusters present; else 1 7. prelim = s_ap · G · s_θ · s_r · s_mis (NaN if G or t is NaN). 8. **GenCall score** GC = prelim^0.35 · gbellmf(prelim; a=0.504, b=1.71, c=1.08); NaN → 0. 9. GC < cutoff (default **0.15** [DOC: Infinium genotyping data analysis tech note]) → NC. 10. Call rate = NumCalls / (NumSNPs − NumIntensityOnly). Note: the 2004 tech note mentions a "DNA score" penalty; current open reimplementation applies none (s_n = 1). [INF] ### 7.6 BAF, LRR, LogR deviation [LIT: Peiffer et al. 2006 Genome Res 16:1136] ``` if θ == θ_AB: BAF = 0.5 ; R_exp = R_AB elif θ < θ_AB: BAF = 0.5·(θ − θ_AA)/(θ_AB − θ_AA) ; R_exp = R_AA + (θ − θ_AA)(R_AA − R_AB)/(θ_AA − θ_AB) else: BAF = 0.5 + 0.5·(θ − θ_AB)/(θ_BB − θ_AB) ; R_exp = R_AB + (θ − θ_AB)(R_BB − R_AB)/(θ_BB − θ_AB) intensity-only loci: R_exp = R̄ (count-weighted mean, §6.4) ; BAF = NaN LRR = log2(R / R_exp) (R == 0 → −FLT_MAX in bit-exact mode) ; θ NaN → BAF = LRR = NaN BAF clamped to [0, 1] LogRDev = population SD of LRR over loci that are autosomal (chrom not X/Y/XY/MT), not intensity-only, finite ``` Optional (clean mode / CN module): cluster re-centering on the batch with shrinkage toward EGT means (weight 0.2, as gtc2vcf `--adjust-clusters`), GC-wave correction (§10). ### 7.7 Sample statistics - GC10 / GC50 = P_10 / P_50 of GenCall scores over called, non-intensity-only loci (Illumina percentile). - Intensity percentiles P05/P50/P95 of raw X and raw Y over all loci (stored as uint16). - Autosomal call rate = call rate restricted to autosomes (definition to confirm, §20). ### 7.8 Sex estimate [OSS: idat2gtc `estimate_gender` :3416-3482] Only loci with GenTrain score ≠ 0. - v1 (legacy): if > 20 X loci called: F if X-het-rate > 0.1 else M. - **v2 (default):** require autosomal call rate > threshold (0.97 historically; 0.0 in current gtc2vcf default), estimated on up to 10 000 random autosomal loci; need ≥ 20 Y loci and ≥ 100 autosomal loci. Then: median R(Y) > 0.3 → **M**; else if X loci < 20 → **F**; else median R(X) > 0.9 → **F**; else **U**. - SimurgArray adds an independent literature method as a second opinion: X-heterozygosity F-statistic + Y intensity (Laurie et al. 2010), reported separately. ### 7.9 Locus summary statistics [OSS: BeadArrayFiles `examples/locus_summary.py`] Call_Freq = calls/(calls+nocalls); AA/AB/BB freq over calls; p_A = (2AA+AB)/(2·calls); Minor_Freq = min(p_A, 1−p_A); Het excess: q = Minor_Freq, p = 1−q, k = (0.013/q)^4, dh = (AB/calls + k)/(2pq + k) − 1; ChiTest_P100 = 2Φ(dh;0,0.1) if dh<0 else 2(1−Φ(dh;0,0.1)); no calls → (0,0); monomorphic → (1,0). Per-locus θ/R means and SDs per genotype; GC10/GC50 per locus. Suggested locus filters [DOC: Infinium genotyping data analysis tech note, Table 1]: Call_Freq < 0.97; Cluster_Sep ≤ 0.3; AA/AB/BB R mean ≤ 0.2; 10%GC ≤ 0.3; Het excess > 0.2 or < −0.3; AB θ mean < 0.2 or > 0.8. --- ## 8. Strand, alleles & VCF export ### 8.1 Strand conventions [DOC: Illumina "TOP/BOT strand and A/B allele" tech note] - **TOP/BOT**: unambiguous SNP — if alleles are {A, C|G} → TOP (A=A, B=C/G); if {T, C|G} → BOT (A=T, B=C/G). Ambiguous [A/T]/[C/G] — "sequence walking": compare pairs (n−1,n+1), (n−2,n+2)… until one pair is unambiguous (one of A/T, the other C/G); if the A/T base is 5′ of the SNP → TOP, if 3′ → BOT. [A/T] TOP: A=A,B=T; BOT: A=T,B=A. [C/G] TOP: A=C,B=G; BOT: A=G,B=C. - **PLUS/MINUS**: relative to the reference genome (+ = 5′ end at p-arm tip). **RefStrand** (BPM v8/CSV) gives it directly; else derive by aligning SourceSeq/TopGenomicSeq flanks to the reference FASTA. - **Indels**: SNP column `[I/D]`/`[D/I]`; IlmnStrand is PLUS/MINUS; resolve inserted sequence and REF/ALT from the CSV SourceSeq against the FASTA (left-normalize). Without CSV → skip indels with a warning. - Base calls in GTC are on the TOP strand. ### 8.2 Harmonization engine For every locus produce an **audit row**: `locus, original_alleles, strand_in, transformed_alleles, strand_out, ref, alt, flipped(bool), swapped(bool), reason, build, confidence`. Never flip silently. Ambiguous A/T–C/G loci without RefStrand are marked `confidence=low` and optionally dropped (`--drop-ambiguous`). Build liftover is out of scope (users supply a manifest for the right build); chromosome naming mapped via a configurable alias table (`1`↔`chr1`, `XY`→PAR, `MT`↔`chrM`). ### 8.3 VCF writer - One VCF per sample (default) or multi-sample BCF (`--merge`), bgzip + tabix/CSI. - REF from FASTA; ALT from A/B alleles mapped to PLUS strand; genotype indices accordingly. - FORMAT: `GT` (ploidy-aware: haploid on Y, MT and non-PAR X for males), `GS` (GenCall score), `BAF`, `LRR` (BAF/LRR re-oriented so that BAF refers to the ALT allele → BAF_vcf = 1 − BAF when A allele = ALT), optional `X,Y,NORMX,NORMY,THETA,R`. - INFO: `ASSAY_TYPE`, `NORM_BIN`, `GENTRAIN`, `CLUSTER_SEP`, `ILMN_STRAND`, `PROBE_IDS`. - **Multi-probe sites** (duplicates, multi-allelic): enumerate candidate genotypes from all queried alleles; remove candidates inconsistent with any probe's call, respecting design limits (Infinium II cannot distinguish A/T or C/G); exactly one survivor → call; otherwise NC. Merged `GS` = min of probe scores; merged records carry no BAF/LRR. Probe IDs joined in `ID` (`,` SNPs; `;` indels). Option `--unsquash-duplicates` writes each probe separately. (Reference behaviour: Illumina GTCtoVCF, Apache-2.0.) - Header: `##source=SimurgArray `, `##reference`, `##simurg_manifest`, `##simurg_cluster_file`, provenance hash. ### 8.4 PLINK `.bed/.bim/.fam` (plink 1.9) and `.pgen/.pvar/.psam` (plink 2, via pgenlib if needed). Sex from sheet or estimate. --- ## 9. QC module ### 9.1 Control-probe metrics [DOC: Illumina BeadArray Controls Reporter / Infinium controls training guide] bkg_grn = max green intensity of Extension(A|T) controls; bkg_red = max red intensity of Extension(C|G); x = 3000. | Metric | Formula | |---|---| | Restoration | Restoration green / (bkg_grn + x) | | Staining green / red | Biotin(High)/Biotin(Bgnd) ; DNP(High)/DNP(Bgnd) | | Extension green | min Ext(C,G) green / max Ext(A,T) green | | Extension red | min Ext(A,T) red / max Ext(C,G) red | | Hybridization High/Med, Med/Low | green ratios | | Target removal | (bkg_grn + x)/TargetRemoval green | | Stringency | String(PM) red / String(MM) red | | Non-specific binding green / red | (bkg + x)/max NSB(Bgnd) | | Non-polymorphic green / red | min NP(C,G)/max NP(A,T) green ; min NP(A,T)/max NP(C,G) red | Example thresholds (configurable YAML): staining 5, extension 5, target removal 5, hyb 1.25, NP 5, stringency 1.5, NSB 5. ### 9.2 Sample QC Call rate, autosomal call rate, LogRDev, GC10/GC50, sex estimate vs reported sex, heterozygosity rate (±3 SD flag), contamination (§9.3), intensity percentiles. Default flags: call rate < 0.98, LogRDev > 0.20. ### 9.3 Contamination — BAF regression [LIT: Jun et al. 2012 AJHG 91:839; OSS: gtc2vcf `BAFregress.c` MIT] Sites with population B-allele frequency p ∈ [0.1, 0.9] (from the cohort or a reference panel). AA calls: y = BAF − median_AA(BAF), x = p. BB calls: y = 1 − (BAF − (median_BB(BAF) − 1)), x = 1 − p. Contamination α̂ = OLS slope m = (nΣxy − ΣxΣy)/(nΣx² − (Σx)²) (median adjustment needs ≥ 5 genotypes/cluster). Oracle: verifyIDintensity (GPL, subprocess). ### 9.4 Identity & relatedness KING-robust kinship (Manichaikul 2010) on high-quality autosomal SNPs; duplicate/replicate detection (concordance ≥ 0.99); pedigree check vs sample sheet Parent1/Parent2 (Mendelian error rate). Oracles: plink2, somalier (MIT), peddy (MIT). ### 9.5 QC report Self-contained HTML (inline JS/CSS, no network): per-sample table with PASS/FAIL, control-probe heatmap, call-rate vs LogRDev scatter, sex plot (X vs Y median R), intensity distributions, per-chip layout. Plus TSV/XLSX export. --- ## 10. Copy number, LOH & mosaicism module (v2) ### 10.1 Pre-processing - GC-wave correction [LIT: Diskin et al. 2008 NAR 36:e126]: regress LRR on probe-local GC content (1 kb / 100 kb / 1 Mb windows computed from FASTA) + optional replication-timing covariate; report GCWF. - Batch reference: per-probe median LRR across a reference set (≥ 24 samples) subtracted ("batch-normalized LRR"). - Population B-allele frequency (PFB) per locus from the cohort or reference panel. ### 10.2 Genome-wide CNV/LOH HMM [LIT: Wang et al. 2007 Genome Res (PennCNV); DOC: Illumina CNV algorithms tech note (cnvPartition)] States: CN0, CN1, CN2, CN2-LOH, CN3, CN4(+) and mosaic states (CN1/CN3/LOH with cell fraction f). Emissions: LRR ~ state-specific Gaussian mixed with uniform outlier component; BAF ~ mixture over genotypes with weights from PFB (expected BAF b/n for n copies, b B-alleles; mosaic: 0.5 ± f/2·…). Transitions depend on inter-marker distance d: P(change) ∝ 1 − exp(−d/D). Reference parameters (cnvPartition): LRR means CN0 −5 (SD 2), CN1 −0.45, CN2 0, CN3 0.3, CN4 0.75 (SD 0.18); BAF SD 0.03. Post-processing: minimum probes/size per event type (defaults: DEL 10 probes / 25 kb; DUP 10 / 50 kb; LOH 500 / 3 Mb; mosaic same), merge nearby segments, QUAL = −10·log10(1 − P(state)) capped at 60, effective size excludes gaps > 150× median probe spacing, sample-level flag when LogRDev > 0.2, only whole-chromosome events when > 0.3. Alternative segmentation for validation: cnvPartition Z-statistic (window 4/8/16/32, |Z| > 6) and homozygosity detector (χ² = N·E_het/(N−E_het), cut-off 23.5, ≥ 50 SNPs). ### 10.3 Mosaicism [LIT: Loh et al. 2018 Nature 559:350 (MoChA); OSS: mocha MIT] Phase genotypes (external phasing: SHAPEIT5/Eagle), then HMM over phased BAF imbalance |BAF−0.5| with consistent haplotype sign; classify loss/gain/CN-LOH with LRR; estimate cell fraction. ### 10.4 Annotation Cytoband, genes (RefSeq/Ensembl GTF), ISCN-2020-like string, overlap with known CNV databases (DGV) — all from public files fetched by script. --- ## 11. Pharmacogenomics module (v2) - **Allele definitions**: PharmVar + PharmGKB/CPIC allele definition tables (public downloads); build a versioned SQLite `pgx_db` mapping array probes → HGVS/genomic variants (via VCF positions). - **Star-allele calling**: VCF (GRCh38) → candidate haplotypes; score diplotypes by matched/mismatched/missing defining positions + genotype quality (GS) + optional population frequency prior; return ranked solutions with probability-like score and "raw score" without prior. Oracles: **PharmCAT** (MPL-2.0, subprocess) and **PyPGx** (MIT, supports chip VCFs). - **CYP2D6 (and others) copy number** [HYP design]: per-sub-region (e.g. CYP2D6 exon 9, intron 2, 5′ region) median of batch- and GC-normalized R/LRR; Gaussian mixture over CN 0–5 trained on ≥ 24 samples (≥ 96 for training a reusable model); combine sub-region CNs into structural haplotypes (*5 deletion, xN, *36/*13 hybrids); check CN-consistent BAF (1/3, 2/3). Literature: Asterix (Lanting et al. 2025 AJHG), Stargazer (Lee et al. 2019). - **Phenotype**: CPIC / DPWG lookup tables (activity scores for CYP2D6, CYP2C9, DPYD). - Truth for validation: star alleles and CN derived from 1000 Genomes 30x WGS (PyPGx/PharmCAT on WGS) for the same individuals genotyped on arrays. --- ## 12. Cluster-file training module — ArrayTrain (v2) Goal: generate our own EGT-format cluster files from ≥ 100–2000 samples (no vendor EGT needed). Explicitly **not** a reproduction of the vendor's clustering; documented and validated on its own. Model per locus in (θ, R) [LIT: Illuminus (Teo 2007), CRLMM (Carvalho 2007), GenTrain 2/3 tech notes]: 1. Candidate models: 1, 2 or 3 clusters (AA/AB/BB assignments ordered by θ), bivariate Student-t (robust) mixture fitted by EM with θ-ordering constraints; sex-aware on X (males haploid) and Y/MT haploid. 2. Model selection score = log-likelihood − BIC penalty + HWE log-likelihood term + separation term; choose best. 3. Missing clusters: impute positions from empirical relationships learned across well-clustered loci of the same assay type / norm bin (regression θ_AB ≈ f(θ_AA, θ_BB), R ratios) — replaces the vendor's neural network. 4. Scores: our own `cluster_sep` (min standardized θ gap between adjacent clusters, squashed to 0–1) and `train_score` (posterior-weighted fit quality) written into the EGT fields; document the formulas. 5. Writer: EGT v3, data-block version 9, readable by BeadArrayFiles/Picard/gtc2vcf. Validation: call concordance vs vendor EGT calls (Omni2.5, 2141 samples) and vs WGS truth (non-ref discordance), per MAF bin; rare-variant performance compared with zCall-style rescue. --- ## 13. SimurgArray output specification > **See `SIMURGARRAY_FUNCTIONAL_SPEC.md` §4 and §6–§9** (items OUT-*, QC-*, CYTO-*, PCN-*, PGX-*, MET-*) for the > exact schemas, headers and column orders. Where file names below differ from the spec, the spec wins: > - `.snv.vcf.gz`, `.cyto.cnv.vcf.gz`, `.pgx_cn.vcf.gz`; > - `.cyto.json` and `.pgx.json`; > - `controls_raw.tsv` and `controls_qc.tsv`, written by `genotype` itself; > - `qc_report_YYYY_MM_DD.html` and `qc_table_YYYY_MM_DD.xlsx`; > - the PGx batch TSVs. Per sample (in `/samples//`): | File | Content | |---|---| | `.gtc` | GTC v5 (interoperable) | | `.vcf.gz` + `.tbi` | GT:GS:BAF:LRR (§8.3) | | `.baf.bedgraph.gz`, `.lrr.bedgraph.gz` | IGV tracks | | `.calls.parquet` | full per-locus table (§5) | | `.cnv.vcf.gz` (v2) | `//` with CN, MF, QUAL | | `.pgx.json` (v2) | diplotypes, CN, phenotype | Per batch (in `/`): | File | Content | |---|---| | `sample_summary.tsv` / `.json` | sample_id, name, call_rate, autosomal_call_rate, logr_dev, gc10, gc50, sex_est, sex_reported, contamination, flags | | `controls_raw.tsv`, `controls_qc.tsv` | §9.1 | | `locus_summary.tsv` | §7.9 | | `cohort.bcf` / plink files | optional merged outputs | | `qc_report.html` | §9.5 | | `run_manifest.json` | provenance (§15) | | `audit/*.parquet` | harmonization + QC audit trail | --- ## 14. CLI specification (`simurg`) > **See `SIMURGARRAY_FUNCTIONAL_SPEC.md` §2–§3** (items CLI-*, IN-*) for: > - the full command tree, and every option with its type, default and validation; > - precedence (CLI > YAML > profile > defaults); > - the `--log-json` schema and the message-code catalogue; > - overwrite/resume policy and exit codes (adds 4 = no sample succeeded, and 130/143 = interrupted); > - sample-sheet and input validation rules. > > The summary below is kept for orientation. ``` simurg inspect # print header/TOC of IDAT/BPM/EGT/GTC simurg ingest --idats DIR|--sheet F --bpm F [--csv F] --out DIR # → sample_intensity parquet simurg normalize --in DIR --bpm F [--norm-version 1.2.0|1.1.2] [--exact] simurg genotype --idats DIR|--sheet F --bpm F --egt F [--cutoff 0.15] [--threads N] --out DIR # IDAT→GTC+calls simurg gtc-to-vcf --gtcs DIR --bpm F --csv F --fasta F [--merge] [--unsquash-duplicates] [--drop-ambiguous] simurg export plink|bedgraph|finalreport ... simurg qc controls|samples|loci|contamination|relatedness|report ... simurg train --gtcs DIR|--calls DIR --bpm F --out cluster.egt # v2 simurg cn call|train|annotate ... # v2 simurg pgx copy-number|star-allele|phenotype ... # v2 simurg pipeline run config.yaml [--resume] [--dry-run] simurg compare a.gtc b.gtc | a.vcf.gz b.vcf.gz [--tolerance ...] # field-level diff (oracle testing) simurg synth --n-samples ... --out DIR # synthetic data ``` Common flags: `--threads`, `--log-json`, `--log-level`, `--seed`, `--config`, `--version`. Exit codes: 0 ok, 2 input validation error, 3 partial failure (some samples failed; listed in summary), 1 internal. Behaviour: never silently overwrite (require `--force`), resumable via content-addressed step cache, progress bars with ETA, per-sample failure isolation. Example `config.yaml`: ```yaml inputs: {sheet: data/gda_demo/SampleSheet.csv, idat_dir: data/gda_demo/idats} reference: {bpm: ..., csv: ..., egt: ..., fasta: data/ref/GRCh38.fa} genotype: {cutoff: 0.15, norm_version: "1.2.0", exact: true} outputs: {gtc: true, vcf: {per_sample: true, merge: true}, plink: true, bedgraph: true} qc: {call_rate_min: 0.98, logr_dev_max: 0.20, controls_yaml: configs/controls_genotyping.yaml} runtime: {threads: 16, seed: 42} ``` --- ## 15. Provenance & reproducibility `run_manifest.json` records: SimurgArray version + git commit, Python/Rust/library versions, OS/CPU, command line, full resolved config, random seed, SHA-256 of every input (IDAT, BPM, CSV, EGT, FASTA) and output, manifest name + genome build, cluster-file name, per-step start/end timestamps, thread count, warnings. Determinism: fixed sort orders, no unordered-set iteration in numerics, deterministic parallel reductions (per-sample independence; cohort reductions in fixed chunk order), float policy documented per field. Regression test: two runs with different thread counts must give byte-identical outputs (excluding timestamps, which are isolated in the manifest). --- ## 16. Testing & validation strategy ### 16.1 Test layers unit (every function) · property-based (Hypothesis: round-trips, invariants such as θ∈[0,1], BAF∈[0,1], permutation invariance of per-sample steps) · malformed-input (truncated/garbage files, wrong magic, bad TOC offsets) · integration (IDAT→VCF on small fixtures) · oracle (field diffs vs reference tools) · regression (golden hashes) · performance (pytest-benchmark) · reproducibility (thread-count invariance). ### 16.2 Oracles (compatibility references) | Oracle | What it checks | Access | |---|---|---| | **Illumina-produced GTCs** shipped with 1000 Genomes Omni2.5 (2141 samples) | normalization transforms, genotypes, scores, BAF/LRR, sample stats — the vendor ground truth | open FTP | | GenomeStudio Final Report (1000G Sanger set) | X/Y, θ/R, GType, BAF/LRR | open FTP | | gtc2vcf `idat2gtc` + `gtc2vcf` (MIT) | full pipeline, any array | build locally | | BeadArrayFiles, Picard | format parsing, normalized intensities | open source | | illuminaio (R) | IDAT parsing | open source (GPL, subprocess) | | plink2, somalier, verifyIDintensity, PennCNV, PharmCAT, PyPGx | QC, CN, PGx | subprocess | ### 16.3 Comparison categories & tolerances Each field compared per locus class (Inf I / Inf II, SNP / indel, autosome / X / Y / MT / PAR, intensity-only): **EXACT** (bitwise), **NEAR** (within tolerance), **DISCREPANCY**, **UNKNOWN** (oracle lacks field). Tolerances (bit-exact mode targets): raw X/Y exact; transform floats ≤ 2 ULP float32; θ, R, BAF, LRR ≤ 1e-5 abs; GenCall score ≤ 1e-5 abs; genotypes exact (target 100%, gate ≥ 99.999%); call rate/LogRDev/GC10/GC50 ≤ 1e-6; sex exact. Report every discrepancy with locus/sample IDs — never hide. ### 16.4 Synthetic data generator (`simurg synth`) Generates canonical intensities *and* valid IDAT pairs + matching toy BPM/EGT with known truth: genotype clusters with configurable separation/noise, missing beads, low-quality samples, outliers, dye bias (channel scale/offset/ crosstalk → tests normalization recovery), batch effects, allele imbalance, contamination fraction α, CN segments (del/dup/LOH/mosaic f), sex chromosomes. Because truth is exact, it validates algorithms independently of vendor data and is the only data safe to commit as fixtures. ### 16.5 Scientific validation metrics Genotype concordance, non-reference discordance, no-call rate, sensitivity/specificity/precision/F1 per genotype class and MAF bin vs NYGC 30x WGS and GIAB v4.2.1; Mendelian error rate (1000G trios, GDA demo trios); replicate concordance (GSE165845 triplicates); strand/allele error rate; batch effect (per-chip call-rate/LogRDev); CNV precision/recall vs 1000G WGS CNV/SV callsets; PGx diplotype concordance vs WGS-derived calls. --- ## 17. Datasets (with URLs) and fetch plan Vendor files are for personal research use only and are **never committed**. `scripts/fetch_data.sh ` downloads, records SHA-256 in `data/CHECKSUMS.sha256`, and unpacks. ### Tier P0 — small, start here (~1.4 GB; verified 2026-09-27) | Dataset | URL | Notes | |---|---|---| | **GSA v4 PGx demo IDATs (36 samples, 72 IDATs)** | https://support.illumina.com/content/dam/illumina-support/documents/documentation/chemistry_documentation/infinium_assays/infinium-gsa-pgx/idats_Demo_36.zip | 872 MB; primary small fixture. BPM/CSV/EGT from the same folder (`GSA-PGx-48v4-0_20079540_E2_{bpm,csv,ClusterFile}.zip`) | | Bioconductor IlluminaDataTestFiles | https://bioconductor.org/packages/IlluminaDataTestFiles | Artistic-2.0 → **may be committed** as IDAT parser fixtures; includes 2 encrypted IDATs (negative tests) | | GDA v1 BPM / CSV / EGT + demo sample sheet | https://support.illumina.com/array/array_kits/infinium-global-diversity-array/product-files.html | second design (1.9 M loci). Note: the "GDA Demo 36" bundle is a **GenomeStudio project, not IDATs** | | GRCh37 & GRCh38 FASTA | e.g. GRCh38 no-alt analysis set (NCBI) | contigs named per manifest | (Bioconductor `hapmap370k` is no longer in the current release.) Oracle GTCs for P1–P4 are generated locally with gtc2vcf `idat2gtc` (MIT) from the GSA demo IDATs — see `docs/validation/P1.md`. ### Tier P1 — validation core (~170 GB) | Dataset | URL | Notes | |---|---|---| | 1000G Omni2.5 Broad: IDAT **+ Illumina GTC**, 2141 samples + BPM/EGT/CSV | ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/supporting/hd_genotype_chip/broad_intensities/ | 107 GB; main bit-exact oracle; 337 trios | | 1000G Omni2.5-8 Sanger, 177 samples + GenomeStudio FinalReport | …/hd_genotype_chip/sanger_intensities/ | 6.2 GB + 11 GB report | | 1000G combined Omni VCF | …/hd_genotype_chip/ALL.chip.omni_broad_sanger_combined.20140818.snps.genotypes.vcf.gz | reference calls | | GEO GSE165845 — GSA, 120 1000G samples × 3 replicates | https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE165845 | 45 GB IDAT + VCF | | NYGC 1kGP 30x (3202 samples) | https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/20220422_3202_phased_SNV_INDEL_SV/ | truth | | GIAB v4.2.1 HG001 | https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/release/NA12878_HG001/NISTv4.2.1/GRCh38/ | truth | ### Tier P2 — breadth & scale GEO GSE96909 (NA12878 on 9 designs, 4.3 GB, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE96909) · GDA PGx demo 480 samples (https://s3.amazonaws.com/webdata.illumina.com/downloads/productfiles/global-diversity-array-pgx/GDAPGx_Demo_idats.zip, 16.5 GB) · OneK1K GSA 1104 samples (GEO GSE196829) · GSE209804 low-input DNA · GSE19247 CytoSNP pedigrees · GEO sweep: 254 genotyping series with IDATs (E-utilities query `IDAT[Supplementary Files] AND gse[Entry Type] AND ("SNP genotyping by SNP array"[DataSet Type] OR "genome variation profiling by SNP array"[DataSet Type])`) for parser fuzzing · 1000G phase 3 SV map (https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/integrated_sv_map/) · 1000 Genomes 30x SV/CNV callsets (NYGC release above) as CNV truth · dbSNP VCF (https://ftp.ncbi.nlm.nih.gov/snp/latest_release/VCF/) · strand files (https://www.chg.ox.ac.uk/~wrayner/strand/ — manual browser download). ### Reference repositories (`scripts/bootstrap_refs.sh`, shallow clones into `refs/`) ``` https://github.com/Illumina/BeadArrayFiles # BSD-2 — formats, GTC spec PDF in docs/ https://github.com/Illumina/GTCtoVCF # Apache-2.0 — VCF conversion logic https://github.com/freeseek/gtc2vcf # MIT — idat2gtc, gtc2vcf, BAFregress https://github.com/freeseek/mocha # MIT — mosaic CNV https://github.com/broadinstitute/picard # MIT — src/main/java/picard/arrays https://github.com/WGLab/PennCNV # Unlicense — CNV HMM https://github.com/sbslee/pypgx # MIT — PGx https://github.com/PharmGKB/PharmCAT # MPL-2.0 — PGx oracle https://github.com/HenrikBengtsson/illuminaio # GPL-2 — oracle only, do not copy https://github.com/gjun/verifyIDintensity # GPL-3 — oracle only, do not copy ``` --- ## 18. Phased roadmap with acceptance criteria | Phase | Work | Acceptance criteria (gate) | |---|---|---| | **P0 Setup** (≈1 wk) | repo skeleton, uv project, CI (ruff, mypy, pytest), fetch scripts, refs bootstrap, `THIRD_PARTY_NOTICES.md` | `pytest` green on empty suite; P0 data fetched with checksums | | **P1 Parsers** (≈2–3 wk) | `binutil`, IDAT, BPM, CSV, EGT, GTC read + GTC write, sample sheet, `simurg inspect` | all fields equal to BeadArrayFiles/illuminaio on P0 + Omni data; GTC read→write round-trip byte-identical; fuzz tests never crash (clean errors) | | **P2 Data model** (≈1 wk) | Arrow schemas, raw X/Y assembly, norm-bin rule | raw X/Y **exactly** equal to GTC TOC 1000/1001 of 1000G Illumina GTCs for 100% loci × ≥50 samples | | **P3 Normalization** (≈3–4 wk) | numerics (Appendix A), fit v1.2.0 (+1.1.2), apply, polar | transforms vs Illumina GTC TOC 400 within 2 ULP for ≥99.9% bins; normalized X/Y within 1e-5 | | **P4 Genotyping** (≈3 wk) | GenCall-equivalent caller, score, BAF/LRR, LogRDev, GC10/50, sex, percentiles, GTC v5 writer | vs Illumina GTCs: genotypes ≥99.999% identical, scores ≤1e-5, sample stats ≤1e-6, sex 100%; md5-identical GTC vs gtc2vcf output (timestamps masked) | | **P5 QC** (≈2 wk) | control metrics, sample/locus summaries, contamination, relatedness, HTML report | control metrics equal reference formulas on synthetic controls; 1000G sex 100% vs pedigree; replicate/duplicate detection 100% on demo | | **P6 Outputs** (≈2 wk) | harmonization + audit, VCF/BCF, PLINK, bedGraph, final-report export | `bcftools` validates; vs gtc2vcf/GTCtoVCF VCF: 0 unexplained genotype differences; strand audit complete | | **P7 Differential-testing harness** (≈1 wk) | `simurg compare`, dashboard of EXACT/NEAR/DISCREPANCY/UNKNOWN | nightly run on P0 + Omni subset produces report | | **P8 Scientific validation** (≈2 wk) | vs 30x WGS/GIAB, trios, replicates | written validation report with all metrics of §16.5 | | **v1.0 release** | docs, examples | | | **P9 Performance** (≈3–4 wk) | Rust core (IDAT, normalization, calling) via PyO3, mmap, rayon threads | outputs identical to Python reference; ≥ 20 samples/min/core for GSA-size arrays; linear scaling to 32 threads | | **P10 Packaging** | pip wheel, conda recipe, Docker image, user guide | clean install + tutorial works on fresh VM | | **v2 modules** | ArrayTrain (§12), CN/LOH/mosaic (§10), PGx (§11), methylation control QC | each with own validation report vs truth sets | Rule: **do not proceed while a fundamental assumption is untested** — log it in `docs/known_unknowns.md`. --- ## 19. Performance plan Expected per-sample cost (GSA ≈ 650k loci; Omni2.5 ≈ 2.4M): IDAT parse (2 × ~10 MB, zero-copy) ≪ 0.1 s; normalization (robust fits on a few thousand points per bin, NN sweeps) ~0.1–0.5 s; calling O(n_loci) ~0.2 s vectorized; VCF writing often dominant (compression) → use htslib multithreaded bgzf. Target: ≥ 1 sample/s/core in Rust. Measure: samples/s, loci/s, GB/h, peak RSS/sample, CPU utilisation, scaling 1→32 threads, disk I/O. Workloads: 10, 100, 1 000, 10 000 samples (synthetic replicate expansion for 100 000). Acceleration verdict: CPU SIMD + multithreading is sufficient for everything except possibly ArrayTrain at ≥ 10k samples (EM over millions of loci) — consider GPU (CuPy/JAX) there only after profiling. FPGA not justified. --- ## 20. Known unknowns & planned experiments | ID | Unknown | Experiment | |---|---|---| | U1 | Raw X is red (open implementations) vs "green" (GTC PDF) for Infinium II | compare assembled X/Y with Illumina GTC raw arrays (P2) | | U2 | IDAT MidBlock element type (int32 vs uint8) | inspect real files from many chips (GEO sweep) | | U3 | EGT intensity-threshold version gate (dbv ≥ 7 vs == 9) | parse EGTs of different vintages; check calls with/without τ | | U4 | Exact quirk reproduction in normalization (NN binned search, robustfit0 adjfactor, zero-padding) | ablation: toggle each quirk, measure transform agreement with Illumina GTCs | | U5 | Which normalization version produced each public GTC set (1.1.2 vs 1.2.0) | fit both; pick best match per dataset (GTC autocall-version string) | | U6 | Autosomal-call-rate and LogRDev locus sets (PAR, XY, MT handling) | recompute from Illumina GTC arrays; compare to stored values | | U7 | Sex estimate thresholds per software generation | sweep on 1000G known-sex samples + synthetic X/Y R | | U8 | Score edge cases: zero intensities, NaN clusters, one-cluster loci, Omni5 bin overflow | synthetic single-locus perturbations | | U9 | GenTrain score / Cluster_Sep formulas (unpublished) | not required for v1 (read from EGT); ArrayTrain defines its own | | U10 | Best CN model for targeted PGx genes | compare GMM vs HMM vs Asterix-style on 1000G arrays with WGS CN truth | | U11 | Floating-point reproducibility across CPUs (x86 vs ARM, AVX) | run golden tests on both architectures | --- ## 21. Risk register | # | Risk | Likelihood | Impact | Mitigation | |---|---|---|---|---| | R1 | Bit-exactness fails on some quirk | M | M | keep "clean" + "exact" modes; document residual differences quantitatively | | R2 | Dependence on vendor EGT files | H | M | ArrayTrain (v2) | | R3 | Vendor data file terms | L (personal use) | M | never commit/redistribute; fetch script only | | R4 | GPL contamination of codebase | M | M | GPL tools only as subprocesses; licence check in CI (`pip-licenses`) | | R5 | Large data exceeds disk/RAM | M | M | streaming per sample, Parquet chunking, P1 subsets first | | R6 | Strand/allele errors in VCF | M | H | audit trail, WGS concordance by strand class, A/T–C/G flagged | | R7 | CN/PGx accuracy lower than expected | M | M | v2 scope, validate vs WGS truth, report honestly | | R8 | Non-determinism across threads/CPUs | M | M | deterministic reductions, golden tests on two architectures | | R9 | Scope creep | H | M | phase gates; v1 = genotyping + QC + VCF only | --- ## 22. References Illumina documentation (public PDFs): GenCall Data Analysis Software tech note (https://www.illumina.com/Documents/products/technotes/technote_gencall_data_analysis_software.pdf) · Infinium Genotyping Data Analysis (https://www.illumina.com/Documents/products/technotes/technote_infinium_genotyping_data_analysis.pdf) · GenTrain 2.0 (https://www.illumina.com/Documents/products/technotes/technote_gentrain2.pdf) · GenTrain 3.0 M-GL-01258 (https://www.illumina.com/content/dam/illumina/gcs/assembled-assets/marketing-literature/gentrain-tech-note-m-gl-01258/gentrain-tech-note-m-gl-01258.pdf) · TOP/BOT strand (https://www.illumina.com/documents/products/technotes/technote_topbot.pdf) · DNA copy number algorithms / cnvPartition (https://www.illumina.com/documents/products/technotes/technote_cnv_algorithms.pdf) · X-chromosome cluster generation (https://www.illumina.com/documents/products/technotes/technote_x_chromosome_snp_cluster_generation.pdf) · Custom cluster file creation M-GL-02142 · iSelect design (https://www.illumina.com/Documents/products/technotes/technote_iselect_design.pdf) · Infinium controls training guide · GTC format v5 PDF (BeadArrayFiles/docs). Patents (expired): US7035740, US7467117 (Kermani, normalization & clustering). Papers: Gunderson 2005 Nat Genet 37:549 · Steemers 2006 Nat Methods 3:31 · Peiffer 2006 Genome Res 16:1136 · Teo 2007 Bioinformatics 23:2741 (Illuminus) · Giannoulatou 2008 Bioinformatics 24:2209 (GenoSNP) · Carvalho 2007 Biostatistics 8:485, Ritchie 2009 Bioinformatics 25:2621 (CRLMM) · Ritchie 2011 BMC Bioinformatics 12:68 · Goldstein 2012 Bioinformatics 28:2543 (zCall) · Shah 2012 Bioinformatics 28:1598 (optiCall) · Wang 2007 Genome Res 17:1665 (PennCNV) · Diskin 2008 NAR 36:e126 · Colella 2007 NAR 35:2013 (QuantiSNP) · Loh 2018 Nature 559:350, Loh 2020 Nature 584:136 (MoChA) · Jun 2012 AJHG 91:839 (contamination) · Laurie 2010 Genet Epidemiol 34:591 (array QC) · Anderson 2010 Nat Protoc 5:1564 · Manichaikul 2010 Bioinformatics 26:2867 (KING) · Gaedigk 2018 CPT 103:399 (PharmVar) · Sangkuhl 2020 CPT 107:203 (PharmCAT) · Lee 2019 Genet Med 21:361 (Stargazer) · Lanting 2025 AJHG 112:1015 (Asterix) · Li 2021 Genome Res (GSE165845) · Haraksingh 2017 BMC Genomics 18:321 (GSE96909). --- ## Appendix A — numeric helper definitions (bit-exact mode) **Illumina percentile** (used for normalization, GC10/GC50): sorted v of length n, target p ∈ [0,100]: ``` i1 = floor(n·p/100); f = n·p/100 − i1; if f < 0.5: i1 −= 1 clamp i1 to [0, n−2]; x1 = 100(i1 + 0.5)/n; x2 = 100(i1 + 1.5)/n P = v[i1] + (v[i1+1] − v[i1])·(p − x1)/(x2 − x1) (clamp at ends) ``` (not numpy's default percentile!) **trimmeanK(v)** = mean of {v : P_{K/2}(v) ≤ v ≤ P_{100−K/2}(v)}. **Robust regression** (MATLAB `robustfit` equivalent): iteratively reweighted least squares, bisquare weights w = (1 − u²)² for |u| < 1 else 0, u = r_adj/(4.685·s), leverage-adjusted residuals r_adj = r/√(1−h), s = MAD-sigma = median(|r| of the p smallest-leverage-excluded residuals)/0.6745 (v1.2.0 "new madsigma"; v1.1.2 uses the median-centred version), max 50 iterations, relative tolerance 1e-6. `robustfit1` = with intercept; `robustfit0` = through origin; bit-exact quirk: its leverage adjfactor is effectively constant ≈ 100. (Reference: gtc2vcf idat2gtc.c:1861-1942.) **Fuzzy memberships** (MATLAB): ``` smf(x; a, b) = 0 x ≤ a = 2((x−a)/(b−a))² a < x ≤ (a+b)/2 = 1 − 2((x−b)/(b−a))² (a+b)/2 < x ≤ b = 1 x > b zmf(x; a, b) = 1 − smf(x; a, b) gbellmf(x; a, b, c) = 1 / (1 + |(x − c)/a|^(2b)) ``` **Nearest-neighbour sweep**: for each of 400 sweep points along an axis, the closest data point (Euclidean); exact brute force for < 2000 points; for ≥ 2000 the vendor's binned approximation (reproduce from gtc2vcf `nearest_neighbor.c` in bit-exact mode; exact KD-tree in clean mode).