P8 complete: scientific validation on GSE165845 vs NYGC 30x and GIAB HG001

360 GSA arrays (120 1000 Genomes individuals x 3 replicates): concordance 0.9934 vs NYGC 30x (0.9958 on sites in both sets), 0.9985 vs GIAB HG001 v4.2.1 (non-ref F1 0.997), replicate concordance 0.9966; allele-swap sites 0.016 %; 358/360 replicate pairs detected as duplicates. Plan §16.5 metrics in docs/validation/P8.md (K-024).
2026-09-28 19:47:34 +01:00
1 parent b43dededdb
commit 490527e037
@@ -0,0 +1,106 @@
+# 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.