From ee4ebaaa1e4f81b54ca0e231c284d355ed01895c Mon Sep 17 00:00:00 2001 From: Kemal Yaylali Date: Sun, 27 Sep 2026 15:20:06 +0100 Subject: [PATCH] SimurgArray: master plan, repo scaffolding and data/ref fetch scripts 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 Claude-Session: https://claude.ai/code/session_01814GKVzZVL3FuF3DeiZXhv --- .gitignore | 18 + CLAUDE.md | 17 + README.md | 11 + SIMURGARRAY_MASTER_PLAN.md | 1004 ++++++++++++++++++++++++++++++++++++ scripts/bootstrap_refs.sh | 24 + scripts/fetch_data.sh | 65 +++ 6 files changed, 1139 insertions(+) create mode 100644 .gitignore create mode 100644 CLAUDE.md create mode 100644 README.md create mode 100644 SIMURGARRAY_MASTER_PLAN.md create mode 100755 scripts/bootstrap_refs.sh create mode 100755 scripts/fetch_data.sh diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..3f56cfc --- /dev/null +++ b/.gitignore @@ -0,0 +1,18 @@ +# vendor / large data — never commit +data/ +refs/ +_phase0_archive/ +*.idat +*.idat.gz +*.bpm +*.egt +*.gtc +# python +__pycache__/ +.venv/ +*.egg-info/ +.pytest_cache/ +.mypy_cache/ +.ruff_cache/ +# rust +target/ diff --git a/CLAUDE.md b/CLAUDE.md new file mode 100644 index 0000000..345f096 --- /dev/null +++ b/CLAUDE.md @@ -0,0 +1,17 @@ +# 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 diff --git a/README.md b/README.md new file mode 100644 index 0000000..420dfdf --- /dev/null +++ b/README.md @@ -0,0 +1,11 @@ +# SimurgArray + +Personal research toolkit for Illumina Infinium genotyping BeadChip data +(IDAT → normalization → genotype calling → BAF/LRR → QC → VCF/PLINK; later CN/LOH, PGx, cluster training). + +Start with **SIMURGARRAY_MASTER_PLAN.md** (full specification and roadmap). + +```bash +bash scripts/bootstrap_refs.sh # clone reference open-source repos into refs/ +bash scripts/fetch_data.sh p0 # download small test datasets into data/ +``` diff --git a/SIMURGARRAY_MASTER_PLAN.md b/SIMURGARRAY_MASTER_PLAN.md new file mode 100644 index 0000000..541c4ec --- /dev/null +++ b/SIMURGARRAY_MASTER_PLAN.md @@ -0,0 +1,1004 @@ +# 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 + +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 — `CLAUDE.md` for 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) +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 +# 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**. +```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` +├── 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: `__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 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: +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 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. +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 C); θ = 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 + +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`) + +``` +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 (~3 GB) +| Dataset | URL | Notes | +|---|---|---| +| Bioconductor IlluminaDataTestFiles | https://bioconductor.org/packages/IlluminaDataTestFiles | Artistic-2.0 → **may be committed** as IDAT parser fixtures | +| Bioconductor hapmap370k | https://bioconductor.org/packages/hapmap370k | Artistic-2.0; 40 HapMap IDATs | +| GDA v1 Demo 36 (IDATs) | https://webdata.illumina.com/downloads/productfiles/global-diversity-array/infinium-global-diversity-array-8-v1-0-D1-demo-data-36.zip | 465 MB; 3 CEU trios incl. NA12878 + replicates | +| GDA v1 sample sheet | https://support.illumina.com/content/dam/illumina-support/documents/downloads/productfiles/global-diversity-array/infinium-global-diversity-array-8-v1-0_D1_SampleSheet_Demo_36.csv | has Parent1/Parent2/Replicate | +| GDA v1 BPM / CSV / EGT | https://webdata.illumina.com/downloads/productfiles/global-diversity-array/infinium-global-diversity-array-8-v1-0-D1-manifest-file-bpm.zip (and `…-csv.zip`, `…-D1-cluster-file.zip`; page: https://support.illumina.com/array/array_kits/infinium-global-diversity-array/product-files.html) | | +| GSA v4 PGx demo (36 samples) | https://support.illumina.com/downloads/infinium-gsa-v4-pgx-demo-data-set-iscan.html | 872 MB; product files at https://support.illumina.com/array/array_kits/infinium-global-screening-array-v4-pgx/product-files.html | +| GRCh37 & GRCh38 FASTA | e.g. GRCh38 no-alt analysis set (NCBI) | contigs named per manifest | + +### 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 + +```markdown +# 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.md` fully. We are starting **Phase P0 (Setup)** then **P1 (Parsers)**. +> 1. 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. +> 2. Write `scripts/bootstrap_refs.sh` (shallow clones of §17 repos into `refs/`) and `scripts/fetch_data.sh` +> (tiers p0/p1/p2, sha256 recording, resumable downloads with aria2c/wget). +> 3. Run bootstrap + fetch p0. +> 4. Implement `formats/binutil.py` and `formats/idat.py` per §6.1 with tests (Bioconductor fixtures, malformed-file +> tests, property tests), then `simurg 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). diff --git a/scripts/bootstrap_refs.sh b/scripts/bootstrap_refs.sh new file mode 100755 index 0000000..509569c --- /dev/null +++ b/scripts/bootstrap_refs.sh @@ -0,0 +1,24 @@ +#!/usr/bin/env bash +# Shallow-clone reference open-source repositories into refs/ (git-ignored). +# GPL repos are for use as external oracles only — do not copy their code into SimurgArray. +set -euo pipefail +cd "$(dirname "$0")/.." +mkdir -p refs +repos=( + "Illumina/BeadArrayFiles" # BSD-2 + "Illumina/GTCtoVCF" # Apache-2.0 + "freeseek/gtc2vcf" # MIT + "freeseek/mocha" # MIT + "broadinstitute/picard" # MIT (see src/main/java/picard/arrays) + "WGLab/PennCNV" # Unlicense + "sbslee/pypgx" # MIT + "PharmGKB/PharmCAT" # MPL-2.0 (oracle) + "HenrikBengtsson/illuminaio" # GPL-2 (oracle only) + "gjun/verifyIDintensity" # GPL-3 (oracle only) +) +for r in "${repos[@]}"; do + d="refs/${r#*/}" + if [[ -d "$d/.git" ]]; then echo "skip $d"; continue; fi + git clone --depth 1 "https://github.com/$r.git" "$d" +done +echo "done: $(ls refs | wc -l) repos in refs/" diff --git a/scripts/fetch_data.sh b/scripts/fetch_data.sh new file mode 100755 index 0000000..c94b31e --- /dev/null +++ b/scripts/fetch_data.sh @@ -0,0 +1,65 @@ +#!/usr/bin/env bash +# Download SimurgArray datasets into data// and record sha256 checksums. +# Vendor files are for personal research use only; never commit them. +# Usage: bash scripts/fetch_data.sh p0|p1|p2 +set -euo pipefail +cd "$(dirname "$0")/.." +tier="${1:?usage: fetch_data.sh p0|p1|p2}" +mkdir -p "data/$tier" +sums="data/CHECKSUMS.sha256"; touch "$sums" + +get() { # get + local url="$1" dir="$2" f + f="$dir/$(basename "${url%%\?*}")" + if [[ -s "$f" ]]; then echo "have $f"; else + if command -v aria2c >/dev/null; then aria2c -x8 -s8 -c -d "$dir" "$url" + else wget -c -P "$dir" "$url"; fi + fi + grep -q " $f\$" "$sums" || sha256sum "$f" >> "$sums" +} + +GDA=https://webdata.illumina.com/downloads/productfiles/global-diversity-array +G1K=https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/supporting/hd_genotype_chip + +case "$tier" in +p0) + d=data/p0/gda_demo36; mkdir -p "$d" + get "$GDA/infinium-global-diversity-array-8-v1-0-D1-demo-data-36.zip" "$d" + get "https://support.illumina.com/content/dam/illumina-support/documents/downloads/productfiles/global-diversity-array/infinium-global-diversity-array-8-v1-0_D1_SampleSheet_Demo_36.csv" "$d" + get "$GDA/infinium-global-diversity-array-8-v1-0-D1-manifest-file-bpm.zip" "$d" + # CSV manifest and EGT: file names can change; if these fail, download manually from + # https://support.illumina.com/array/array_kits/infinium-global-diversity-array/product-files.html + get "$GDA/infinium-global-diversity-array-8-v1-0-D1-manifest-file-csv.zip" "$d" || echo "!! fetch CSV manifest manually" + get "$GDA/infinium-global-diversity-array-8-v1-0-D1-cluster-file.zip" "$d" || echo "!! fetch EGT manually" + if command -v Rscript >/dev/null; then + Rscript -e 'if(!requireNamespace("BiocManager",quietly=TRUE)) install.packages("BiocManager",repos="https://cloud.r-project.org"); BiocManager::install(c("IlluminaDataTestFiles","hapmap370k"),ask=FALSE,update=FALSE)' + Rscript -e 'cat(system.file("extdata",package="IlluminaDataTestFiles"),"\n")' > data/p0/bioc_fixture_path.txt + else + echo "!! R not found: install IlluminaDataTestFiles + hapmap370k from Bioconductor for parser fixtures" + fi + echo "NOTE: GSA v4 PGx demo (36 samples) + product files: download manually from" + echo " https://support.illumina.com/downloads/infinium-gsa-v4-pgx-demo-data-set-iscan.html" + echo " https://support.illumina.com/array/array_kits/infinium-global-screening-array-v4-pgx/product-files.html" + ;; +p1) + d=data/p1/1kg_omni_broad; mkdir -p "$d" + for f in HumanOmni2.5-4v1-Multi_B.bpm HumanOmni2.5-4v1_B.csv; do get "$G1K/broad_intensities/$f" "$d" || true; done + echo "Listing broad_intensities/ (pick EGT + the 107 GB Omni25_idats_gtcs_2141_samples.tgz):" + curl -s "$G1K/broad_intensities/" | grep -o 'href="[^"]*"' | cut -d'"' -f2 | grep -v '^?' || true + echo "Then: bash scripts/fetch_data.sh p1-big (downloads the 107 GB archive)" + d=data/p1/1kg_omni_sanger; mkdir -p "$d" + echo "Sanger set index:"; curl -s "$G1K/sanger_intensities/" | grep -o 'href="[^"]*"' | cut -d'"' -f2 | grep -v '^?' || true + d=data/p1/gse165845; mkdir -p "$d" + get "https://ftp.ncbi.nlm.nih.gov/geo/series/GSE165nnn/GSE165845/suppl/GSE165845_RAW.tar" "$d" + ;; +p1-big) + get "$G1K/broad_intensities/Omni25_idats_gtcs_2141_samples.tgz" data/p1/1kg_omni_broad + ;; +p2) + d=data/p2; mkdir -p "$d/gse96909" "$d/gda_pgx480" + get "https://ftp.ncbi.nlm.nih.gov/geo/series/GSE96nnn/GSE96909/suppl/GSE96909_RAW.tar" "$d/gse96909" + get "https://s3.amazonaws.com/webdata.illumina.com/downloads/productfiles/global-diversity-array-pgx/GDAPGx_Demo_idats.zip" "$d/gda_pgx480" + ;; +*) echo "unknown tier $tier"; exit 2;; +esac +echo "checksums in $sums"