kemal
d11226ca31
- vendor comparison (1000 Genomes Omni2.5, AutoCall 1.6.2.2): every computed field identical in 2138/2141 samples; the 3 others trace back to the single-locus raw-intensity rule (K-009); sex rule 1 matches all - AutoCall 1.6.2.x leaves GC10/GC50/#calls/#no-calls at 0 in GTC v3: new `--preset autocall-1.6.2` (1.1.2 normalization, sex rule 1, GTC v3, legacy statistics); the sample summary keeps the real statistics - NaN sign canonicalised to the x86 pattern (K-014); report P4.md
56 lines
3.1 KiB
Markdown
56 lines
3.1 KiB
Markdown
# Phase P4 validation report — genotype calling and GTC output
|
||
|
||
Date: 2026-09-28. Development and oracle tests on the Pi; vendor comparison on a self-destructing AWS spot instance
|
||
(`run-20260927T233126Z`, 44 min, ≈ $0.19). Summary: `P4_omni2141_summary.json`.
|
||
|
||
## Delivered
|
||
|
||
| Item | Module |
|
||
|---|---|
|
||
| Cluster file aligned to manifest order (float32 mean R) | `genotype/clusters.py` |
|
||
| GenCall-equivalent calling: SD flair, six-case half-flattened likelihoods, priors/posteriors, fuzzy sub-scores, score map | `genotype/gencall.py` |
|
||
| BAF / LRR / LogR dev (Peiffer 2006), canonical NaN pattern | `genotype/baf_lrr.py` |
|
||
| Sex estimate, rules v1 and v2 | `genotype/sex.py` |
|
||
| IDAT pair → complete GTC (v5 or v3), TOP-strand base calls, sample stats, presets | `genotype/sample.py` |
|
||
| `simurg genotype` (GTC per sample, `sample_summary.tsv/json` per OUT-060…063, run manifest) | `pipeline/genotype.py` |
|
||
|
||
## Acceptance (plan §18 P4)
|
||
|
||
### Against the oracle (gtc2vcf built without FMA), P0 GSA data
|
||
|
||
| Check | Result |
|
||
|---|---|
|
||
| 36 samples × normalization 1.2.0 and 1.1.2: whole GTC files after aligning the four free-text fields (file names, AutoCall date/version) | **72/72 byte-identical** |
|
||
|
||
Every data field matches bit for bit: raw X/Y, transforms, genotypes, GenCall scores, TOP base calls, BAF, LRR, call
|
||
rate, sex, LogR dev, GC10/GC50, call counts, percentiles and controls. The only adjustment is the sign of NaNs that
|
||
libm computes (K-014).
|
||
|
||
### Against 2141 vendor GTCs (1000 Genomes Omni2.5, AutoCall 1.6.2.2, cluster file `HumanOmni2.5-4v1-Multi_B.egt`)
|
||
|
||
| Field | Samples that differ | Explanation |
|
||
|---|---|---|
|
||
| genotypes | 3 (5 loci in total) | the 3 samples of K-009, whose raw intensity differs in one locus; that difference shifts their normalization |
|
||
| GenCall scores | 3 (18 229 loci) | same 3 samples (the transform change moves every score of the affected bin) |
|
||
| normalization transforms | 2 | same K-009 samples |
|
||
| raw X / raw Y | 1 / 2 | K-009 |
|
||
| call rate, LogR dev | 2 / 3 | follow from the above |
|
||
| base calls | 3 | follow from genotypes |
|
||
| **sex (rule 1)** | **0 of 2141** | rule 2 would differ in 12 samples, so AutoCall 1.6.2.2 uses rule 1 |
|
||
| GC10, GC50, #calls, #no-calls | 2141 | AutoCall 1.6.2.2 leaves these at 0 in GTC v3 (K-015); now emulated by `--preset autocall-1.6.2` |
|
||
|
||
So in 2138 of 2141 samples every computed field agrees with the vendor; the remaining 3 samples trace back to the
|
||
single-locus raw-intensity rule of K-009. Whole-file byte identity was 0/2141 in this run only because of the four
|
||
legacy-empty statistics; the preset that reproduces them was added after the run (unit-tested, not yet re-run on
|
||
the cloud — expected 2138/2141).
|
||
|
||
Not covered here: BAF and LRR against vendor output (GTC v3 has no BAF/LRR fields). The 1000 Genomes Sanger set
|
||
ships a GenomeStudio Final Report with BAF/LRR for 177 samples — optional follow-up run (≈ $0.10).
|
||
|
||
## Performance
|
||
|
||
| Machine | Result |
|
||
|---|---|
|
||
| Raspberry Pi 5, GSA (683 k loci) | 3.6 s per sample (1 core); 36 samples in 2 min with 4 workers |
|
||
| AWS m7g.4xlarge, Omni2.5 (2.45 M loci), 16 workers | 2141 samples genotyped and compared in 19.6 min; peak RSS 3.2 GB |
|