SimurgArray/docs/validation/P8.md

Phase P8 validation — scientific accuracy against whole-genome truth

Date: 2026-09-28. One AWS spot run (m7g.4xlarge, 16 vCPU, 36 min; plus an aborted first attempt, see Notes). Cost ≈ $0.33 in total. Full numbers: docs/validation/P8_gse165845_summary.json.

Data and method

Item	Source
Arrays	GEO GSE165845: 120 1000 Genomes individuals × 3 replicates = 360 samples (GSA-MD v1 BeadChips, 27 chips), IDATs from GEO
Manifest / clusters	GSA-24v1-0 C2 build 38 BPM + CSV, GSA v1 C1 cluster file (the only public GSA v1 EGT; 99.996 % of the manifest's bead addresses are on the GSA-MD chips)
Calls	simurg genotype (defaults, 360 samples in 1 min 55 s) → simurg gtc-to-vcf --merge --csv (cohort BCF, 4 min)
Truth 1	NYGC 1kGP 30x high-coverage phased panel (20220422, SNV/indel, GRCh38), all 120 individuals
Truth 2	GIAB HG001 (NA12878) v4.2.1 GRCh38, inside its confident-region BED, chromosomes 1–22
Reference	Ensembl GRCh38 release 110 primary assembly

Scripts: scripts/p8/prepare.py (GEO download, IDATs laid out by the barcode/position in their headers), scripts/cloud/tasks/p8.sh (the whole run), scripts/p8/validate_p8.py (metrics). Both call sets are split to biallelic records and left-normalised against the same FASTA (bcftools norm -m-any) and matched on (CHROM, POS, REF, ALT). Chromosomes 1–22 and X (GIAB: 1–22). An array site absent from the WGS call set is counted as homozygous reference (both call sets are complete for their individuals); a site where the WGS set has different alleles at the same position is excluded and counted (871 of 617 144). Genotype classes are hom-ref / het / hom-alt; a variant-level true positive needs the exact genotype. MAF comes from the NYGC AF over 3202 genomes.

Results (plan §16.5)

Against NYGC 30x — 360 samples, 617 144 array records, 219 M called genotypes

Metric	All sites	Sites in both sets	Excluding the one failed array (359 samples)
Genotype concordance	0.99338	0.99580	0.99478
Non-reference discordance (NRD)	0.0245	0.0150	—
No-call rate	0.0108	0.0110	0.0102
Non-ref sensitivity / precision / F1	0.9956 / 0.9768 / 0.9861	0.9956 / 0.9864 / 0.9910	—

Per genotype class (sites in both sets):

Class	Truth genotypes	Sensitivity	Precision	Specificity	F1
hom-ref	151.6 M	0.9959	0.9995	0.9986	0.9977
het	36.9 M	0.9943	0.9905	0.9980	0.9924
hom-alt	21.2 M	0.9979	0.9792	0.9976	0.9885

Per MAF bin (sites in both sets; non-ref precision is low for rare variants, the usual array behaviour — most "alt" calls of a very rare variant are false positives):

MAF	Genotypes	Concordance	NRD	Non-ref sensitivity	Non-ref precision	F1
< 0.1 %	3.9 M	0.9959	0.799	0.993	0.201	0.334
0.1–1 %	49.0 M	0.9975	0.152	0.997	0.849	0.917
1–5 %	38.9 M	0.9959	0.074	0.995	0.927	0.960
5–10 %	21.9 M	0.9947	0.030	0.995	0.972	0.983
10–20 %	32.8 M	0.9951	0.013	0.995	0.988	0.992
20–50 %	63.2 M	0.9952	0.008	0.996	0.994	0.995

SNVs: concordance 0.99336 (218.2 M genotypes). Indels: concordance 0.99754 (1.24 M genotypes), non-ref sensitivity 0.990, precision 0.828. Array sites absent from NYGC (27 510 records, 9.8 M genotypes): 94.2 % hom-ref calls; the remaining alt calls there are the main source of the lower "all sites" precision.

Per sample: median concordance 0.99485, median NRD 0.0189. Lowest: 0.3369 (the failed array below), then 0.9866.

Against GIAB HG001 v4.2.1 — NA12878 × 3 replicates, confident regions

Metric	Value
Array records inside the BED	577 617 (150 671 in the GIAB VCF, 426 658 hom-ref by absence)
Genotype concordance	0.99853 (replicates 0.99911, 0.99727, 0.99918)
NRD	0.0056
No-call rate	0.0165
Non-ref sensitivity / precision / F1	0.9992 / 0.9946 / 0.9969
hom-ref / het / hom-alt F1	0.9991 / 0.9962 / 0.9981
Indels (9371 genotypes)	concordance 0.99989, sensitivity 1.000, precision 0.991

Replicates, site errors, batch, relatedness

Metric	Result
Replicate concordance (360 pairs, both called)	pooled 0.99660; median 0.99963; without the failed array (358 pairs) 0.99942
Median non-ref discordance between replicates	0.0014
Systematically discordant sites (wrong in ≥ 50 % of ≥ 10 called samples)	2419 of 616 221 (0.39 %); 1444 of them absent from NYGC
Allele-swap pattern (strand/allele error: array = 2 − truth at ≥ 80 % of the discordant calls)	100 sites (0.016 %), 5 of them A/T or C/G
Per-chip batch (27 chips)	median call rate 0.9911, median LogR deviation 0.246; 4 chips flagged as call-rate outliers (3 × MAD): 203672270124 (holds the failed array; mean concordance 0.830), 203660990106 and 203660990145 (hold the 2nd and 3rd lowest-concordance arrays, chip means 0.992/0.993), 203660990061 (one array, call rate 0.980, concordance 0.994)
Replicates detected as duplicates by simurg genotype cohort relatedness	358/360 pairs; the 2 missed pairs involve the failed array; 0 duplicates between different individuals
Mendelian error rate	not measurable here (no complete trio among the 120 individuals); P5: median 0.031 % over 335 Omni2.5 trios (docs/validation/P5.md)

Lists: per_sample.tsv, replicate_pairs.tsv.gz, systematic_sites.tsv.gz in the run's results (s3://…/run-20260928T180818Z/p8/, kept out of the repository).

Notes

• One failed array. GSM5051309 (NA20362 replicate 2, 203672270124_R12C02) has call rate 0.762 and LogR deviation 1.38; simurg genotype flags it low_call_rate;high_lrr_sd;het_outlier. Its genotypes agree with its own replicates at only 0.336, so the sample itself (not the pipeline) is bad. The run used no QC profile, so the flag did not set a FAIL status; with --qc-profile genotyping it would have been excluded. The headline numbers keep it.
• Cluster file. GSE165845 was run on GSA-MD v1 chips, for which no public cluster file exists; the GSA v1 C1 EGT was used. The high median LogR deviation (0.246) reflects that mismatch, not intensity noise, and call rate and concordance could only improve with a matching EGT.
• Absent-as-hom-ref. 27 510 array records (4.5 %) have no NYGC record. Most are called hom-ref in all samples; the 1444 systematically "wrong" ones among them are array alt calls at sites NYGC does not report (filtered or truly monomorphic). The "sites in both sets" column shows the numbers without this assumption.
• First attempt aborted. The first run downloaded the 40 GB NYGC VCFs before using them; EBI served ≈ 15 MB/s, the CPU sat idle for 30 minutes and the idle alarm terminated the instance (as designed; ≈ $0.17). The task now streams each chromosome through bcftools view -T (manifest positions ± 100 bp, 120 individuals) while genotyping runs, which keeps the CPU busy and needs no disk for the truth set.
• Deferred parts of §16.5. CNV precision/recall and PGx diplotype concordance need the CN/LOH and PGx modules (plan §18 "v2 modules"); they are measured when those modules exist.