# R3 — Science, File Formats and Algorithms behind Illumina Infinium Genotyping (for SimurgArray)

> **Archived Phase-0 study (sanitized copy, see `README.md`).** Written 2026-09-27, before the decisions in
> `SIMURGARRAY_MASTER_PLAN.md` §2. The analysed vendor product is called "the reference product" and its CLI
> `refcli`; its documentation URLs are replaced by `[ref-docs:…]` tokens. Where this file speaks of a clean room,
> legal clearance, or running/benchmarking the vendor binary, plan §2 supersedes it: personal research; MIT/BSD/Apache
> code may be reused with attribution; GPL tools are subprocess oracles only; vendor binaries are never run for
> comparison, inspected or decompiled; oracles are published vendor-produced GTCs, gtc2vcf and BeadArrayFiles.

Author role: statistical-genetics research pass for an independent re-implementation (SimurgArray) of reference-product-like functionality.
Date of research: 2026-09-27.

## Evidence tags

| Tag | Meaning |
|---|---|
| **[PUBLIC DOCUMENTATION]** | Illumina technical notes, support bulletins, format specs, reference-product online help, Illumina patents |
| **[SCIENTIFIC LITERATURE]** | Peer-reviewed papers (DOIs given; bibliographic data from PubMed) |
| **[OPEN-SOURCE IMPLEMENTATION]** | Behavior read from public source code (Illumina/BeadArrayFiles, illuminaio, freeseek/gtc2vcf + idat2gtc, BAFregress) |
| **[INFERRED]** | My deduction from combining the sources above |
| **[HYPOTHESIS]** | Plausible but unverified; needs testing against real Illumina output |

### Clean-room and licensing notes (read first)

* **Illumina/BeadArrayFiles** (GTC/BPM/EGT parsers, GTC v5 format PDF) is published by Illumina. It is the primary format reference used here. https://github.com/Illumina/BeadArrayFiles
* **illuminaio** (R/Bioconductor) is GPL-licensed. Its encrypted-IDAT path uses a DES routine (`src/des.c`, GPL-2+) with a key embedded in the package. Do not copy code. Re-derive only the structural facts. https://github.com/HenrikBengtsson/illuminaio
* **freeseek/gtc2vcf `idat2gtc.c`** is MIT-licensed. Its header says it "reimplements functionalities and ideas present in AutoConvert (v1.6.3.1), GTCtoVCF, BeadArrayFiles". It also cites Illumina's released *GenCall Source Code* (C#: `NormalizationInfinium.cs`, `GenTrain60.cs`, `Utils.cs`; http://support.illumina.com/downloads/gencall_software.html). It deliberately reproduces Illumina bugs "for interoperability". **This is the most complete public description of AutoCall/GenCall behavior. Its provenance includes Illumina binaries and Illumina-released source, so get a legal opinion before porting it line by line.** The safest clean-room use is: (a) describe the algorithms from the patents, papers and tech notes, and (b) use idat2gtc and the Illumina binaries only as black-box **test oracles**. The key patents US7035740 (normalization) and US7467117 (clustering/calling) are marked **expired** on Google Patents. [PUBLIC DOCUMENTATION] https://patents.google.com/patent/US7467117B2/en ; https://patents.google.com/patent/US7035740

---

## 1. Infinium chemistry

### 1.1 Assay principle
* Whole-genome amplification, fragmentation, and hybridisation to 50-mer probes on beads (bead type = "address"). Then single-base extension (SBE) with labelled terminators, staining, and two-colour imaging on iScan. [SCIENTIFIC LITERATURE] Gunderson et al. 2005, Nat Genet 37:549, https://doi.org/10.1038/ng1547 ; Steemers et al. 2006, Nat Methods 3:31, https://doi.org/10.1038/nmeth842
* Each bead type is present in about 15–30 copies. IDATs store per-address summary statistics (mean, SD, bead count), not per-bead data. [PUBLIC DOCUMENTATION]/[OPEN-SOURCE IMPLEMENTATION] (IDAT fields below)

### 1.2 Infinium I vs Infinium II
* **Infinium II** (one bead type per locus, `AddressA_ID` only, `AddressB_ID`=0; BPM assay_type 0).
  * The probe ends one base 5′ of the SNP. The SNP base itself is extended with a labelled ddNTP.
  * A and T terminators carry one hapten, detected in the **red** channel. C and G carry the other, detected in the **green** channel.
  * So the two alleles appear in different colours. [SCIENTIFIC LITERATURE] Steemers 2006.
  * Consequence: [A/T] and [C/G] SNPs cannot be assayed with Infinium II, because both alleles have the same colour. [INFERRED from chemistry; consistent with manifests]
* **Infinium I** (two bead types per locus, `AddressA_ID` and `AddressB_ID`; assay_type 1 or 2).
  * The two probes are allele-specific at their 3′ end. Extension adds the *next* (non-polymorphic) base.
  * Both alleles therefore read in **one** channel, determined by the next base: A/T → red (assay_type 1), C/G → green (assay_type 2).
  * [OPEN-SOURCE IMPLEMENTATION] `get_assay_type()` in idat2gtc.c infers assay type from the bases flanking `[x/y]` in SourceSeq for old manifests (locus-entry version < 6); newer BPMs store it explicitly.
* **Intensity assembly (IDAT → raw X/Y per locus).** [OPEN-SOURCE IMPLEMENTATION] idat2gtc `fill_array` (matches AutoCall `SampleData.fillArray`):
  * Infinium II: X = Red.Mean[AddressA], Y = Grn.Mean[AddressA]
  * Infinium I A/T: X = Red.Mean[AddressA], Y = Red.Mean[AddressB]
  * Infinium I C/G: X = Grn.Mean[AddressA], Y = Grn.Mean[AddressB]
  * A locus is set to (0,0) unless **both** addresses have NBeads ≥ 2.
  * Controls: controls_x = Red, controls_y = Grn for each control address listed in the BPM control config.
  * **Discrepancy:** the GTC v5 PDF describes Raw X as "green" and Raw Y as "red". The open implementations use X = red for Infinium II. Treat the PDF wording as unreliable and validate against real GTCs. [INFERRED]

### 1.3 A/B alleles, TOP/BOT, PLUS/MINUS, Forward/Reverse
Source: Illumina technical note "TOP/BOT Strand and A/B Allele". [PUBLIC DOCUMENTATION] https://www.illumina.com/documents/products/technotes/technote_topbot.pdf

**TOP/BOT algorithm:**
1. **Unambiguous SNPs.**
   * If one allele is A and the other is C or G → strand = **TOP**, Allele A = A, Allele B = C/G.
   * If one allele is T and the other is C or G → strand = **BOT**, Allele A = T, Allele B = C/G.
2. **Ambiguous SNPs ([A/T], [C/G]) — "sequence walking".**
   * Let the SNP be position n. Compare the pairs (n−1, n+1), (n−2, n+2), … until the first *unambiguous* pair, meaning one base is A or T and the other is C or G.
   * If the A or T of that pair lies **5′** of the SNP → TOP. If it lies **3′** → BOT.
   * Alleles: [A/T] on TOP gives A=A, B=T; on BOT gives A=T, B=A. [C/G] on TOP gives A=C, B=G; on BOT gives A=G, B=C.
3. dbSNP adopted TOP/BOT in 2005 (per the same note).

**Other strand designations** [PUBLIC DOCUMENTATION] https://knowledge.illumina.com/microarray/general/microarray-general-reference_material-list/000001489 :
* **PLUS/MINUS**: relative to the reference FASTA. The + strand's 5′ end is at the tip of the p arm.
* **Forward/Reverse** (`SourceStrand`): Forward means "the alleles in the SNP column match the RefSNP alleles displayed in dbSNP".
* **IlmnStrand**: the design strand. It is TOP/BOT for SNPs and PLUS/MINUS for indels.
* **RefStrand**: the +/− strand of the SNP column. It is only present in newer manifests (BPM locus-entry version 8).

**Base calls in GTC** are always reported on the TOP strand. [PUBLIC DOCUMENTATION] GTC v5 spec. idat2gtc derives them from the SNP column `[a/b]` and IlmnStrand: TOP → (a, b); BOT → complement(a), complement(b). [OPEN-SOURCE IMPLEMENTATION]
* PLUS and Forward calls: BeadArrayFiles `get_base_calls_plus_strand` / `_forward_strand` take the design-strand SNP `[a/b]`. They complement each allele when the strand annotation (RefStrand/SourceStrand) differs from the requested report strand. A '-' is emitted for NC/NULL/unknown strand. [OPEN-SOURCE IMPLEMENTATION] GenotypeCalls.py

**Indels:**
* The SNP column is `[I/D]` or `[D/I]`. IlmnStrand is PLUS/MINUS.
* The complement of I is I and of D is D (BeadArrayUtility.complement).
* Mapping I/D to actual REF/ALT sequence needs the CSV manifest (SourceSeq / TopGenomicSeq). gtc2vcf states "Indel alleles cannot be inferred and will be skipped without the CSV manifest file." [OPEN-SOURCE IMPLEMENTATION] https://github.com/freeseek/gtc2vcf
* idat2gtc also normalises a SourceSeq of the form `x[-/y]` by swapping the dash position. [OPEN-SOURCE IMPLEMENTATION]

---

## 2. File formats

General conventions: little-endian. Strings are .NET `BinaryWriter` strings, i.e. a **7-bit varint (LEB128) length prefix** followed by UTF-8 bytes (BeadArrayUtility.read_string). The GTC PDF says the first byte is the length; that is only true for strings shorter than 128 bytes. [OPEN-SOURCE IMPLEMENTATION]

### 2.1 IDAT (binary, version 3)
[OPEN-SOURCE IMPLEMENTATION] illuminaio `readIDAT_nonenc.R`; idat2gtc.c (field-name comments such as `ID_N_CORES` appear to come from Illumina/Broad naming).

**Header:**
```
char[4]  magic = "IDAT"
int64    version (must be 3 for the unencrypted binary format)
int32    nFields
nFields × { uint16 fieldCode; int64 byteOffset }   # TOC
```

**Field codes:**

| code | name (illuminaio / idat2gtc) | type at offset |
|---|---|---|
| 1000 | nSNPsRead / ID_N_CORES | int32 N (lowest offset in file) |
| 102 | IlluminaID / ID_BEAD_TYPES | int32[N] (address IDs) |
| 103 | SD / ID_DEVS | uint16[N] |
| 104 | Mean / ID_MEANS (a trimmed mean in practice) | uint16[N] |
| 107 | NBeads / ID_N_GOOD_BEADS | uint8[N] |
| 200 | MidBlock / ID_ILLUMICODES | int32 count + array (illuminaio reads int32[]; idat2gtc reads uint8[] — **discrepancy**, verify) |
| 300 | RunInfo / ID_PROCESS_HISTORY | int32 n; n × 5 strings (RunTime, BlockType, BlockPars, BlockCode, CodeVersion) |
| 400 | RedGreen / ID_TENTH_PERCENTILE | 4 bytes |
| 401 | Manifest / ID_SAMPLE_BEADSET | string |
| 402 | Barcode (Sentrix ID) | string |
| 403 | ChipType / ID_SENTRIX_FORMAT | string |
| 404 | SentrixPosition / ID_SECTION_LABEL (e.g. R01C01) | string |
| 405 | Beadset / "opa" | string |
| 406 | SampleID / ID_DNA | string |
| 407 | Description / "descr" | string |
| 408 | Plate | string |
| 409 | Well | string |
| 410 | SampleCount | int32 (illuminaio treats as unknown bytes) |
| 510 | VLN / "unknown" | string |

Codes 100, 101, 105, 106, 108 (backgrounds, background devs, medians, n_beads, trimmed means) are named in idat2gtc comments but not used.

* One IDAT exists per colour per section: `<barcode>_<R##C##>_Grn.idat` and `..._Red.idat`.
* The IlluminaID arrays of the two files index the same addresses. Join by address, not by position (idat2gtc builds a hash address→index). [OPEN-SOURCE IMPLEMENTATION]

**Encrypted IDAT ("version 1", older BeadScan/GoldenGate era)** [OPEN-SOURCE IMPLEMENTATION] illuminaio `readIDAT_enc.R` + `des.c`:
* The file is DES-encrypted. Decrypted, it is XML carrying base64-encoded binary blocks named ` __<Name>="..."`.
* The data blocks are six float32 arrays followed by four int32 arrays, plus `<ProcessEntry>` run info and `<BarCode>`, `<SentrixFormat>`, `<SectionLabel>`.
* The DES key is embedded in illuminaio (GPL). SimurgArray should treat encrypted IDATs as **out of scope** or support them only via a user-supplied external converter. [INFERRED]

### 2.2 BPM (binary bead pool manifest)
[OPEN-SOURCE IMPLEMENTATION] BeadArrayFiles `BeadPoolManifest.py`; idat2gtc `bpm_init`/`locusentry_read`.
```
char[3] "BPM"; uint8 = 1
int32 version  (may carry flag 0x1000 -> XOR it off; supported 3..5)
string manifest_name
string control_config            (if version>1; CSV-like lines "addr[:addr...],Category,Color,Name")
int32  num_loci
int32[num_loci] indexes          (1-based order index of each locus entry)
string[num_loci] names           (sorted order)
uint8[num_loci] norm_ids         (0..99, sub-bead-pool / normalization bin)
num_loci × LocusEntry
int32 m_header; string[m_header] header lines   (idat2gtc reads this trailer; BeadArrayFiles ignores it)
```

**LocusEntry** (int32 version 4, 6, 7 or 8):
* Strings, in order: IlmnID, Name, 3 unused strings, **int32 index**, unused string, IlmnStrand, SNP, Chr, Ploidy, Species, MapInfo (string), TopGenomicSeq (v4 only), CustomerStrand.
* Then int32 AddressA, int32 AddressB.
* Then strings: AlleleA_ProbeSeq, AlleleB_ProbeSeq (v4 only), GenomeBuild, Source, SourceVersion, SourceStrand, SourceSeq (v4 only).
* **v≥6:** 1 byte unused, int8 Exp_Clusters, int8 Intensity_Only, uint8 AssayType (0 = Inf II, 1 = Inf I A/T, 2 = Inf I G/C). Consistency rule: AddressB==0 ⇔ AssayType==0.
* **v≥7:** 4 × float (frac_A, frac_C, frac_T, frac_G).
* **v8:** string RefStrand.

**Normalization-bin lookup** (critical for exact reproduction):
```
norm_id'  = (norm_id + 100 * assay_type) mod 256      # BeadArrayFiles comment: mimics byte wrap in GenomeStudio/AutoCall/IAAP
lookup    = rank of norm_id' among the sorted set of distinct norm_id' values
XForm used for locus i = GTC.normalization_transforms[lookup[i]]
```
idat2gtc notes that Omni5 has overflow cases where some Infinium I G/C probes then share a bin with Infinium II probes. [OPEN-SOURCE IMPLEMENTATION]

### 2.3 CSV manifest (human-readable)
[PUBLIC DOCUMENTATION] https://knowledge.illumina.com/microarray/general/microarray-general-reference_material-list/000001565
* 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.
* Column meanings:
  * Exp_Clusters: 1 = non-polymorphic, 2 = mtDNA/Y, 3 = other diploid.
  * Intensity_Only: 1 = "out of stats"; always no-call; used for CN.
  * TopGenomicSeq: CSV only, not in BPM.
* The CSV (not the BPM) is needed for indel allele resolution and for VCF conversion in the reference product (`gtc-to-vcf`). [PUBLIC DOCUMENTATION] [ref-docs:input-files]

### 2.4 EGT cluster file
[OPEN-SOURCE IMPLEMENTATION] BeadArrayFiles `ClusterFile.py`; idat2gtc `egt_init`; Picard `InfiniumEGTFile.java` is cited there.
```
int32 version (=3)
string gencall_version, cluster_version, call_version, normalization_version, date_created
uint8 is_wgt (must be 1)
string manifest_name
int32 data_block_version  (5..9 seen)
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
   if dbv >= 7 (idat2gtc) / ==9 (BeadArrayFiles): float intensity_threshold + 14 unused floats
   elif dbv >= 6: 15 unused floats ; elif dbv < 5: 2 unused floats
num_records × ClusterScore: float cluster_separation, float total_score (= GenTrain score), float original_score, uint8 edited
num_records × string genotype   (unused)
num_records × string locus name
num_records × int32 address
num_records × int32[3] cluster counts (must equal N_AA, N_AB, N_BB)
if dbv == 9: num_records × float (unknown; skipped by idat2gtc)
```
* EGT records are matched to BPM loci **by name**.
* The derived quantity R̄ = (N_AA·R_AA + N_AB·R_AB + N_BB·R_BB)/(N_AA+N_AB+N_BB) is used as the expected R for Intensity_Only loci. [OPEN-SOURCE IMPLEMENTATION]
* The unused 14–15 floats per record are **undocumented** (possibly cluster covariances or GenTrain internals). [HYPOTHESIS]

### 2.5 GTC genotype file (versions 3, 4, 5)
[PUBLIC DOCUMENTATION] Illumina, `docs/GTC_File_Format_v5.pdf` in BeadArrayFiles; [OPEN-SOURCE IMPLEMENTATION] GenotypeCalls.py.

**Header:** `char[3] "gtc"; uint8 version; int32 M;` then M TOC entries `{int16 id; uint32 offset}`.
* For int-typed IDs (1, 2, 3) the "offset" **is the value**.
* Arrays are `int32 count` followed by the elements.

**TOC IDs:**

| ID | Name | Type |
|---|---|---|
| 1 | NumSNPs | int32 (in TOC) |
| 2 | Ploidy | int32 (in TOC) |
| 3 | PloidyType (1 diploid, 2 autopolyploid, 3 allopolyploid) | int32 (in TOC) |
| 10/11/12 | Sample name / plate / well | string |
| 100 | Cluster file name | string |
| 101 | SNP 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[] ((sections/sample) × nControls) |
| 1000/1001 | Raw X / Raw Y | uint16[] |
| 1002 | Genotypes | uint8[] |
| 1003 | BaseCalls | char[2][] (TOP strand; '-' = no call) |
| 1004 | GenCall scores | float[] |
| 1005 | Scanner data: name, pmt_green (int32), pmt_red (int32), scanner version, imaging user | struct |
| 1006 | Call rate | float |
| 1007 | Gender | char 'M'/'F'/'U' |
| 1008 | LogR dev | float |
| 1009 | p10GC | float |
| 1010 | DX flag | int32 |
| 1011 | Extended sample data: P50GC float, NumCalls int32, NumNoCalls int32, NumIntensityOnly int32 | struct |
| 1012 | B allele freqs | float[] (v≥4) |
| 1013 | LogR ratios | float[] (v≥4) |
| 1014/1015 | Intensity percentiles X/Y: P05, P50, P95 | 3 × uint16 (v5) |
| 1016 | Sentrix ID | string (v5) |

* **Genotype codes:** 0 NC, 1 AA, 2 AB, 3 BB, 4 NULL, 5 A, 6 B, 7 AAA … up to 45 BBBBBBBB (polyploid series).
* **XForm (52 B):** int32 version; float offset_x, offset_y, scale_x, scale_y, shear, theta; then 6 reserved floats. idat2gtc writes nn12, rr12, taa, tbb into these reserved slots.
* **Write-complete check:** v3 → NumIntensityOnly readable; v4 → LogR ratios readable; v5 → Sentrix ID readable.
* **The GTC does not store normalized X/Y, theta or R.** They are recomputed from raw X/Y + XForm + BPM norm lookup.

**Normalized intensities (exact, from the spec):**
```
tx  = x_raw − offset_x ;           ty  = y_raw − offset_y
tx2 =  cosθ·tx + sinθ·ty ;         ty2 = −sinθ·tx + cosθ·ty
tx3 = tx2 − shear·ty2 ;            ty3 = ty2
xn  = tx3 / scale_x ;              yn  = ty3 / scale_y
```
* Negative results are clamped to 0 (BeadArrayFiles `threshold=True`; idat2gtc clamps before scaling).
* If x_raw ≤ 0 and y_raw ≤ 0, the result is NaN.

### 2.6 Illumina (GenomeStudio) Final Report
[PUBLIC DOCUMENTATION] GenomeStudio Genotyping Module v2 User Guide (11319113), https://support.illumina.com/content/dam/illumina-support/documents/documentation/software_documentation/genomestudio/genomestudio-2-0/genomestudio-genotyping-module-v2-user-guide-11319113-01.pdf ; [OPEN-SOURCE IMPLEMENTATION] BeadArrayFiles `examples/gtc_final_report.py`.
* Tab-delimited text.
* `[Header]` block: `GSGT Version`, `Processing Date`, `Content` (manifest), `Num SNPs`, `Total SNPs`, `Num Samples`, `Total Samples`.
* `[Data]` block, "Standard" (long) format: one row per SNP × sample. Columns are user-selectable, e.g. `SNP Name, Sample ID, Allele1 - Top, Allele2 - Top, Allele1 - Forward, Allele1 - AB, Allele1 - Plus, GC Score, Chr, Position, Theta, R, X, Y, X Raw, Y Raw, B Allele Freq, Log R Ratio, GT Score, Cluster Sep`.
* A "Matrix" format also exists (SNP × sample).
* Exact column naming and ordering vary with the selected fields. [INFERRED]

### 2.7 Locus/DNA summary reports
[OPEN-SOURCE IMPLEMENTATION] `examples/locus_summary.py`, `dna_summary.py`.
* Locus columns: Call_Freq, A/A_Freq, A/B_Freq, B/B_Freq, Minor_Freq, Gentrain_Score, 50%_GC_Score, 10%_GC_Score, Het_Excess_Freq, ChiTest_P100, Cluster_Sep, AA/AB/BB T_Mean, T_Std, R_Mean, R_Std, Plus/Minus.
* Formulas are in §8.

---

## 3. Normalization

### 3.1 Published description
* [PUBLIC DOCUMENTATION] US7035740 (Kermani, "Artificial intelligence and global normalization methods for genotyping") and US7467117 (divisional). The global affine model has "6 degrees of freedom — Xoffset, Yoffset, Xscale, Yscale, Xshear, and Yshear". Control points are the data points nearest to "sweep points" along the axes, found via Delaunay triangulation. Transform order: translation, then rotation, then shear, then scale. https://patents.google.com/patent/US7467117B2/en
* [SCIENTIFIC LITERATURE] Peiffer et al. 2006 (https://doi.org/10.1101/gr.5402306) list five steps:
  1. outlier removal;
  2. translation — asymptotes fitted to candidate AA and BB homozygotes; their intersection is the new origin;
  3. rotation — the angle of the AA asymptote to the x-axis;
  4. shear — the angle of the BB asymptote to the rotated y-axis;
  5. scaling — robust centroids of candidate AA (x) and BB (y), giving a scale of about 1.
* [PUBLIC DOCUMENTATION] Illumina, "Illumina's Genotyping Data Normalization Methods", Pub. 970-2006-010 (mirror: http://dnatech.genomecenter.ucdavis.edu/wp-content/uploads/2013/06/illumina_gt_normalization.pdf — the fetch failed with ECONNREFUSED; content known via quotes in idat2gtc):
  * **400** sweep points.
  * A SNP is an outlier if x, y or x/(x+y) is below min(5th smallest, 1st percentile) or above max(5th largest, 99th percentile).
* [PUBLIC DOCUMENTATION] GenTrain 3.0 tech note 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):
  * Infinium I and II loci go in separate normalization bins.
  * Bins with **< 192 loci** degrade normalization.
  * GenTrain 3.0 fits a **2-parameter** model for Infinium I bins with < 192 loci, because both alleles are in one colour.
  * GenTrain 3.0 ignores zero-intensity samples when estimating the BB mean in small bins.

### 3.2 Full algorithm as reproduced
[OPEN-SOURCE IMPLEMENTATION] idat2gtc; normalization versions 1.1.2 = AutoConvert/GenTrain 2, and 1.2.0 = AutoConvert 2.0 / IAAP / ACLI / GenTrain 3.

For each bin b (loci grouped by norm_id', raw uint16 X,Y as float; missing pairs enter as zeros — a reproduced quirk):

**GenTrain-3 small Infinium I bin** (norm_id' ≥ 100 and n < 192) uses the single-channel 2-parameter model:
1. Mirror: swap so that x ≥ y.
2. c = P2(y) (2nd percentile). Subtract c from x and y.
3. Polar coordinates with r = |x|+|y| and t = (2/π)·atan2(y,x).
4. AA set = { t < 0.1 }.
5. ȳ = trimmean20(y_AA) − 50.
6. r̄ = trimmean20(r_nonNaN, zero-padded) − 2ȳ.
7. Result: offset_x = offset_y = c + ȳ, θ = 0, shear = 0, scale_x = scale_y = r̄.

**Otherwise, six-parameter fit:**
1. **Outliers:** skipped if n < 400. Remove points with x ≤ xcut1, x ≥ xcut2, y ≤ …, or t = y/(x+y) outside (tcut1, tcut2), where cut1 = min(5th smallest, P1) and cut2 = max(5th largest, P99).
2. **Offset** (skipped if n < 192, giving 0):
   * Shift by (min x, min y).
   * Build 400 equally spaced points on the x-axis, [0, max x], and 400 on the y-axis.
   * Candidate AA = unique nearest data points (Euclidean) to the x-axis points. Candidate BB = the same for the y-axis points.
   * Robust line fits: y = b_a + m_a·x on the AA set; x = b_b + m_b·y on the BB set.
   * Intersection: o_x = (b_b + m_b·b_a)/(1 − m_b·m_a), o_y = b_a + m_a·o_x.
   * offset = o + (min x, min y).
3. **nn12/rr12** (diagnostics stored in reserved slots):
   * nn12 = fraction of points with x < 0.2·P95(x) and y < 0.2·P95(y).
   * rr12 = mean(x+y over those points) / mean(x+y over all).
4. **Rotation:**
   * New 400-point x-sweep over [min x, max x]; AA set = union with the previous set.
   * Robust through-origin fit y = m·x (bisquare).
   * θ = atan(m). Rotate: x′ = cosθ·x + sinθ·y, y′ = −sinθ·x + cosθ·y.
5. **Shear:**
   * y-sweep [0, max y]; BB set = union.
   * Robust fit x = m·y through the origin. shear = tan(atan m) = m. Apply x′ = x − shear·y.
6. **Scale:**
   * n ≥ 192: AA = { x > 5y and x > 0.1·P99(x) }, BB = { y > 5x and y > 0.1·P99(y) }. scale_x = trimmean50(x_AA), scale_y = trimmean50(y_BB).
   * n < 192: clamp negatives to 0. AA = { t° < 10° }, BB = { t° > 80° } using degree angles (GenTrain 2 also counts zero points as AA). scale = trimmean20.
   * Base variant used in the older code path: scale = 0.5·trimmean50 + 0.7413·IQR.

**Robust fit details:**
* MATLAB `robustfit` equivalent: bisquare weights, tuning constant 4.685, at most 50 iterations, relative tolerance 1e-6, leverage-adjusted residuals.
* Scale estimate: `madsigma` = median(|r| sorted, from the p-th element)/0.6745. GenTrain 3 uses the new madsigma; GenTrain 2 the old median-centred version.
* The no-intercept fit uses a constant adjfactor ≈ 100 (a reproduced Illumina bug).

**Percentile definition (used everywhere, also for GC10/GC50):** for sorted v of length n and target p:
```
i1 = floor(n·p/100); f = n·p/100 − i1; if f < 0.5: i1 −= 1
clamp to ends; x1 = 100(i1+0.5)/n; x2 = 100(i1+1.5)/n
P = v[i1] + (v[i1+1]−v[i1])·(p − x1)/(x2 − x1)
```
Trimmed mean trimmeanK = mean of values v with P(K/2) ≤ v ≤ P(100−K/2).

### 3.3 Theta and R
* [PUBLIC DOCUMENTATION] GenTrain2 tech note (Pub. 037-2009-015, https://www.illumina.com/documents/products/technotes/technote_gentrain2.pdf): θ = (2/π)·arctan(B/A), R = A + B.
* [OPEN-SOURCE IMPLEMENTATION] Edge cases:
  * θ = (2/π)·atan2(yn, xn), R = xn + yn.
  * (0,0) gives NaN for both.
  * idat2gtc uses R = |xn|+|yn|, and if both are negative sets R = FLT_MIN·FLT_EPSILON.

---

## 4. GenCall / GenTrain

### 4.1 What Illumina has published
[PUBLIC DOCUMENTATION] Tech notes: GenCall (Pub. 370-2004-009, https://www.illumina.com/Documents/products/technotes/technote_gencall_data_analysis_software.pdf), GenTrain2 (037-2009-015), GenTrain 3.0 (M-GL-01258).

**GenTrain clustering** (builds the EGT from many samples):
1. **Preliminary clustering** on θ, with no fixed number of clusters (1 to about 6).
2. **Secondary clustering:** candidate models assign preliminary clusters to AA/AB/BB (all 3, 2, or 1 cluster present). Each cluster gets mean and SD of θ and R, so 12 parameters per locus.
3. **Final clustering:** each model is scored on compactness, spread between clusters, and the likelihood of the assignment under HWE. The highest-scoring model is kept.
* Missing clusters are placed by **neural networks** (patent: 3-layer perceptron, trained with a genetic algorithm; 7 cluster models for a bi-allelic locus).
* The GenTrain score comes from a fuzzy/ANN score "evolved using top/bottom strand data", combined with penalty terms.
* GenTrain 3.0 caps the negative penalty for homozygous clusters lying on the axis.

**GenCall calling:**
* A Bayesian classification against the EGT clusters.
* GenCall score = GenTrain score × data-to-model fit score, then "further penalized by the DNA Score".
* It is "not a probability". The sensitive range is 0.2–0.7; scores < 0.2 usually mean failed genotypes.
* **No-call threshold 0.15** (FastTrack standard, chosen from trio and replicate data). [PUBLIC DOCUMENTATION] technote_infinium_genotyping_data_analysis.pdf
* Call rate = #genotypes with score ≥ threshold / #genotypes.

### 4.2 Calling algorithm in full
[OPEN-SOURCE IMPLEMENTATION] idat2gtc `compute_score_call_prelim`, `gen_std_flair`, `modilik`, `gencall_score_map`, `make_calls`. These cite GenTrain60.cs / GenTrain62. The README claims output "equivalent to the output of the Illumina GenCall algorithm".

Given a locus with EGT cluster k ∈ {AA, AB, BB} = (μθ_k, σθ_k, μR_k, σR_k, N_k), GenTrain score G, intensity threshold τ, and sample point (t, r):

1. If r < τ → NC, score NaN (then 0).
2. **Pre-flooring:** σR_k ← max(σR_k, 0.1).
3. **"gen_std_flair" dispersion adjustment:**
   * z1 = ½(μθ_AB−μθ_AA)/(σθ_AB+σθ_AA), z2 = ½(μθ_BB−μθ_AB)/(σθ_BB+σθ_AB), mz = min(z1, z2).
   * η = max(mz/3, 1)·min(mz/3, 1). Set σθ_k ← max(η·σθ_k, 0.02).
   * m = max(0.2, median(σθ_AA, σθ_AB, σθ_BB) of the original values). Set σR_k ← max(σR_k, max(μR_k/7, m)).
4. **Likelihoods ("modilik"):**
   * d_k = 100·|t−μθ_k| + |r−μR_k|.
   * Choose the two nearest clusters (by the d ordering and R ordering). For the pair, use **half-flattened Gaussians** on the outward sides: `normpdf_vleft(x;μ,σ)` equals the peak density for x<μ, and `vright` equals the peak for x>μ.
   * Each L_k = φ*(t;μθ_k,σθ_k)·φ*(r;μR_k,σR_k), with the flattened side chosen by the relative positions of clusters in θ and R. See the source for the six cases.
   * The third cluster keeps L=0. If no case applies, use plain Gaussians for all three.
5. **Priors:** P_k = (N_k/N + ω)/(1+3ω), with ω = 1.
6. **Posteriors:** AP_k = L_k P_k / Σ L_j P_j. Genotype = argmax.
7. **Sub-scores:**
   * s_ap = smf(log10(AP_max/AP_2nd); 0, 2), where smf is MATLAB's S-shaped function.
   * 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 is present and N_AB = 0; 0.95 if 2 clusters are present; else 1.
   * s_n = 1.
8. prelim = s_ap · G · s_θ · s_r · s_n · s_mis. NaN if G is NaN or t is NaN.
9. **Final GenCall score:** GC = prelim^0.35 · gbellmf(prelim; a=0.504, b=1.71, c=1.08), where gbellmf(x) = 1/(1+|(x−c)/a|^{2b}). NaN becomes 0.
10. Genotype ← NC if GC < cutoff (default **0.15**). Intensity_Only loci are forced to NC and counted separately.
11. Call rate = NumCalls / (NumSNPs − NumIntensityOnly).

**Standard MATLAB fuzzy memberships used** (restated):
```
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 = 1 − smf
```

### 4.3 Published vs unknown for GenCall/GenTrain

| Item | Status |
|---|---|
| θ/R transform, 12-parameter cluster model | Published (tech note) |
| 0.15 no-call threshold; GC10/GC50 definitions | Published |
| Per-sample calling given an EGT | Reproduced in open source (idat2gtc), MIT, patent expired; "equivalent" per its author [OPEN-SOURCE IMPLEMENTATION]; needs our own concordance test |
| "DNA score" penalty mentioned in the 2005 tech note | **Unknown / apparently absent** in modern AutoCall (s_n = 1 in idat2gtc) [INFERRED] |
| GenTrain clustering (preliminary clustering method, candidate-model scoring, HWE term, ANN weights for missing clusters) | **Not public beyond qualitative descriptions + patent**. GenTrain score and Cluster Sep formulas are **unknown**. |
| Meaning of the extra 14–15 floats per EGT record | Unknown |
| The reference product's genotyper = ACLI GenCall (norm 1.2.0/GenTrain 3)? | [HYPOTHESIS] strongly suggested because the reference product consumes the same EGT and BPM. Validate by comparing GTCs. |

**SimurgArray implication** [INFERRED]:
* Calling from a *given* EGT can be implemented to near-bit-exact fidelity.
* Producing new EGTs (reclustering) cannot match GenTrain. Use a documented alternative, e.g. a GMM/t-mixture in (θ,R) with HWE-aware model selection, and compute our own separation score. Label it as non-GenTrain.

---

## 5. BAF and LRR

[SCIENTIFIC LITERATURE] Peiffer et al. 2006 (https://doi.org/10.1101/gr.5402306): LRR = log2(R_obs/R_exp). R_exp is interpolated from the canonical cluster positions of about 120 reference samples. BAF is the θ position relative to the cluster θ means.
[OPEN-SOURCE IMPLEMENTATION] idat2gtc `get_baf_lrr` (matching AutoCall `CalculateBAlleleFreq`/`CalculateLogRDev`).

With cluster θ means (θ_AA, θ_AB, θ_BB) and R means (R_AA, R_AB, R_BB) from the EGT, and sample (θ, R):
```
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 (θ > θ_AB):     BAF = 0.5 + 0.5·(θ − θ_AB)/(θ_BB − θ_AB)
                     R_exp = R_AB + (θ − θ_AB)·(R_BB − R_AB)/(θ_BB − θ_AB)
(θ NaN → BAF = LRR = NaN)
Intensity_Only loci: R_exp = N-weighted mean R of the three clusters (R̄ above)
LRR = log2(R / R_exp)     (R == 0 → −FLT_MAX)
BAF stored in GTC clamped to [0,1]; θ outside [θ_AA, θ_BB] extrapolates linearly before clamping
```
**LogR dev** = population SD of LRR over loci that are autosomal (chromosome not starting with X, Y or M), not Intensity_Only, and finite: sqrt(E[LRR²] − E[LRR]²). Missing or NaN EGT cluster means propagate NaN.

**Notes:**
* The EGT sometimes lacks one cluster (monomorphic loci). GenomeStudio then uses ANN-predicted positions stored in the EGT. [INFERRED]
* gtc2vcf `--adjust-clusters` re-centres the clusters on the study data to reduce LRR noise. MoChA also applies median BAF/LRR adjustments. [OPEN-SOURCE IMPLEMENTATION] https://github.com/freeseek/gtc2vcf
* cnvPartition expected BAF for copy number n with b B-alleles is b/n; LRR means: CN0 −5 (SD 2), CN1 −0.45, CN2 0, CN3 0.3, CN4 0.75 (SD 0.18); BAF SD 0.03. [PUBLIC DOCUMENTATION] technote_cnv_algorithms.pdf

---

## 6. Alternative genotype callers

| Method | Model | Key facts |
|---|---|---|
| **Illuminus** — Teo et al. 2007, Bioinformatics 23:2741, https://doi.org/10.1093/bioinformatics/btm443 | [SCIENTIFIC LITERATURE] Model-based, population-level, no training data. Per-SNP mixture of bivariate (truncated) **t-distributions** in transformed contrast/strength space, fitted by EM; optimal initialisation handles shifted clouds; perturbation analysis gives a stability-based quality metric. | Needs large n: poorer at ≤50 samples (Ritchie 2011). Has a sex-aware X model. |
| **GenoSNP** — Giannoulatou et al. 2008, Bioinformatics 24:2209, https://doi.org/10.1093/bioinformatics/btn386 | [SCIENTIFIC LITERATURE] **Within-sample**: clusters all SNPs of one sample (per bead pool) with a 4-component Student-t mixture (AA, AB, BB, null) on log intensities, fitted by **variational Bayes**. | >99% accuracy on HapMap; no reference population; works for n=1. Also useful as independent QC. |
| **CRLMM (Illumina)** — Carvalho et al. 2007, Biostatistics 8:485, https://doi.org/10.1093/biostatistics/kxl042 ; Ritchie et al. 2009, Bioinformatics 25:2621, https://doi.org/10.1093/bioinformatics/btp470 | [SCIENTIFIC LITERATURE] Preprocessing: IDAT summaries, then quantile normalisation, then log-ratio M = log2(A/B). Hierarchical two-level mixture with SNP-specific cluster centres shrunk towards HapMap-derived priors, fitted by EM; per-call confidence; SNR sample metric. | Bioconductor `crlmm`. |
| **Birdseed** (Birdsuite) — Korn et al. 2008, Nat Genet 40:1253, https://doi.org/10.1038/ng.237 | [SCIENTIFIC LITERATURE] Affymetrix SNP6. 2-D Gaussian mixture EM with SNP-specific priors; integrates copy-number-aware genotypes (A-null, AAB …). | Affy-only, but the CN-aware genotype idea transfers. |
| **zCall** — Goldstein et al. 2012, Bioinformatics 28:2543, https://doi.org/10.1093/bioinformatics/bts479 | [SCIENTIFIC LITERATURE] Post-processor for GenCall no-calls at rare variants. Uses the mean and SD of the common homozygote cluster in normalized X and Y to set z-score thresholds (T_x = μ_x + z·σ_x, likewise for y) that define the other two genotype regions. | Improved rare-allele detection vs exome sequencing on HumanExome. |
| **optiCall** — Shah et al. 2012, Bioinformatics 28:1598, https://doi.org/10.1093/bioinformatics/bts180 | [SCIENTIFIC LITERATURE] Combines SNP-wise and sample-wise calling (mixture-model based) for rare, low-frequency and common variants. | On Immunochip (n=4537), SNPs lost to false QC failure or misclassified as monomorphic: 1.38% (optiCall) vs 3.87% (Illuminus), 7.85% (GenoSNP), 4.09% (GenCall). |

**Head-to-head comparison:** Ritchie et al. 2011, BMC Bioinformatics 12:68, https://doi.org/10.1186/1471-2105-12-68 [SCIENTIFIC LITERATURE]
* GenCall, Illuminus, GenoSNP and CRLMM differ only slightly in accuracy.
* CRLMM and GenoSNP consistently beat GenCall, especially at low MAF.
* Sex-aware models (Illuminus, CRLMM) are better on chrX.
* Sample QC flags agree well across the four methods.

---

## 7. CNV, LOH and mosaicism

* **PennCNV** — Wang et al. 2007, Genome Res 17:1665, https://doi.org/10.1101/gr.6861907 [SCIENTIFIC LITERATURE]
  * First-order HMM over markers ordered by position. Hidden states: CN0, CN1, CN2, CN2-LOH (copy-neutral LOH), CN3, CN4.
  * Emissions combine LRR (a state-specific Gaussian mixed with a small uniform component for outliers) and BAF. BAF emission uses the population **B allele frequency (PFB)** to weight genotype-specific BAF components.
  * Transitions depend on inter-marker distance d, roughly P(change) ∝ 1 − exp(−d/D).
  * Supports trio/joint calling.
  * GC-wave correction: Diskin et al. 2008, NAR 36:e126, https://doi.org/10.1093/nar/gkn556 — regress LRR on local GC content; the GC-wave factor (GCWF) tracks DNA quantity.
* **QuantiSNP** — Colella et al. 2007, NAR 35:2013, https://doi.org/10.1093/nar/gkm076 [SCIENTIFIC LITERATURE]
  * Objective-Bayes HMM on LRR+BAF. Hyperparameters are calibrated to a fixed type-I error by resampling.
  * Reports Bayes factors per segment.
* **cnvPartition** (Illumina GenomeStudio plug-in) [PUBLIC DOCUMENTATION] https://www.illumina.com/documents/products/technotes/technote_cnv_algorithms.pdf
  * **Per-locus likelihoods** for 14 genotypes as bivariate Gaussians: L_g = N(LRR; μ_g, σ_g)·N(BAF; β_g, 0.03), with μ, σ from the table in §5. For DD: L_DD = N(LRR; −5, 2) with BAF uniform.
  * **Copy-number likelihoods:** L_k = Σ_{g∈CN k} L_g. Preliminary CN: X = (L1+2L2+3L3+4L4)/(L0+…+L4).
  * **Segmentation:** maximise |Z_ij| with Z_ij = [1/(j−i) + 1/(n−j+i)]^{−1/2} · [(S_j−S_i)/(j−i) − (S_n−S_j+S_i)/(n−j+i)], where S = cumulative sum of X.
    * Sliding windows w = 4, 8, 16, 32, then extended in both directions.
    * Threshold |Z| > 6; recurse on the pieces.
  * **Region CN** = argmax_k Σ log L_k. **Confidence** = Σ log L_k − Σ log L_2 (report threshold 35). Probe-gap filter 1 Mb.
  * Filters: CN1 on X/Y in males is removed; Y is ignored in females; BAF is ignored for Intensity_Only probes.
  * **CN-LOH:** homozygosity detector on CN2 regions, MinHomozygousRegionSize 10 Mb.
  * **Homozygosity Detector:** E_het = Σ 2f_i(1−f_i) and χ² = N·E_het/(N−E_het), with cut-off 23.5 and ≥ 50 SNPs.
  * **LOH Score** = log10 Π P(gt|LOH)/P(gt|noLOH), with P_error = 0.001.
* **MoChA** — Loh et al. 2018, Nature 559:350, https://doi.org/10.1038/s41586-018-0321-x ; Loh et al. 2020, Nature 584:136, https://doi.org/10.1038/s41586-020-2430-6 [SCIENTIFIC LITERATURE]
  * Detects **mosaic** chromosomal alterations using *phased* BAF: at heterozygous sites, allelic imbalance appears as |BAF−0.5| with a consistent haplotype sign.
  * An HMM over phased hets tolerates phase switch errors; LRR classifies loss / gain / CN-LOH.
  * Sensitive to cell fractions of a few percent. Open-source: https://github.com/freeseek/mocha [OPEN-SOURCE IMPLEMENTATION]
* **Reference-product CNV** [PUBLIC DOCUMENTATION] [ref-docs:home]
  * Inputs: a "Cytogenetics model file (.dat) … contains probe GC content information"; a "PGx CN model file (.dat)" paired to the manifest; a mask file (.msk) for CN training.
  * PGx CN analysis needs ≥ 24 samples, with ≥ 22 passing LogR dev < 0.2.
  * Output VCF QUAL is Phred-scaled and capped at 60.
  * The **model internals are not published.** [HYPOTHESIS] GC regression of LRR plus a per-probe reference (baseline) model, then HMM/segmentation.

---

## 8. QC metrics

**Per sample** [OPEN-SOURCE IMPLEMENTATION] / [PUBLIC DOCUMENTATION]:
* **Call rate** = NumCalls/(NumSNPs − NumIntensityOnly). Expectation > 0.99 for human products (tech note).
* **p10GC / p50GC:** the Illumina percentile (§3.2) of GenCall scores over loci with score > cutoff and not Intensity_Only. Sample GC10 < ~0.2 flags a poor DNA (GenCall tech note).
* **LogR dev:** SD of autosomal LRR (§5). Reference-product PGx CNV QC uses < 0.2.
* **Intensity percentiles:** P05/P50/P95 of raw X and raw Y over all loci (GTC v5, uint16).
* **Sex estimation** [OPEN-SOURCE IMPLEMENTATION] idat2gtc, describing the AutoCall versions. Only loci with GenTrain score ≠ 0 count.
  * *v1* (AutoConvert 1.6): with more than 20 called X loci, male if X-het-rate ≤ 0.1, else female.
  * *v2* (AutoConvert 2.0 / IAAP / ACLI):
    * Requires autosomal call rate > 0.97, estimated on up to 10,000 random autosomal loci, and ≥ 20 Y loci / ≥ 100 autosomal loci.
    * Then: median R of Y loci > 0.3 → M. Else if X loci < 20 → F. Else median R of X loci > 0.9 → F, otherwise U.
    * IAAP 1.1 had a bug counting null-cluster loci as missing.
  * Output 'M', 'F' or 'U'.
* **Sex from genotypes** (literature): X heterozygosity / PLINK F on X (F > 0.8 male, < 0.2 female), plus Y intensity. Laurie et al. 2010 combine X and Y intensities with X-het to separate sex mix-ups from aneuploidy. [SCIENTIFIC LITERATURE] https://doi.org/10.1002/gepi.20516
* **Heterozygosity rate / inbreeding F:** flag samples beyond ±3 SD of the mean het rate. Anderson et al. 2010, Nat Protoc 5:1564, https://doi.org/10.1038/nprot.2010.116 [SCIENTIFIC LITERATURE]
* **Contamination:**
  * **VerifyIDintensity** — Jun et al. 2012, AJHG 91:839, https://doi.org/10.1016/j.ajhg.2012.09.004 [SCIENTIFIC LITERATURE].
    * Models the observed allele intensity ratio as a mixture of (1−α) of the true sample and α of an unknown contaminant drawn from population allele frequencies.
    * Per-genotype Gaussian intensity components; α is estimated by maximum likelihood over pairs of genotypes weighted by HWE priors.
    * Detects contamination down to about 1%.
  * **BAFRegress** (same group; http://genome.sph.umich.edu/wiki/BAFRegress; reimplemented in gtc2vcf `BAFregress.c`) [OPEN-SOURCE IMPLEMENTATION]:
    * Use sites with population B-allele frequency p ∈ [0.1, 0.9].
    * AA-called sites: regress y = BAF − median_AA(BAF) on x = p.
    * BB-called sites: regress y = 1 − (BAF − (median_BB(BAF) − 1)) on x = 1 − p.
    * OLS slope m = (nΣxy − ΣxΣy)/(nΣx² − (Σx)²) is the contamination estimate (E[BAF|AA] ≈ α·p). The median adjustment needs ≥ 5 genotypes per cluster.
* **Relatedness / identity:** KING robust kinship (Manichaikul et al. 2010, https://doi.org/10.1093/bioinformatics/btq559); PLINK IBD π̂. [SCIENTIFIC LITERATURE]

**Per locus** [OPEN-SOURCE IMPLEMENTATION] locus_summary.py, reproducing GenomeStudio:
* Call_Freq = calls/(calls+nocalls). Genotype frequencies are over calls.
* Minor_Freq = min(p_A, 1−p_A), where p_A = (2·AA+AB)/(2·calls).
* **Het_Excess:**
  * q = Minor_Freq, p = 1−q, k = (0.013/q)⁴.
  * dh = (AB/calls + k)/(2pq + k) − 1.
  * ChiTest_P100 = 2Φ(dh; 0, 0.1) if dh < 0, else 2(1 − Φ(dh; 0, 0.1)).
  * If there are no calls → (0,0). If monomorphic → (1,0).
* GC10/GC50 per locus use the same percentile rule.

**Locus filter thresholds** (hard cut-off / grey zone) [PUBLIC DOCUMENTATION] technote_infinium_genotyping_data_analysis.pdf, Table 1:
* Call_Freq < 0.97 / ≤ 0.99
* Cluster Sep ≤ 0.3 / < 0.45
* AA/AB/BB R Mean ≤ 0.2 / < 0.4
* 10%GC ≤ 0.3 (grey)
* Het Excess > 0.2 (grey; the text also says < −0.3)
* AB T Mean < 0.2 or > 0.8 (grey)
* AA T Mean > 0.3 when A/A_Freq = 1, etc.
* **Cluster Sep definition:** "separation between the three genotype clusters in the theta dimension", range 0–1. The exact formula is **unpublished**.

---

## 9. Pharmacogenomics on arrays

* **Reference-product PGx** [PUBLIC DOCUMENTATION] [ref-docs:home] ; [ref-docs:dev-article-1.0-pgx]
  * Star-allele calling from array genotypes, with definitions from **PharmVar/PharmGKB**. Phenotype is looked up in **CPIC or DPWG**.
  * Star-allele-specific copy numbers; CYP2D6 call rate > 97%.
  * CYP2D6 CN is measured at sub-gene regions (e.g. CYP2D6.exon.9, CYP2D6.intron.2, CYP2D6.p5), which enables detection of hybrids and deletions/duplications.
  * Inputs: PGx database .zip (variant mapping), PGx CN model .dat, mask .msk.
  * Supported products: GDA-ePGx, GSAv4-ePGx, GCRA-ePGx.
  * The internal diplotype-scoring and CN algorithms are **not published**.
* **Star-allele nomenclature:** PharmVar — Gaedigk et al. 2018, CPT 103:399, https://doi.org/10.1002/cpt.910 [SCIENTIFIC LITERATURE]
* **Diplotype inference from genotypes:** PharmCAT — Sangkuhl et al. 2020, CPT 107:203, https://doi.org/10.1002/cpt.1568. It matches observed alleles against allele-definition tables and scores candidate diplotypes by the number of matched defining positions. [SCIENTIFIC LITERATURE]
* **CYP2D6 SV/CN** (NGS read-depth paradigm, adaptable to array LRR): Stargazer — Lee et al. 2019, Genet Med 21:361, https://doi.org/10.1038/s41436-018-0054-0. Uses paralog-specific CN (CYP2D6 vs CYP2D7) and handles hybrids. [SCIENTIFIC LITERATURE]
* **Array-based CYP2D6 CNV caller + PGx passport pipeline:** Asterix — Lanting et al. 2025, AJHG 112:1015, https://doi.org/10.1016/j.ajhg.2025.03.003. Clinically validated calling, phasing and imputation for 12 genes plus a CYP2D6 CNV tool on biobank arrays. [SCIENTIFIC LITERATURE]
* **Suggested SimurgArray approach** [HYPOTHESIS]:
  1. Per-probe LRR normalised against a batch reference, restricted to probes in CYP2D6 sub-regions and flanking CYP2D7/REP regions, with a GC regression.
  2. Classify CN per sub-region with a Gaussian mixture on the median LRR. Train on ≥ 24 samples, mirroring the reference product's minimum.
  3. Combine region CNs into structural haplotypes (*5 deletion, xN duplications, *36/*13 hybrids).
  4. Assign SNV-defined star alleles PharmCAT-style, with CN-consistent BAF checks (e.g. BAF ≈ 1/3 or 2/3 for a duplication).

---

## 10. Summary: specified vs unknown

**Fully specified publicly** (can be implemented and verified bit-for-bit or near):
* IDAT v3, BPM (locus entry v4–v8), EGT (v3, data blocks 5–9) and GTC v3–5 binary layouts.
* The normalization transform application, θ/R, BAF/LRR interpolation and LogR dev.
* Genotype codes, TOP/BOT/A-B rules, base-call strand conversion.
* Normalization-bin lookup including the mod-256 wrap.
* Sample stats (call rate, GC10/GC50 percentile rule, intensity percentiles).
* Locus summary statistics including Het Excess.
* cnvPartition likelihood model and segmentation statistic.
* The 0.15 no-call threshold.

**Reproduced only in third-party open source** (idat2gtc, MIT; derived from Illumina binaries and source):
* The normalization *fit* (robust fits, sweeps, trims, 192/400 constants).
* The GenCall Bayesian call and score (flair, modilik, sub-scores, gbell map).
* Sex-estimation rules.
* Use it as a test oracle; legal review is recommended before direct porting. *(Superseded: plan §2 allows MIT reuse with attribution.)*

**Unknown / proprietary:**
* GenTrain clustering internals (preliminary clustering, model scoring weights, ANN weights for missing clusters).
* GenTrain score and Cluster Sep formulas.
* The EGT's unused per-record floats.
* The "DNA score" penalty.
* The reference product's CNV and PGx-CN models (.dat) and star-allele diplotype scoring.
* Whether the reference product's genotyping is byte-identical to ACLI [HYPOTHESIS: yes, GenTrain 3 / normalization 1.2.0].

## Key URLs (formats)
* GTC spec PDF: https://github.com/Illumina/BeadArrayFiles/blob/develop/docs/GTC_File_Format_v5.pdf
* BeadArrayFiles modules: https://github.com/Illumina/BeadArrayFiles/tree/develop/module
* illuminaio IDAT readers: https://github.com/HenrikBengtsson/illuminaio/tree/master/R
* gtc2vcf / idat2gtc / BAFregress: https://github.com/freeseek/gtc2vcf
* Public Illumina manifests and EGTs (FTP list): https://github.com/freeseek/gtc2vcf/blob/master/Illumina.md