SimurgArray/docs/validation/V2_arraytrain.md

V2 — ArrayTrain (plan §12)

simurg train builds an EGT cluster file from a cohort's GTCs. It does not reproduce the vendor's clustering, so the benchmark is genotype accuracy against sequencing truth: an EGT trained by SimurgArray versus the vendor's published cluster file, used by the same GenCall emulation on the same arrays.

Set-up

• Data. GSE165845: 120 1000 Genomes individuals, each on 3 GSA v1 arrays. The C2 build-38 BPM is used; the public GSA v1 C1 cluster file is the vendor baseline (K-024).
• Training. Replicate 3, 120 GTCs (59 F, 60 M, 1 U), all 618 540 loci. Pi 5, 69 min, peak 1.9 GB.
• Evaluation. Replicates 1–2 (240 arrays) go from IDAT through simurg genotype with the trained EGT, then simurg gtc-to-vcf. The same arrays are also genotyped with the vendor EGT.
• Truth. The NYGC 30x phased 1000G panel on chr22: 8278 array SNV sites (scripts/train/concordance_chr.py). Concordance counts called genotypes only. NRD is the non-reference discordance.

The evaluation individuals are also in the training set, but on different arrays. The cluster file never sees their truth; the same holds for the vendor file.

Results (240 arrays, chr22)

Cluster file	Call rate	Concordance	NRD
ArrayTrain v1 (fixed θ windows)	99.47 %	99.27 %	1.89 %
ArrayTrain v2 (BIC + HWE labelling)	99.61 %	99.44 %	1.30 %
Vendor GSA v1 C1 EGT	99.15 %	99.46 %	1.23 %

V2_arraytrain_gse165845.json holds the numbers. v2 calls 0.46 % more genotypes than the vendor file at nearly the same concordance (−0.02 points).

What changed from v1, and why

The v1 training fitted one three-component θ mixture. Its means were bounded to fixed windows (AA ≤ 0.4, 0.2 ≤ AB ≤ 0.8, BB ≥ 0.6), and close clusters were merged by a fixed θ gap.

A per-locus diagnosis on chr22 compared our trained clusters, the vendor clusters and the truth for 80 arrays:

• Errors concentrate in few loci. The top 20 loci held 43 % of the errors.
• 2-cluster loci were the problem. There, v1 made 1418 errors against the vendor's 186.
• The cause was the fixed θ windows. Many loci are compressed or asymmetric in θ; for example the vendor file has AA 0.05 / AB 0.27 / BB 0.45, or AA 0.37 / AB 0.68 / BB 0.96. With fixed windows such loci were mislabelled: rs2003814 (truth 14 AB + 66 BB) was clustered as AA + AB.

v2 (src/simurg/train/clusters.py) replaces this with:

1. Unbounded 1-, 2- and 3-component mixtures on θ, plus a uniform outlier component. Each is fitted by EM from several starts: canonical positions and θ quantiles.

   • An SD floor of 0.02 inside the EM stops a component collapsing onto the θ = 0 / 1 spikes.
   • BIC chooses among admissible fits. A fit is admissible if every component has ≥ max(3, 1 %) of the samples and adjacent components are ≥ 1.5 standardised gaps apart.

2. Labelling of the θ-ordered components by the best non-decreasing map onto AA/AB/BB. A map's score adds up:

   • the HWE goodness of fit of the implied genotype counts, with an inbreeding allowance F = 0.1. This is the fit relative to the saturated multinomial, so monomorphic labellings are not favoured automatically;
   • weak priors on the cluster positions: AA ~ N(0.03, 0.2), AB ~ N(0.5, 0.25), BB ~ N(0.97, 0.2);
   • a weak prior on the θ gaps between genotype groups: 0.42 ± 0.12 per step.

   Adjacent components may share a genotype (a split cluster) only when their standardised gap is < 3.5.

Fast chr22 loop (train on replicate 3, call 80 evaluation arrays, score against truth):

Step	Errors (80 arrays)
v1	2562
Unbounded mixtures + HWE labelling of 2-cluster fits	1793
+ MIN_SEP 1.5 and position prior	1611
+ general labelling with split-cluster merges (HWE goodness of fit)	1380
Vendor EGT	1237

• Parameter sweep. MERGE_SEP {2.5, 3.5, 5} × HWE F {0.05, 0.1, 0.2} × gap SD {0.12, 0.2}. Every setting lands at 99.70–99.79 % accuracy; MERGE_SEP 5 is clearly worse. The chosen values are in the middle of a flat optimum.
• Where the errors are now. By our cluster count, v2's errors are 123 at 1-cluster loci (vendor 157), 217 at 2-cluster loci (vendor 113) and 1040 at 3-cluster loci (vendor 967).
• Many 3-cluster errors are shared with the vendor. Both files see three clear clusters where the truth is monomorphic, e.g. rs191067607 and rs5758036. These are probe or truth problems, not clustering errors.

Limits (K-030)

• Only one cohort and one chromosome were used for accuracy. The tuning used this evaluation, so the numbers are optimistic for other chips and cohorts.
• 120 training arrays is close to the 100-sample minimum. Rare-allele clusters are imputed from cross-locus relations, not observed.
• The trained EGT scores (GenTrain, cluster separation) are ours. GenCall scores from a trained file are not comparable with scores from vendor files. At a matched call rate, GenCall score ranks v1's errors worse than the vendor's.
• Training takes 69 min for 120 × 618 540 on the Pi (8 EM runs per locus); a Rust kernel would help at cohort scale.