diff --git a/README.md b/README.md index 5d12e7e..08e680d 100644 --- a/README.md +++ b/README.md @@ -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 [--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 [--json]`, `simurg ingest`, `simurg normalize`, `simurg genotype`. ## Quick start diff --git a/docs/known_unknowns.md b/docs/known_unknowns.md index 2c38c4a..647549c 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-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 | diff --git a/docs/validation/P4.md b/docs/validation/P4.md new file mode 100644 index 0000000..b77026b --- /dev/null +++ b/docs/validation/P4.md @@ -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 | diff --git a/docs/validation/P4_omni2141_summary.json b/docs/validation/P4_omni2141_summary.json new file mode 100644 index 0000000..cbe598f --- /dev/null +++ b/docs/validation/P4_omni2141_summary.json @@ -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 +} diff --git a/scripts/validate_p4_large.py b/scripts/validate_p4_large.py index f2ccf63..321ddb1 100644 --- a/scripts/validate_p4_large.py +++ b/scripts/validate_p4_large.py @@ -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 diff --git a/src/simurg/cli/main.py b/src/simurg/cli/main.py index 694e47b..9c2ecb1 100644 --- a/src/simurg/cli/main.py +++ b/src/simurg/cli/main.py @@ -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, diff --git a/src/simurg/genotype/sample.py b/src/simurg/genotype/sample.py index 1fdbe19..b698f04 100644 --- a/src/simurg/genotype/sample.py +++ b/src/simurg/genotype/sample.py @@ -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, + } diff --git a/src/simurg/pipeline/genotype.py b/src/simurg/pipeline/genotype.py index 1becc69..d3acdef 100644 --- a/src/simurg/pipeline/genotype.py +++ b/src/simurg/pipeline/genotype.py @@ -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}") diff --git a/tests/unit/test_genotype.py b/tests/unit/test_genotype.py index 3e3ac16..8fc5244 100644 --- a/tests/unit/test_genotype.py +++ b/tests/unit/test_genotype.py @@ -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