kemal
1eb7ea8154
Fixed θ windows mislabelled compressed/asymmetric loci (2-cluster loci: 1418 vs 186 errors). Now 1/2/3-component mixtures chosen by BIC, labelled by HWE goodness of fit plus weak position and gap priors, split clusters merged. GSE165845 chr22 vs 1000G 30x (240 arrays): concordance 99.44 % (v1 99.27 %, vendor EGT 99.46 %), NRD 1.30 % (1.89 %, 1.23 %), call rate 99.61 %. docs/validation/V2_arraytrain.md, K-030; plan §12, CLI-038.
87 lines
5.1 KiB
Markdown
87 lines
5.1 KiB
Markdown
# 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.
|