P4 validated on 2141 vendor GTCs; legacy AutoCall 1.6.2 preset
- 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
This commit is contained in:
kemal committed
1 parent
7032074bfe
commit
5d3e334ab7
9 files changed
+792
-29
No files matched your search
@@ -9,10 +9,10 @@ Personal research toolkit for Illumina Infinium genotyping BeadChip data
|
||||
|
||||
## Status
|
||||
|
||||
Phases **P0 (setup)**, **P1 (file parsers)**, **P2 (canonical model, raw X/Y)** and **P3 (normalization)** are
|
||||
complete — see `docs/validation/`. Each is validated on 2141 vendor GTCs of the 1000 Genomes Omni2.5 set
|
||||
(normalization transforms bit-identical).
|
||||
Available commands: `simurg version`, `simurg inspect <files…> [--json]`, `simurg ingest`, `simurg normalize`.
|
||||
Phases **P0 (setup)**, **P1 (file parsers)**, **P2 (canonical model, raw X/Y)**, **P3 (normalization)** and **P4 (genotype
|
||||
calling, GTC output)** are complete — see `docs/validation/`. Each is validated on 2141 vendor GTCs of the 1000 Genomes Omni2.5 set
|
||||
(normalization bit-identical; genotypes/scores identical in 2138/2141 samples).
|
||||
Available commands: `simurg version`, `simurg inspect <files…> [--json]`, `simurg ingest`, `simurg normalize`, `simurg genotype`.
|
||||
|
||||
## Quick start
|
||||
|
||||
|
||||
@@ -4,6 +4,9 @@ Plan §20 lists the open questions; this file records what the implementation ha
|
||||
|
||||
| ID | Date | Topic | Finding | Status |
|
||||
|---|---|---|---|---|
|
||||
| K-016 | 2026-09-28 | Sex rule of AutoCall 1.6.2.2 | Rule 1 (X heterozygosity) matches all 2141 vendor GTCs; rule 2 (X/Y intensity) differs in 12. | resolved |
|
||||
| K-015 | 2026-09-28 | Legacy GTC v3 statistics | AutoCall 1.6.2.2 writes GC10 = GC50 = 0 and #calls = #no-calls = 0 in GTC v3 (it does fill LogR dev and #intensity-only). gtc2vcf's AutoConvert 1.6.3.1 preset fills them, so this is specific to 1.6.2.x. Emulated by `--preset autocall-1.6.2`; the sample summary still reports the real statistics. | resolved; cloud re-run pending |
|
||||
| K-014 | 2026-09-28 | NaN sign | NaNs produced by libm (e.g. log2 of a negative extrapolated R) take the platform's default sign: + on aarch64, − on x86-64. SimurgArray always writes 0xFFC00000 (the reference's `-NAN`, x86 behaviour). The aarch64-built oracle differs in 1 locus of 9 GSA samples; oracle tests compare with NaNs canonicalised. | resolved |
|
||||
| K-013 | 2026-09-27 | Normalization version of public vendor GTCs | The 1000 Genomes Omni2.5 GTCs (AutoCall 1.6.2.2) are reproduced bit-for-bit by normalization 1.1.2 (109 191/109 191 bins, reserved fields included). | resolved |
|
||||
| K-012 | 2026-09-27 | FMA and bit-exactness | gcc on aarch64 contracts a*b+c into FMA by default; the gtc2vcf oracle then differs from the vendor by up to 4 ULP in some bins. Oracles are built with `-ffp-contract=off`. Any future C/Rust port (P9) must also forbid FMA contraction in exact mode. | resolved; constraint for P9 |
|
||||
| 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 |
|
||||
|
||||
@@ -0,0 +1,55 @@
|
||||
# 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 |
|
||||
@@ -0,0 +1,615 @@
|
||||
{
|
||||
"n_gtc": 2141,
|
||||
"n_compared": 2141,
|
||||
"errors": [],
|
||||
"byte_identical": 0,
|
||||
"samples_with_no_field_diff": 0,
|
||||
"samples_with_field_diff": {
|
||||
"gc10": 2141,
|
||||
"gc50": 2141,
|
||||
"num_calls": 2141,
|
||||
"num_no_calls": 2141,
|
||||
"slide_identifier": 2141,
|
||||
"raw_x": 1,
|
||||
"genotypes": 3,
|
||||
"genotype_scores": 3,
|
||||
"base_calls": 3,
|
||||
"call_rate": 2,
|
||||
"logr_dev": 3,
|
||||
"raw_y": 2,
|
||||
"normalization_transforms": 2
|
||||
},
|
||||
"genotype_diff_loci_total": 5,
|
||||
"score_diff_loci_total": 18229,
|
||||
"max_score_abs_diff": 0.9426311254501343,
|
||||
"cluster_files": [
|
||||
"HumanOmni2.5-4v1-Multi_B.egt"
|
||||
],
|
||||
"sex_agreement_by_version": {
|
||||
"1": 2141,
|
||||
"2": 2129
|
||||
},
|
||||
"ref_sex_counts": {
|
||||
"M": 1054,
|
||||
"F": 1087,
|
||||
"U": 0
|
||||
},
|
||||
"examples": [
|
||||
{
|
||||
"gtc": "/work/omni/data/5426511011/5426511011_R01C01.gtc",
|
||||
"cluster_file": "HumanOmni2.5-4v1-Multi_B.egt",
|
||||
"field_diffs": {
|
||||
"gc10": [
|
||||
"0.6847535967826843",
|
||||
"0.0"
|
||||
],
|
||||
"gc50": [
|
||||
"0.8497611880302429",
|
||||
"0.0"
|
||||
],
|
||||
"num_calls": [
|
||||
"2433210",
|
||||
"0"
|
||||
],
|
||||
"num_no_calls": [
|
||||
"16790",
|
||||
"0"
|
||||
],
|
||||
"slide_identifier": [
|
||||
"'5426511011'",
|
||||
"''"
|
||||
]
|
||||
},
|
||||
"sex_by_version": {
|
||||
"1": "F",
|
||||
"2": "F"
|
||||
},
|
||||
"sex_ref": "F",
|
||||
"byte_identical": false,
|
||||
"ref_layout_canonical": true,
|
||||
"metadata": {
|
||||
"autocall_version": "1.6.2.2",
|
||||
"cluster_file": "HumanOmni2.5-4v1-Multi_B.egt",
|
||||
"imaging_date": "Wednesday, August 25, 2010 7:11:03 PM"
|
||||
}
|
||||
},
|
||||
{
|
||||
"gtc": "/work/omni/data/5426511011/5426511011_R02C01.gtc",
|
||||
"cluster_file": "HumanOmni2.5-4v1-Multi_B.egt",
|
||||
"field_diffs": {
|
||||
"gc10": [
|
||||
"0.6842713356018066",
|
||||
"0.0"
|
||||
],
|
||||
"gc50": [
|
||||
"0.8497611880302429",
|
||||
"0.0"
|
||||
],
|
||||
"num_calls": [
|
||||
"2436202",
|
||||
"0"
|
||||
],
|
||||
"num_no_calls": [
|
||||
"13798",
|
||||
"0"
|
||||
],
|
||||
"slide_identifier": [
|
||||
"'5426511011'",
|
||||
"''"
|
||||
]
|
||||
},
|
||||
"sex_by_version": {
|
||||
"1": "M",
|
||||
"2": "M"
|
||||
},
|
||||
"sex_ref": "M",
|
||||
"byte_identical": false,
|
||||
"ref_layout_canonical": true,
|
||||
"metadata": {
|
||||
"autocall_version": "1.6.2.2",
|
||||
"cluster_file": "HumanOmni2.5-4v1-Multi_B.egt",
|
||||
"imaging_date": "Wednesday, August 25, 2010 7:24:37 PM"
|
||||
}
|
||||
},
|
||||
{
|
||||
"gtc": "/work/omni/data/5426511011/5426511011_R03C01.gtc",
|
||||
"cluster_file": "HumanOmni2.5-4v1-Multi_B.egt",
|
||||
"field_diffs": {
|
||||
"gc10": [
|
||||
"0.6842615604400635",
|
||||
"0.0"
|
||||
],
|
||||
"gc50": [
|
||||
"0.8497611880302429",
|
||||
"0.0"
|
||||
],
|
||||
"num_calls": [
|
||||
"2435069",
|
||||
"0"
|
||||
],
|
||||
"num_no_calls": [
|
||||
"14931",
|
||||
"0"
|
||||
],
|
||||
"slide_identifier": [
|
||||
"'5426511011'",
|
||||
"''"
|
||||
]
|
||||
},
|
||||
"sex_by_version": {
|
||||
"1": "F",
|
||||
"2": "F"
|
||||
},
|
||||
"sex_ref": "F",
|
||||
"byte_identical": false,
|
||||
"ref_layout_canonical": true,
|
||||
"metadata": {
|
||||
"autocall_version": "1.6.2.2",
|
||||
"cluster_file": "HumanOmni2.5-4v1-Multi_B.egt",
|
||||
"imaging_date": "Wednesday, August 25, 2010 7:37:41 PM"
|
||||
}
|
||||
},
|
||||
{
|
||||
"gtc": "/work/omni/data/5426511011/5426511011_R04C01.gtc",
|
||||
"cluster_file": "HumanOmni2.5-4v1-Multi_B.egt",
|
||||
"field_diffs": {
|
||||
"gc10": [
|
||||
"0.6863988637924194",
|
||||
"0.0"
|
||||
],
|
||||
"gc50": [
|
||||
"0.8497668504714966",
|
||||
"0.0"
|
||||
],
|
||||
"num_calls": [
|
||||
"2435856",
|
||||
"0"
|
||||
],
|
||||
"num_no_calls": [
|
||||
"14144",
|
||||
"0"
|
||||
],
|
||||
"slide_identifier": [
|
||||
"'5426511011'",
|
||||
"''"
|
||||
]
|
||||
},
|
||||
"sex_by_version": {
|
||||
"1": "F",
|
||||
"2": "F"
|
||||
},
|
||||
"sex_ref": "F",
|
||||
"byte_identical": false,
|
||||
"ref_layout_canonical": true,
|
||||
"metadata": {
|
||||
"autocall_version": "1.6.2.2",
|
||||
"cluster_file": "HumanOmni2.5-4v1-Multi_B.egt",
|
||||
"imaging_date": "Wednesday, August 25, 2010 7:45:50 PM"
|
||||
}
|
||||
},
|
||||
{
|
||||
"gtc": "/work/omni/data/5426511013/5426511013_R01C01.gtc",
|
||||
"cluster_file": "HumanOmni2.5-4v1-Multi_B.egt",
|
||||
"field_diffs": {
|
||||
"gc10": [
|
||||
"0.6843520402908325",
|
||||
"0.0"
|
||||
],
|
||||
"gc50": [
|
||||
"0.8497611880302429",
|
||||
"0.0"
|
||||
],
|
||||
"num_calls": [
|
||||
"2434354",
|
||||
"0"
|
||||
],
|
||||
"num_no_calls": [
|
||||
"15646",
|
||||
"0"
|
||||
],
|
||||
"slide_identifier": [
|
||||
"'5426511013'",
|
||||
"''"
|
||||
]
|
||||
},
|
||||
"sex_by_version": {
|
||||
"1": "M",
|
||||
"2": "M"
|
||||
},
|
||||
"sex_ref": "M",
|
||||
"byte_identical": false,
|
||||
"ref_layout_canonical": true,
|
||||
"metadata": {
|
||||
"autocall_version": "1.6.2.2",
|
||||
"cluster_file": "HumanOmni2.5-4v1-Multi_B.egt",
|
||||
"imaging_date": "Thursday, August 26, 2010 12:08:26 AM"
|
||||
}
|
||||
},
|
||||
{
|
||||
"gtc": "/work/omni/data/5426511013/5426511013_R02C01.gtc",
|
||||
"cluster_file": "HumanOmni2.5-4v1-Multi_B.egt",
|
||||
"field_diffs": {
|
||||
"gc10": [
|
||||
"0.685818076133728",
|
||||
"0.0"
|
||||
],
|
||||
"gc50": [
|
||||
"0.8497965335845947",
|
||||
"0.0"
|
||||
],
|
||||
"num_calls": [
|
||||
"2439588",
|
||||
"0"
|
||||
],
|
||||
"num_no_calls": [
|
||||
"10412",
|
||||
"0"
|
||||
],
|
||||
"slide_identifier": [
|
||||
"'5426511013'",
|
||||
"''"
|
||||
]
|
||||
},
|
||||
"sex_by_version": {
|
||||
"1": "F",
|
||||
"2": "F"
|
||||
},
|
||||
"sex_ref": "F",
|
||||
"byte_identical": false,
|
||||
"ref_layout_canonical": true,
|
||||
"metadata": {
|
||||
"autocall_version": "1.6.2.2",
|
||||
"cluster_file": "HumanOmni2.5-4v1-Multi_B.egt",
|
||||
"imaging_date": "Thursday, August 26, 2010 12:22:41 AM"
|
||||
}
|
||||
},
|
||||
{
|
||||
"gtc": "/work/omni/data/5426511013/5426511013_R03C01.gtc",
|
||||
"cluster_file": "HumanOmni2.5-4v1-Multi_B.egt",
|
||||
"field_diffs": {
|
||||
"gc10": [
|
||||
"0.6856454014778137",
|
||||
"0.0"
|
||||
],
|
||||
"gc50": [
|
||||
"0.8497673273086548",
|
||||
"0.0"
|
||||
],
|
||||
"num_calls": [
|
||||
"2442131",
|
||||
"0"
|
||||
],
|
||||
"num_no_calls": [
|
||||
"7869",
|
||||
"0"
|
||||
],
|
||||
"slide_identifier": [
|
||||
"'5426511013'",
|
||||
"''"
|
||||
]
|
||||
},
|
||||
"sex_by_version": {
|
||||
"1": "M",
|
||||
"2": "M"
|
||||
},
|
||||
"sex_ref": "M",
|
||||
"byte_identical": false,
|
||||
"ref_layout_canonical": true,
|
||||
"metadata": {
|
||||
"autocall_version": "1.6.2.2",
|
||||
"cluster_file": "HumanOmni2.5-4v1-Multi_B.egt",
|
||||
"imaging_date": "Thursday, August 26, 2010 12:37:24 AM"
|
||||
}
|
||||
},
|
||||
{
|
||||
"gtc": "/work/omni/data/5426511013/5426511013_R04C01.gtc",
|
||||
"cluster_file": "HumanOmni2.5-4v1-Multi_B.egt",
|
||||
"field_diffs": {
|
||||
"gc10": [
|
||||
"0.6856228113174438",
|
||||
"0.0"
|
||||
],
|
||||
"gc50": [
|
||||
"0.8497611880302429",
|
||||
"0.0"
|
||||
],
|
||||
"num_calls": [
|
||||
"2440841",
|
||||
"0"
|
||||
],
|
||||
"num_no_calls": [
|
||||
"9159",
|
||||
"0"
|
||||
],
|
||||
"slide_identifier": [
|
||||
"'5426511013'",
|
||||
"''"
|
||||
]
|
||||
},
|
||||
"sex_by_version": {
|
||||
"1": "F",
|
||||
"2": "F"
|
||||
},
|
||||
"sex_ref": "F",
|
||||
"byte_identical": false,
|
||||
"ref_layout_canonical": true,
|
||||
"metadata": {
|
||||
"autocall_version": "1.6.2.2",
|
||||
"cluster_file": "HumanOmni2.5-4v1-Multi_B.egt",
|
||||
"imaging_date": "Thursday, August 26, 2010 12:46:33 AM"
|
||||
}
|
||||
},
|
||||
{
|
||||
"gtc": "/work/omni/data/5426511028/5426511028_R01C01.gtc",
|
||||
"cluster_file": "HumanOmni2.5-4v1-Multi_B.egt",
|
||||
"field_diffs": {
|
||||
"gc10": [
|
||||
"0.6833689212799072",
|
||||
"0.0"
|
||||
],
|
||||
"gc50": [
|
||||
"0.8496665954589844",
|
||||
"0.0"
|
||||
],
|
||||
"num_calls": [
|
||||
"2435388",
|
||||
"0"
|
||||
],
|
||||
"num_no_calls": [
|
||||
"14612",
|
||||
"0"
|
||||
],
|
||||
"slide_identifier": [
|
||||
"'5426511028'",
|
||||
"''"
|
||||
]
|
||||
},
|
||||
"sex_by_version": {
|
||||
"1": "M",
|
||||
"2": "M"
|
||||
},
|
||||
"sex_ref": "M",
|
||||
"byte_identical": false,
|
||||
"ref_layout_canonical": true,
|
||||
"metadata": {
|
||||
"autocall_version": "1.6.2.2",
|
||||
"cluster_file": "HumanOmni2.5-4v1-Multi_B.egt",
|
||||
"imaging_date": "Thursday, August 26, 2010 8:13:23 AM"
|
||||
}
|
||||
},
|
||||
{
|
||||
"gtc": "/work/omni/data/5426511028/5426511028_R02C01.gtc",
|
||||
"cluster_file": "HumanOmni2.5-4v1-Multi_B.egt",
|
||||
"field_diffs": {
|
||||
"gc10": [
|
||||
"0.6857903599739075",
|
||||
"0.0"
|
||||
],
|
||||
"gc50": [
|
||||
"0.8497611880302429",
|
||||
"0.0"
|
||||
],
|
||||
"num_calls": [
|
||||
"2435889",
|
||||
"0"
|
||||
],
|
||||
"num_no_calls": [
|
||||
"14111",
|
||||
"0"
|
||||
],
|
||||
"slide_identifier": [
|
||||
"'5426511028'",
|
||||
"''"
|
||||
]
|
||||
},
|
||||
"sex_by_version": {
|
||||
"1": "F",
|
||||
"2": "F"
|
||||
},
|
||||
"sex_ref": "F",
|
||||
"byte_identical": false,
|
||||
"ref_layout_canonical": true,
|
||||
"metadata": {
|
||||
"autocall_version": "1.6.2.2",
|
||||
"cluster_file": "HumanOmni2.5-4v1-Multi_B.egt",
|
||||
"imaging_date": "Thursday, August 26, 2010 8:28:43 AM"
|
||||
}
|
||||
},
|
||||
{
|
||||
"gtc": "/work/omni/data/5426511028/5426511028_R03C01.gtc",
|
||||
"cluster_file": "HumanOmni2.5-4v1-Multi_B.egt",
|
||||
"field_diffs": {
|
||||
"gc10": [
|
||||
"0.6847374439239502",
|
||||
"0.0"
|
||||
],
|
||||
"gc50": [
|
||||
"0.8497611880302429",
|
||||
"0.0"
|
||||
],
|
||||
"num_calls": [
|
||||
"2435026",
|
||||
"0"
|
||||
],
|
||||
"num_no_calls": [
|
||||
"14974",
|
||||
"0"
|
||||
],
|
||||
"slide_identifier": [
|
||||
"'5426511028'",
|
||||
"''"
|
||||
]
|
||||
},
|
||||
"sex_by_version": {
|
||||
"1": "F",
|
||||
"2": "F"
|
||||
},
|
||||
"sex_ref": "F",
|
||||
"byte_identical": false,
|
||||
"ref_layout_canonical": true,
|
||||
"metadata": {
|
||||
"autocall_version": "1.6.2.2",
|
||||
"cluster_file": "HumanOmni2.5-4v1-Multi_B.egt",
|
||||
"imaging_date": "Thursday, August 26, 2010 8:40:16 AM"
|
||||
}
|
||||
},
|
||||
{
|
||||
"gtc": "/work/omni/data/5426511028/5426511028_R04C01.gtc",
|
||||
"cluster_file": "HumanOmni2.5-4v1-Multi_B.egt",
|
||||
"field_diffs": {
|
||||
"gc10": [
|
||||
"0.6849170923233032",
|
||||
"0.0"
|
||||
],
|
||||
"gc50": [
|
||||
"0.8497409820556641",
|
||||
"0.0"
|
||||
],
|
||||
"num_calls": [
|
||||
"2433211",
|
||||
"0"
|
||||
],
|
||||
"num_no_calls": [
|
||||
"16789",
|
||||
"0"
|
||||
],
|
||||
"slide_identifier": [
|
||||
"'5426511028'",
|
||||
"''"
|
||||
]
|
||||
},
|
||||
"sex_by_version": {
|
||||
"1": "F",
|
||||
"2": "F"
|
||||
},
|
||||
"sex_ref": "F",
|
||||
"byte_identical": false,
|
||||
"ref_layout_canonical": true,
|
||||
"metadata": {
|
||||
"autocall_version": "1.6.2.2",
|
||||
"cluster_file": "HumanOmni2.5-4v1-Multi_B.egt",
|
||||
"imaging_date": "Thursday, August 26, 2010 8:50:08 AM"
|
||||
}
|
||||
},
|
||||
{
|
||||
"gtc": "/work/omni/data/5426511029/5426511029_R01C01.gtc",
|
||||
"cluster_file": "HumanOmni2.5-4v1-Multi_B.egt",
|
||||
"field_diffs": {
|
||||
"gc10": [
|
||||
"0.683781087398529",
|
||||
"0.0"
|
||||
],
|
||||
"gc50": [
|
||||
"0.8497233390808105",
|
||||
"0.0"
|
||||
],
|
||||
"num_calls": [
|
||||
"2439124",
|
||||
"0"
|
||||
],
|
||||
"num_no_calls": [
|
||||
"10876",
|
||||
"0"
|
||||
],
|
||||
"slide_identifier": [
|
||||
"'5426511029'",
|
||||
"''"
|
||||
]
|
||||
},
|
||||
"sex_by_version": {
|
||||
"1": "M",
|
||||
"2": "M"
|
||||
},
|
||||
"sex_ref": "M",
|
||||
"byte_identical": false,
|
||||
"ref_layout_canonical": true,
|
||||
"metadata": {
|
||||
"autocall_version": "1.6.2.2",
|
||||
"cluster_file": "HumanOmni2.5-4v1-Multi_B.egt",
|
||||
"imaging_date": "Thursday, August 26, 2010 9:12:02 AM"
|
||||
}
|
||||
},
|
||||
{
|
||||
"gtc": "/work/omni/data/5426511029/5426511029_R02C01.gtc",
|
||||
"cluster_file": "HumanOmni2.5-4v1-Multi_B.egt",
|
||||
"field_diffs": {
|
||||
"gc10": [
|
||||
"0.6832002401351929",
|
||||
"0.0"
|
||||
],
|
||||
"gc50": [
|
||||
"0.8497028946876526",
|
||||
"0.0"
|
||||
],
|
||||
"num_calls": [
|
||||
"2437630",
|
||||
"0"
|
||||
],
|
||||
"num_no_calls": [
|
||||
"12370",
|
||||
"0"
|
||||
],
|
||||
"slide_identifier": [
|
||||
"'5426511029'",
|
||||
"''"
|
||||
]
|
||||
},
|
||||
"sex_by_version": {
|
||||
"1": "M",
|
||||
"2": "M"
|
||||
},
|
||||
"sex_ref": "M",
|
||||
"byte_identical": false,
|
||||
"ref_layout_canonical": true,
|
||||
"metadata": {
|
||||
"autocall_version": "1.6.2.2",
|
||||
"cluster_file": "HumanOmni2.5-4v1-Multi_B.egt",
|
||||
"imaging_date": "Thursday, August 26, 2010 9:26:14 AM"
|
||||
}
|
||||
},
|
||||
{
|
||||
"gtc": "/work/omni/data/5426511029/5426511029_R03C01.gtc",
|
||||
"cluster_file": "HumanOmni2.5-4v1-Multi_B.egt",
|
||||
"field_diffs": {
|
||||
"gc10": [
|
||||
"0.6829559803009033",
|
||||
"0.0"
|
||||
],
|
||||
"gc50": [
|
||||
"0.8496484756469727",
|
||||
"0.0"
|
||||
],
|
||||
"num_calls": [
|
||||
"2435690",
|
||||
"0"
|
||||
],
|
||||
"num_no_calls": [
|
||||
"14310",
|
||||
"0"
|
||||
],
|
||||
"slide_identifier": [
|
||||
"'5426511029'",
|
||||
"''"
|
||||
]
|
||||
},
|
||||
"sex_by_version": {
|
||||
"1": "F",
|
||||
"2": "F"
|
||||
},
|
||||
"sex_ref": "F",
|
||||
"byte_identical": false,
|
||||
"ref_layout_canonical": true,
|
||||
"metadata": {
|
||||
"autocall_version": "1.6.2.2",
|
||||
"cluster_file": "HumanOmni2.5-4v1-Multi_B.egt",
|
||||
"imaging_date": "Thursday, August 26, 2010 9:40:16 AM"
|
||||
}
|
||||
}
|
||||
],
|
||||
"params": {
|
||||
"norm_version": "1.1.2",
|
||||
"sex_version": 1
|
||||
},
|
||||
"seconds": 1174.6,
|
||||
"passed": false
|
||||
}
|
||||
@@ -37,8 +37,8 @@ ALIGN = (
|
||||
_W: dict[str, Any] = {}
|
||||
|
||||
|
||||
def _init(bpm: str, egt_dir: str, norm: str, sex_version: int) -> None:
|
||||
_W.update(bpm=bpm, egt_dir=egt_dir, norm=norm, sex_version=sex_version, ctx={})
|
||||
def _init(bpm: str, egt_dir: str, norm: str, sex_version: int, preset: str) -> None:
|
||||
_W.update(bpm=bpm, egt_dir=egt_dir, norm=norm, sex_version=sex_version, preset=preset, ctx={})
|
||||
|
||||
|
||||
def _ctx(cluster_field: str) -> Any:
|
||||
@@ -69,7 +69,7 @@ def _canon(a: np.ndarray) -> np.ndarray: # type: ignore[type-arg]
|
||||
def check(args: tuple[str, str, str]) -> dict[str, Any]:
|
||||
from simurg.formats.gtc import encode_gtc, read_gtc
|
||||
from simurg.formats.idat import read_idat
|
||||
from simurg.genotype.sample import GenotypeParams, genotype_sample
|
||||
from simurg.genotype.sample import genotype_sample, params_for_preset
|
||||
from simurg.genotype.sex import SexParams
|
||||
|
||||
gtc, grn, red = args
|
||||
@@ -77,8 +77,12 @@ def check(args: tuple[str, str, str]) -> dict[str, Any]:
|
||||
try:
|
||||
ref = read_gtc(gtc, check_write_complete=False)
|
||||
ctx = _ctx(ref.cluster_file)
|
||||
params = GenotypeParams(
|
||||
norm_version=_W["norm"], sex=SexParams(version=_W["sex_version"]), gtc_version=ref.version
|
||||
params = (
|
||||
params_for_preset(_W["preset"], gtc_version=ref.version)
|
||||
if _W["preset"]
|
||||
else params_for_preset(
|
||||
"current", norm_version=_W["norm"], sex=SexParams(version=_W["sex_version"]), gtc_version=ref.version
|
||||
)
|
||||
)
|
||||
g = genotype_sample(read_idat(grn), read_idat(red), ctx, params)
|
||||
res["cluster_file"] = ctx.cluster_file_name
|
||||
@@ -110,12 +114,13 @@ def check(args: tuple[str, str, str]) -> dict[str, Any]:
|
||||
"ploidy",
|
||||
"ploidy_type",
|
||||
"dx",
|
||||
"slide_identifier",
|
||||
):
|
||||
a, b = getattr(g, k), getattr(ref, k)
|
||||
same = a == b or (isinstance(a, float) and np.isnan(a) and np.isnan(b))
|
||||
if not same:
|
||||
diffs[k] = [repr(a), repr(b)]
|
||||
if ref.version >= 5 and g.slide_identifier != ref.slide_identifier: # field 1016 exists only in v5
|
||||
diffs["slide_identifier"] = [g.slide_identifier, ref.slide_identifier]
|
||||
if "genotype_scores" in diffs:
|
||||
d = np.abs(g.genotype_scores.astype(np.float64) - ref.genotype_scores.astype(np.float64))
|
||||
res["max_score_abs_diff"] = float(np.nanmax(d))
|
||||
@@ -155,6 +160,7 @@ def main() -> int:
|
||||
ap.add_argument("--out", required=True, type=Path)
|
||||
ap.add_argument("--norm-version", default="1.1.2")
|
||||
ap.add_argument("--sex-version", type=int, default=1)
|
||||
ap.add_argument("--preset", default="autocall-1.6.2", help="genotype preset; '' = use --norm/--sex-version")
|
||||
ap.add_argument("--workers", type=int, default=os.cpu_count() or 1)
|
||||
ap.add_argument("--limit", type=int, default=0)
|
||||
a = ap.parse_args()
|
||||
@@ -184,7 +190,7 @@ def main() -> int:
|
||||
a.workers,
|
||||
mp_context=mp.get_context("spawn"),
|
||||
initializer=_init,
|
||||
initargs=(str(a.bpm), str(a.egt_dir), a.norm_version, a.sex_version),
|
||||
initargs=(str(a.bpm), str(a.egt_dir), a.norm_version, a.sex_version, a.preset),
|
||||
) as ex:
|
||||
results = list(ex.map(check, triples, chunksize=2))
|
||||
|
||||
@@ -218,7 +224,7 @@ def main() -> int:
|
||||
},
|
||||
"ref_sex_counts": {s: sum(r.get("sex_ref") == s for r in results) for s in ("M", "F", "U")},
|
||||
"examples": [r for r in results if r.get("field_diffs")][:15],
|
||||
"params": {"norm_version": a.norm_version, "sex_version": a.sex_version},
|
||||
"params": {"preset": a.preset, "norm_version": a.norm_version, "sex_version": a.sex_version},
|
||||
"seconds": round(time.time() - t0, 1),
|
||||
}
|
||||
summary["passed"] = summary["byte_identical"] == len(results) > 0
|
||||
|
||||
+19
-5
@@ -284,9 +284,18 @@ def genotype_cmd(
|
||||
sheet: Annotated[Path | None, typer.Option("--sheet", help="Sample sheet (CSV).")] = None,
|
||||
csv: Annotated[Path | None, typer.Option("--csv", help="CSV manifest, checked against the BPM.")] = None,
|
||||
cutoff: Annotated[float, typer.Option("--cutoff", help="GenCall score below which a call is a no-call.")] = 0.15,
|
||||
norm_version: Annotated[str, typer.Option("--norm-version", help="1.2.0 (current) or 1.1.2")] = "1.2.0",
|
||||
sex_version: Annotated[int, typer.Option("--sex-version", min=1, max=2, help="Sex estimate rule (1 or 2).")] = 2,
|
||||
gtc_version: Annotated[int, typer.Option("--gtc-version", help="GTC format version (3 or 5).")] = 5,
|
||||
preset: Annotated[
|
||||
str, typer.Option("--preset", help="current (default) or autocall-1.6.2 (legacy vendor behaviour)")
|
||||
] = "current",
|
||||
norm_version: Annotated[
|
||||
str | None, typer.Option("--norm-version", help="1.2.0 or 1.1.2 (default: from --preset)")
|
||||
] = None,
|
||||
sex_version: Annotated[
|
||||
int | None, typer.Option("--sex-version", min=1, max=2, help="Sex estimate rule 1 or 2 (default: preset)")
|
||||
] = None,
|
||||
gtc_version: Annotated[
|
||||
int | None, typer.Option("--gtc-version", help="GTC format 3 or 5 (default: preset)")
|
||||
] = None,
|
||||
autocall_date: Annotated[
|
||||
str | None, typer.Option("--autocall-date", help="Fixed value for GTC field 201 (reproducible output).")
|
||||
] = None,
|
||||
@@ -303,9 +312,13 @@ def genotype_cmd(
|
||||
|
||||
if layout not in ("nested", "flat"):
|
||||
raise typer.BadParameter("--layout must be nested or flat")
|
||||
if gtc_version not in (3, 5):
|
||||
from simurg.genotype.sample import PRESETS
|
||||
|
||||
if preset not in PRESETS:
|
||||
raise typer.BadParameter(f"--preset must be one of {sorted(PRESETS)}")
|
||||
if gtc_version not in (None, 3, 5):
|
||||
raise typer.BadParameter("--gtc-version must be 3 or 5")
|
||||
if norm_version not in ("1.2.0", "1.1.2"):
|
||||
if norm_version not in (None, "1.2.0", "1.1.2"):
|
||||
raise typer.BadParameter("--norm-version must be 1.2.0 or 1.1.2")
|
||||
rep = Reporter(
|
||||
command="genotype",
|
||||
@@ -327,6 +340,7 @@ def genotype_cmd(
|
||||
layout=layout,
|
||||
on_duplicate_idat=on_duplicate_idat,
|
||||
cutoff=cutoff,
|
||||
preset=preset,
|
||||
norm_version=norm_version,
|
||||
sex_version=sex_version,
|
||||
gtc_version=gtc_version,
|
||||
|
||||
@@ -39,6 +39,27 @@ class GenotypeParams:
|
||||
gtc_version: int = 5
|
||||
autocall_version: str = f"{PRODUCT_NAME} {__version__}"
|
||||
autocall_date: str | None = None # None: run start time (invariant format, OUT-001)
|
||||
legacy_sample_stats: bool = False # AutoCall 1.6.2.x left GC10, GC50, #calls and #no-calls at 0 (K-015)
|
||||
|
||||
|
||||
PRESETS: dict[str, dict[str, object]] = {
|
||||
# default: current vendor software (normalization 1.2.0 / GenTrain 3, sex rule 2, GTC v5)
|
||||
"current": {},
|
||||
# AutoCall 1.6.2.x as used for the public 1000 Genomes Omni2.5 GTCs (validated on 2141 samples, P4)
|
||||
"autocall-1.6.2": {
|
||||
"norm_version": "1.1.2",
|
||||
"sex": SexParams(version=1),
|
||||
"gtc_version": 3,
|
||||
"legacy_sample_stats": True,
|
||||
},
|
||||
}
|
||||
|
||||
|
||||
def params_for_preset(name: str, **overrides: object) -> GenotypeParams:
|
||||
if name not in PRESETS:
|
||||
raise ValueError(f"unknown preset {name!r}; choose from {sorted(PRESETS)}")
|
||||
values = {**PRESETS[name], **{k: v for k, v in overrides.items() if v is not None}}
|
||||
return GenotypeParams(**values) # type: ignore[arg-type]
|
||||
|
||||
|
||||
@dataclass
|
||||
@@ -193,4 +214,23 @@ def genotype_sample(
|
||||
g.sample_data["num_calls"] = calls.num_calls
|
||||
g.sample_data["num_no_calls"] = m.num_loci - calls.num_intensity_only - calls.num_calls
|
||||
g.sample_data["num_intensity_only"] = calls.num_intensity_only
|
||||
if params.legacy_sample_stats:
|
||||
g.gc10 = 0.0
|
||||
g.sample_data["p50gc"] = 0.0
|
||||
g.sample_data["num_calls"] = 0
|
||||
g.sample_data["num_no_calls"] = 0
|
||||
return g
|
||||
|
||||
|
||||
def summary_stats(g: Gtc, intensity_only: npt.NDArray[np.bool_], cutoff: float) -> dict[str, float | int]:
|
||||
"""Sample statistics recomputed from calls and scores (independent of legacy GTC field conventions)."""
|
||||
gs = np.sort(g.genotype_scores[(g.genotype_scores > F32(cutoff)) & ~intensity_only])
|
||||
n_calls = int((g.genotypes > 0).sum())
|
||||
n_io = int(intensity_only.sum())
|
||||
return {
|
||||
"gc10": float(percentile(gs, 10)),
|
||||
"gc50": float(percentile(gs, 50)),
|
||||
"n_calls": n_calls,
|
||||
"n_nocalls": g.num_snps - n_io - n_calls,
|
||||
"n_intensity_only": n_io,
|
||||
}
|
||||
@@ -27,7 +27,14 @@ from simurg.formats.gtc import write_gtc
|
||||
from simurg.formats.idat import check_pair, read_idat
|
||||
from simurg.genotype.baf_lrr import chrom_core
|
||||
from simurg.genotype.clusters import align_clusters
|
||||
from simurg.genotype.sample import GenotypeParams, ManifestContext, genotype_sample, manifest_context
|
||||
from simurg.genotype.sample import (
|
||||
GenotypeParams,
|
||||
ManifestContext,
|
||||
genotype_sample,
|
||||
manifest_context,
|
||||
params_for_preset,
|
||||
summary_stats,
|
||||
)
|
||||
from simurg.genotype.sex import SexParams
|
||||
from simurg.pipeline.inputs import check_paths, resolve_sample_set, sample_dir, sheet_values
|
||||
from simurg.signal.normalization.apply import normalize, polar
|
||||
@@ -77,9 +84,10 @@ class GenotypeOptions:
|
||||
layout: str = "nested"
|
||||
on_duplicate_idat: str = "error"
|
||||
cutoff: float = 0.15
|
||||
norm_version: str = "1.2.0"
|
||||
sex_version: int = 2
|
||||
gtc_version: int = 5
|
||||
preset: str = "current"
|
||||
norm_version: str | None = None # None: from the preset
|
||||
sex_version: int | None = None
|
||||
gtc_version: int | None = None
|
||||
autocall_date: str | None = None
|
||||
assay_control_probes: tuple[str, ...] = ("TGA_Ctrl_5716",)
|
||||
argv: list[str] = field(default_factory=list)
|
||||
@@ -133,11 +141,7 @@ def _process(job: Job) -> dict[str, Any]:
|
||||
call_rate=float(g.call_rate),
|
||||
logr_dev=float(g.logr_dev),
|
||||
sex_est=g.gender_char,
|
||||
gc10=float(g.gc10),
|
||||
gc50=float(g.gc50),
|
||||
n_calls=g.num_calls,
|
||||
n_nocalls=g.num_no_calls,
|
||||
n_intensity_only=g.num_intensity_only,
|
||||
**summary_stats(g, ctx.manifest.intensity_only, params.cutoff),
|
||||
p05_x=int(g.percentiles_x[0]),
|
||||
p50_x=int(g.percentiles_x[1]),
|
||||
p95_x=int(g.percentiles_x[2]),
|
||||
@@ -261,13 +265,17 @@ def run_genotype(opts: GenotypeOptions, rep: Reporter) -> int:
|
||||
auto_1_22 = np.array([chrom_core(c) in {str(i) for i in range(1, 23)} for c in manifest.chroms])
|
||||
name_idx = manifest.name_index()
|
||||
probes = {p: name_idx[p] for p in opts.assay_control_probes if p in name_idx}
|
||||
params = GenotypeParams(
|
||||
params = params_for_preset(
|
||||
opts.preset,
|
||||
cutoff=opts.cutoff,
|
||||
norm_version=opts.norm_version,
|
||||
sex=SexParams(version=opts.sex_version),
|
||||
sex=SexParams(version=opts.sex_version) if opts.sex_version else None,
|
||||
gtc_version=opts.gtc_version,
|
||||
autocall_date=opts.autocall_date,
|
||||
)
|
||||
rm.set_option("effective_norm_version", params.norm_version, f"preset:{opts.preset}")
|
||||
rm.set_option("effective_sex_version", params.sex.version, f"preset:{opts.preset}")
|
||||
rm.set_option("effective_gtc_version", params.gtc_version, f"preset:{opts.preset}")
|
||||
started = _dt.datetime.now()
|
||||
|
||||
workers = max(1, min(opts.threads or (os.cpu_count() or 1), len(jobs) or 1))
|
||||
@@ -335,7 +343,7 @@ def run_genotype(opts: GenotypeOptions, rep: Reporter) -> int:
|
||||
"warnings": rep.warning_counts(),
|
||||
}
|
||||
rm.add_output(summary_tsv, SUMMARY_SCHEMA)
|
||||
rm.schema_versions["*.gtc"] = f"gtc/{opts.gtc_version}"
|
||||
rm.schema_versions["*.gtc"] = f"gtc/{params.gtc_version}"
|
||||
rm.status = "completed"
|
||||
rm.write(opts.out / "run_manifest.json")
|
||||
rep.info(f"done: {ok}/{len(out_rows)} samples genotyped; summary {summary_tsv}")
|
||||
|
||||
@@ -233,3 +233,25 @@ def test_cluster_alignment_errors(tmp_path: Path) -> None:
|
||||
assert e.value.code == "E_EGT_MANIFEST_MISMATCH"
|
||||
c = align_clusters(read_egt(tmp_path / "small.egt"), read_bpm(f["bpm"]), allow_missing=True)
|
||||
assert int(c.present.sum()) == 1
|
||||
|
||||
|
||||
def test_legacy_preset(tmp_path: Path) -> None:
|
||||
from simurg.cli.main import app
|
||||
|
||||
f = _synthetic_run(tmp_path)
|
||||
out = tmp_path / "legacy"
|
||||
r = CliRunner().invoke(app, ["genotype", "--quiet", "--bpm", str(f["bpm"]), "--egt", str(f["egt"]),
|
||||
"--idats", str(f["idats"]), "--out", str(out), "--threads", "1",
|
||||
"--preset", "autocall-1.6.2"])
|
||||
assert r.exit_code == 0, r.output
|
||||
gtc = next(out.rglob("*.gtc"))
|
||||
g = read_gtc(gtc)
|
||||
assert g.version == 3 and g.gc10 == 0 and g.num_calls == 0 and g.gc50 == 0 # vendor 1.6.2.x conventions
|
||||
summary = json.loads((out / "sample_summary.json").read_text())
|
||||
s0 = summary["samples"][0]
|
||||
assert s0["n_calls"] > 0 and s0["gc10"] > 0 # the summary still reports the real statistics
|
||||
run = json.loads((out / "run_manifest.json").read_text())
|
||||
assert run["effective_options"]["effective_norm_version"]["value"] == "1.1.2"
|
||||
bad = CliRunner().invoke(app, ["genotype", "--bpm", str(f["bpm"]), "--egt", str(f["egt"]), "--idats",
|
||||
str(f["idats"]), "--out", str(tmp_path / "x"), "--preset", "nope"])
|
||||
assert bad.exit_code == 2
|
||||
Reference in new issue
Block a user