SimurgArray/SIMURGARRAY_MASTER_PLAN.md

# 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: `<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.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 <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_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 `<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`:
```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)
| 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).