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.