Add sanitized Phase-0 research archive under docs/research
Research notes R1-R4, documentation-study files A-E (source of the functional spec's `src:` tags) and PGx coverage CSVs, rewritten to use neutral naming and current project decisions (plan §2).
This commit is contained in:
kemal committed
1 parent
180b6f160a
commit
91e97424c7
12 files changed
+7748
No files matched your search
@@ -0,0 +1,608 @@
|
||||
# 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
|
||||
Reference in new issue
Block a user