kemalandClaude Opus 5.5
ee4ebaaa1e
Phase 0 deliverables: full technical specification and roadmap (SIMURGARRAY_MASTER_PLAN.md), CLAUDE.md, bootstrap script for reference open-source repos and tiered dataset fetch script with checksums. Co-Authored-By: Claude Opus 5.5 <[email protected]> Claude-Session: https://claude.ai/code/session_01814GKVzZVL3FuF3DeiZXhv
61 KiB
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 (human or Claude Code) 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
- How to use this document
- Scope, goals, non-goals
- Ground rules (licensing & provenance)
- Target platform & environment setup
- Architecture & repository layout
- Canonical data model
- File-format specifications (IDAT, BPM, CSV manifest, EGT, GTC, sample sheet)
- Core algorithms (exact math)
- Strand, alleles & VCF export
- QC module
- Copy number, LOH & mosaicism module
- Pharmacogenomics module
- Cluster-file training module (ArrayTrain)
- SimurgArray output specification
- CLI specification
- Provenance & reproducibility
- Testing & validation strategy
- Datasets (with URLs) and fetch plan
- Phased roadmap with acceptance criteria
- Performance plan
- Known unknowns & planned experiments
- Risk register
- References
- Appendix A —
CLAUDE.mdfor the new repository - Appendix B — first-session prompt for Claude Code
- Appendix C — numeric helper definitions (percentile, trimmean, robustfit, fuzzy memberships)
0. How to use this document
- 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)
- Read IDAT (v3 binary), BPM, CSV manifest, EGT, GTC (v3–v5), sample sheets.
- Assemble per-locus raw X/Y from Red/Green IDATs.
- Fit and apply the per-bin affine intensity normalization.
- Call genotypes from an existing cluster file (EGT) with a GenCall-equivalent score.
- Compute θ, R, BAF, LRR, LogR deviation, call rates, GC10/GC50, intensity percentiles, sex estimate.
- Control-probe QC metrics; sample- and locus-level QC tables; HTML QC report.
- 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.
- Allele/strand harmonization with a full audit trail.
- Deterministic, provenance-tracked, resumable batch pipeline; scalable to ≥10 000 samples.
1.2 Goals (v2.0)
- ArrayTrain: our own cluster training → EGT files (so SimurgArray is not dependent on vendor EGTs).
- Genome-wide CNV/LOH/mosaic calling from LRR/BAF (HMM), with ISCN-like annotation.
- Targeted PGx copy number (e.g. CYP2D6 sub-regions) and star-allele diplotype calling from array VCFs.
- Contamination (BAF regression) and relatedness/sample-identity QC.
- 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.
- 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.
- 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.
- 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
- Never decompile, disassemble or inspect vendor software binaries, and never bypass licence checks.
- 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.
- 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
# 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
# GitHub CLI (authenticate to avoid API rate limits)
sudo apt install -y gh && gh auth login
# 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 private Gitea repo https://gitea.yaylali.uk/kemal/SimurgArray.
# 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`
├── CLAUDE.md # Appendix A
├── 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 C)
│ │ └── 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)
├── bench/
├── 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:
<barcode>_<RxxCxx>_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.frombufferat 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 C); 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:
- Mirror each point so that x ≥ y.
- c = P_2(y); x ← x − c, y ← y − c.
- r = |x|+|y|, t = (2/π)·atan2(y, x).
- AA = {t < 0.1}. ȳ = trimmean20(y_AA) − 50.
- r̄ = trimmean20(r with NaNs removed and array padded back with zeros — quirk) − 2ȳ.
- offset_x = offset_y = c + ȳ; θ = 0; shear = 0; scale_x = scale_y = r̄.
(b) Otherwise — six-parameter fit:
- 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.
- 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 gtc2vcfnearest_neighbor.c.) - robust fits (bisquare, Appendix C): 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.
- 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).
- 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 C); θ = atan(m); rotate x′ = cosθ·x + sinθ·y, y′ = −sinθ·x + cosθ·y. - 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.
- 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.
- 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).
- Intensity-only locus → NC (counted separately). If r < τ → NC, score NaN.
- Floor: σR_k ← max(σR_k, 0.1).
- 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)).
- 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). - 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.
- 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
- prelim = s_ap · G · s_θ · s_r · s_mis (NaN if G or t is NaN).
- GenCall score GC = prelim^0.35 · gbellmf(prelim; a=0.504, b=1.71, c=1.08); NaN → 0.
- GC < cutoff (default 0.15 [DOC: Infinium genotyping data analysis tech note]) → NC.
- 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), optionalX,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 inID(,SNPs;;indels). Option--unsquash-duplicateswrites each probe separately. (Reference behaviour: Illumina GTCtoVCF, Apache-2.0.) - Header:
##source=SimurgArray <version>,##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_dbmapping 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]:
- 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.
- Model selection score = log-likelihood − BIC penalty + HWE log-likelihood term + separation term; choose best.
- 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.
- Scores: our own
cluster_sep(min standardized θ gap between adjacent clusters, squashed to 0–1) andtrain_score(posterior-weighted fit quality) written into the EGT fields; document the formulas. - 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
Per sample (in <out>/samples/<sample>/):
| File | Content |
|---|---|
<sample>.gtc |
GTC v5 (interoperable) |
<sample>.vcf.gz + .tbi |
GT:GS:BAF:LRR (§8.3) |
<sample>.baf.bedgraph.gz, <sample>.lrr.bedgraph.gz |
IGV tracks |
<sample>.calls.parquet |
full per-locus table (§5) |
<sample>.cnv.vcf.gz (v2) |
<DEL>/<DUP>/<LOH> with CN, MF, QUAL |
<sample>.pgx.json (v2) |
diplotypes, CN, phenotype |
Per batch (in <out>/): |
|
| 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)
simurg inspect <file> # 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:
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 <tier> downloads,
records SHA-256 in data/CHECKSUMS.sha256, and unpacks.
Tier P0 — small, start here (~3 GB)
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 C), 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 — CLAUDE.md for the new repository
# SimurgArray
Personal research toolkit for Illumina Infinium genotyping arrays (IDAT→GTC→VCF, QC, CN, PGx).
The full spec and roadmap is SIMURGARRAY_MASTER_PLAN.md — read the relevant section before changing code.
## Rules
- Work phase by phase (§18). Don't start a phase until the previous gate passes.
- Science modules (signal/, genotype/, cn/, qc/, pgx/) are pure functions; no I/O inside them.
- Every algorithm cites its source (§7 tags) in its docstring.
- MIT/BSD/Apache/Unlicense code may be ported with attribution in THIRD_PARTY_NOTICES.md.
GPL tools (plink2, illuminaio, verifyIDintensity) = subprocess oracles only; never copy their code.
- Never commit vendor data (BPM/CSV/EGT/IDAT) — only tests/fixtures/ (Artistic-2.0 Bioconductor files, synthetic data).
- Never hide oracle discrepancies; record them in docs/known_unknowns.md.
- Default mode is bit-exact vs Illumina-produced GTCs; keep a clean-math mode too.
## Commands
- uv sync ; uv run pytest -q ; uv run ruff check . ; uv run mypy src
- bash scripts/fetch_data.sh p0|p1|p2 ; bash scripts/bootstrap_refs.sh
Appendix B — first-session prompt for Claude Code (on the new machine)
Read
SIMURGARRAY_MASTER_PLAN.mdfully. We are starting Phase P0 (Setup) then P1 (Parsers).
- Create the repository skeleton of §4.2 with a uv-managed
pyproject.toml(Python 3.12), ruff, mypy, pytest, and a GitHub Actions CI workflow.- Write
scripts/bootstrap_refs.sh(shallow clones of §17 repos intorefs/) andscripts/fetch_data.sh(tiers p0/p1/p2, sha256 recording, resumable downloads with aria2c/wget).- Run bootstrap + fetch p0.
- Implement
formats/binutil.pyandformats/idat.pyper §6.1 with tests (Bioconductor fixtures, malformed-file tests, property tests), thensimurg inspect. Stop after the IDAT reader passes the P1 criteria for IDAT and report results.
Appendix C — 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).