kemal/SimurgArray
1
0 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
85229fba84d4776ae9c2a0f84cced3c7db6c2d2a
SimurgArray/docs/research/R3_science_formats_algorithms.md
T
kemal 5695119e9c 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).
2026-09-27 16:27:34 +01:00

46 KiB
Raw Blame History

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

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)

Reference in new issue View Git Blame Copy Permalink