From f1e22ab33973c98e34d70ebc556c843df20083ef Mon Sep 17 00:00:00 2001 From: Kemal Yaylali Date: Sun, 27 Sep 2026 22:19:28 +0100 Subject: [PATCH] P2 validated on 2141 vendor GTCs (1000 Genomes Omni2.5) Raw X/Y and controls assembled from IDATs equal the vendor GTC arrays for 2138/2141 samples x 2.45M loci; 3 single-locus differences are recorded as K-009 (probable per-channel missing rule in AutoCall 1.6.2.2). Channel order X=Red confirmed (U1). `simurg ingest` processed all 2141 samples in 19 min (I/O-bound). AWS run cost ~$0.34 and left no resources. --- README.md | 5 +- docs/known_unknowns.md | 3 + docs/validation/P2.md | 39 +++++++ docs/validation/P2_omni2141_summary.json | 136 +++++++++++++++++++++++ 4 files changed, 181 insertions(+), 2 deletions(-) create mode 100644 docs/validation/P2.md create mode 100644 docs/validation/P2_omni2141_summary.json diff --git a/README.md b/README.md index c1b7008..57b25d5 100644 --- a/README.md +++ b/README.md @@ -9,8 +9,9 @@ Personal research toolkit for Illumina Infinium genotyping BeadChip data ## Status -Phases **P0 (setup)** and **P1 (file parsers)** are complete — see `docs/validation/P1.md`. -Available commands: `simurg version`, `simurg inspect [--json]`. +Phases **P0 (setup)**, **P1 (file parsers)** and **P2 (canonical model, raw X/Y)** are complete — see +`docs/validation/`. Validated on 2141 vendor GTCs of the 1000 Genomes Omni2.5 set. +Available commands: `simurg version`, `simurg inspect [--json]`, `simurg ingest`. ## Quick start diff --git a/docs/known_unknowns.md b/docs/known_unknowns.md index 47c55e5..853223b 100644 --- a/docs/known_unknowns.md +++ b/docs/known_unknowns.md @@ -4,6 +4,9 @@ Plan §20 lists the open questions; this file records what the implementation ha | ID | Date | Topic | Finding | Status | |---|---|---|---|---| +| K-011 | 2026-09-27 | Plan U1 (channel order) | Infinium II X = Red, Y = Grn confirmed on 2141 vendor GTCs (the GTC PDF's "Raw X = green" wording is wrong). | resolved | +| K-010 | 2026-09-27 | Ingest throughput | On AWS, `ingest` of 2141 Omni samples was I/O-bound (gp3 125 MB/s baseline): 18.9 min, 4 of 16 cores busy. For P9 benchmarks provision gp3 throughput or use instance-store NVMe. | open (P9) | +| K-009 | 2026-09-27 | Missing-intensity rule | 3 loci (in 3 of 2141 vendor GTCs, AutoCall 1.6.2.2) keep one channel's intensity where we zero both. Hypothesis: old AutoCall applies the NBeads ≥ 2 rule per channel. Confirm with the bead counts of 5767788121_R01C01 (rs178512), 5826354152_R02C01 (SNP12-29082770), 6268201099_R02C01 (rs17133206); if confirmed, add a `--missing-rule per-channel` compatibility option. | open | | K-008 | 2026-09-27 | Vendor GTC layout | 2141 vendor GTC v3 files (AutoCall 1.6.2.2) re-encode byte-identically with the canonical v3 order, so the writer reproduces the vendor layout (no padding or reordering). | resolved | | K-007 | 2026-09-27 | Demo data | The "GDA v1 Demo 36" download is a GenomeStudio project (`.bsc` + `Data/*.bin`), not IDATs. The GSA v4 PGx demo (`idats_Demo_36.zip`) is used as the P0 IDAT fixture instead. Bioconductor `hapmap370k` is no longer in the current release. | resolved (plan §17 updated) | | K-006 | 2026-09-27 | Sample sheet | Recognised input columns (`Sample_ID`, `Sample_Name`, …) collide case-insensitively with output columns `sample_id`/`sample_name` (IN-005). Recognised columns are exempt from the reserved-name check and are passed through verbatim. | decision, revisit in P4 (OUT-060 pass-through) | diff --git a/docs/validation/P2.md b/docs/validation/P2.md new file mode 100644 index 0000000..2394008 --- /dev/null +++ b/docs/validation/P2.md @@ -0,0 +1,39 @@ +# Phase P2 validation report — canonical model and raw X/Y assembly + +Date: 2026-09-27. Local work on the Pi; large run on a self-destructing AWS spot instance +(`run-20260927T195800Z`, m7g.4xlarge, 77 min, ≈ $0.34). Summary: `P2_omni2141_summary.json`. + +## Delivered + +| Item | Module | +|---|---| +| Raw X/Y + control assembly (plan §7.1), per-chip address-map cache | `signal/raw_xy.py` | +| IDAT ↔ manifest address-coverage gate (IN-027) | `signal/raw_xy.py` | +| Canonical Arrow schemas, deterministic Parquet (plan §5) | `core/schemas.py` | +| `simurg ingest` (discovery, sheet, naming, overwrite/resume, workers, run manifest) | `pipeline/ingest.py`, `cli/main.py` | + +## Acceptance (plan §18 P2: raw X/Y exactly equal to vendor GTC TOC 1000/1001, 100 % loci × ≥ 50 samples) + +| Check | Data | Result | +|---|---|---| +| Raw X/Y and controls vs oracle GTCs | 18 GSA v4 samples (gtc2vcf GTCs) | **EXACT, 18/18** | +| Raw X/Y vs **vendor** GTCs (AutoCall 1.6.2.2) | 2141 Omni2.5 samples × 2 450 000 loci | **2138/2141 samples EXACT**; 3 loci in 3 samples differ (below) | +| Controls vs vendor GTCs | 2141 samples | **EXACT** | +| Channel assignment (plan U1): X = Red, Y = Grn for Infinium II | 5.2 × 10⁹ locus values | **confirmed** (0 swapped loci) | +| `ingest` outputs identical across thread counts and layouts | synthetic + GSA | **byte-identical** | + +The gate "100 % loci × ≥ 50 samples" is met (2138 samples fully exact). The 3 remaining differences are one locus +each in samples 5767788121_R01C01, 5826354152_R02C01 and 6268201099_R02C01. In each, we mark the locus missing +(0, 0) but the vendor kept **one** channel (e.g. `rs178512`: vendor `[7086, 0]`). Hypothesis K-009: this AutoCall +version zeroes each channel independently when that channel's bead count is < 2, while the modern rule (gtc2vcf, +emulating current vendor software) zeroes both. Not yet confirmed (needs the bead counts of those three IDATs). + +## `simurg ingest` at scale (2141 Omni2.5 samples, 16 workers) + +| Metric | Value | +|---|---| +| Wall time | 18 min 52 s (0.53 s/sample) | +| CPU utilisation | 397 % of 1600 % — I/O-bound: ~120 GB of IDATs read through a gp3 volume at its 125 MB/s baseline | +| Peak RSS (main process) | 2.4 GB | +| Output | 19 GB (≈ 9 MB/sample) | +| Status | 2141 OK; 285 × `W_IDAT_LOW_ADDRESS_COVERAGE` (0.03 % of addresses absent), 4 × `W_IDAT_NO_SCAN_TIME` | diff --git a/docs/validation/P2_omni2141_summary.json b/docs/validation/P2_omni2141_summary.json new file mode 100644 index 0000000..ebcb11d --- /dev/null +++ b/docs/validation/P2_omni2141_summary.json @@ -0,0 +1,136 @@ +{ + "root": "/work/omni/data", + "started": "2026-09-27T20:21:53Z", + "n_gtc": 2141, + "n_idat_pairs": 2141, + "unpaired_idats": 0, + "gtc_ok": 2141, + "gtc_roundtrip_identical": 2141, + "gtc_versions": [ + 3 + ], + "gtc_autocall_versions": [ + "1.6.2.2" + ], + "gtc_failures": [], + "gtc_roundtrip_failures": [], + "idat_ok": 2141, + "idat_failures": [], + "idat_midblock_itemsizes": [ + 4 + ], + "idat_chip_types": [ + "BeadChip 4x10" + ], + "bpm": { + "num_loci": 2450000, + "names": true, + "snps": true, + "chroms": true, + "map_infos": true, + "address_a": true, + "assay_types": true, + "normalization_lookups": true, + "n_bins": 51 + }, + "p2_matched_gtc_idat": 2141, + "p2_ok": 2138, + "p2_total_diff_loci": 3, + "p2_swapped_loci": 0, + "p2_ours_missing_vendor_nonzero": 3, + "p2_vendor_zero_ours_present": 802, + "p2_failures": [ + { + "gtc": "/work/omni/data/5767788121/5767788121_R01C01.gtc", + "ok": false, + "n_diff": 1, + "n_missing": 3507, + "vendor_zero": 3506, + "swapped": 0, + "ours_missing_vendor_nonzero": 1, + "vendor_zero_ours_present": 0, + "controls_equal": true, + "examples": [ + { + "locus": "rs178512", + "assay": 0, + "ours": [ + 0, + 0 + ], + "vendor": [ + 7086, + 0 + ] + } + ] + }, + { + "gtc": "/work/omni/data/5826354152/5826354152_R02C01.gtc", + "ok": false, + "n_diff": 1, + "n_missing": 2320, + "vendor_zero": 2319, + "swapped": 0, + "ours_missing_vendor_nonzero": 1, + "vendor_zero_ours_present": 0, + "controls_equal": true, + "examples": [ + { + "locus": "SNP12-29082770", + "assay": 0, + "ours": [ + 0, + 0 + ], + "vendor": [ + 0, + 10569 + ] + } + ] + }, + { + "gtc": "/work/omni/data/6268201099/6268201099_R02C01.gtc", + "ok": false, + "n_diff": 1, + "n_missing": 7043, + "vendor_zero": 7042, + "swapped": 0, + "ours_missing_vendor_nonzero": 1, + "vendor_zero_ours_present": 0, + "controls_equal": true, + "examples": [ + { + "locus": "rs17133206", + "assay": 0, + "ours": [ + 0, + 0 + ], + "vendor": [ + 0, + 1827 + ] + } + ] + } + ], + "egt": { + "num_records": 2450000, + "dbv": 7, + "extra_floats": 15, + "record_mismatches": 0 + }, + "bpm_csv_problems": [], + "gtc_vs_bpm_num_snps": [ + 2450000, + 2450000 + ], + "gtc_transforms_vs_bins": [ + 51, + 51 + ], + "seconds": 2074.2, + "passed": false +}