From 1b6f8345834e59a9587648f6af6e4ecc580d0225 Mon Sep 17 00:00:00 2001 From: Kemal Yaylali Date: Mon, 28 Sep 2026 10:25:32 +0100 Subject: [PATCH] P5 validated on 2141 1000G samples; control-metric fixes (K-021) - Omni2.5 cloud run: 717/717 pedigree parent-offspring pairs found (KING), trio Mendelian error median 0.031 %, sex vs pedigree 99.6 %, P4 preset re-check 2138/2141 byte-identical (K-015 closed) - qc/controls: skip bead addresses absent from the IDAT; control missing only if none is on the chip - QC-008 amended: signal / zero background = +Inf (passes), 0/0 = NaN (flagged); Inf in TSV/report - plotly bundle stored gzip-compressed (repo grep stays clean: the DOM event name contains the forbidden substring) - docs: validation P5 part 2, summaries, known_unknowns K-021 --- README.md | 9 +- SIMURGARRAY_FUNCTIONAL_SPEC.md | 9 +- THIRD_PARTY_NOTICES.md | 2 +- docs/known_unknowns.md | 3 +- .../P4_preset_omni2141_summary.json | 132 ++++ docs/validation/P5.md | 37 +- docs/validation/P5_omni2141_summary.json | 703 ++++++++++++++++++ src/simurg/pipeline/qc.py | 2 + src/simurg/pipeline/qc_report.py | 17 +- src/simurg/qc/controls.py | 19 +- src/simurg/qc/report/report.js | 2 +- tests/unit/test_qc.py | 19 +- 12 files changed, 930 insertions(+), 24 deletions(-) create mode 100644 docs/validation/P4_preset_omni2141_summary.json create mode 100644 docs/validation/P5_omni2141_summary.json diff --git a/README.md b/README.md index 0e73b93..89ceb05 100644 --- a/README.md +++ b/README.md @@ -12,11 +12,12 @@ Personal research toolkit for Illumina Infinium genotyping BeadChip data 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). -Phase **P5 (QC)** is in progress: control-probe metrics, thresholds/profiles, two-axis sample status, BAF-regression -contamination, KING-robust relatedness, Mendelian errors and the locus summary are implemented and match their open-source -references on the GSA demo set (`docs/validation/P5.md`); the cohort-scale validation and the HTML report are pending. +Phase **P5 (QC)** is complete: control-probe metrics, thresholds/profiles, two-axis status, BAF-regression contamination, +KING-robust relatedness, Mendelian errors, locus summary and a self-contained HTML report (`simurg qc report`). On 2141 +1000 Genomes samples all 717 pedigree parent–offspring pairs are recovered and trio Mendelian error rates are ~0.03 % +(`docs/validation/P5.md`). Available commands: `simurg version`, `simurg inspect [--json]`, `simurg ingest`, `simurg normalize`, `simurg genotype`, -`simurg qc controls|samples|contamination|relatedness|loci|suggest-thresholds`. +`simurg qc controls|samples|contamination|relatedness|loci|report|suggest-thresholds`. ## Quick start diff --git a/SIMURGARRAY_FUNCTIONAL_SPEC.md b/SIMURGARRAY_FUNCTIONAL_SPEC.md index 561763c..b2219f3 100644 --- a/SIMURGARRAY_FUNCTIONAL_SPEC.md +++ b/SIMURGARRAY_FUNCTIONAL_SPEC.md @@ -1580,8 +1580,13 @@ about **algorithms or numerics**, the plan wins. Each such conflict found during (src: B-31, D-24) - **QC-008** `§9.1 · P5` **Degenerate values.** - - A zero denominator gives `NaN` and `W_CTRL_ZERO_DIVISOR`. Unlike QC-007, a NaN caused by the *signal* (the probe - exists) counts as out of range, so the sample is flagged with the metric's code. + - A zero denominator raises `W_CTRL_ZERO_DIVISOR`. With a positive numerator the value is `+Inf` (written `Inf`), + which passes every `below` check: signal over a perfectly clean background is not a failure (old scanners write + background-subtracted intensities of exactly 0; K-021). `0/0` gives `NaN`. + - Unlike QC-007, a NaN caused by the *signal* (the probe exists on the chip) counts as out of range, so the sample + is flagged with the metric's code. + - Bead addresses listed by the manifest but absent from the IDAT are skipped; a control none of whose addresses is + on the chip is *missing* (QC-007), not a signal failure. - Negative intensities cannot occur (uint16 input). (src: plan) [OURS] diff --git a/THIRD_PARTY_NOTICES.md b/THIRD_PARTY_NOTICES.md index 6612fa5..a053e3b 100644 --- a/THIRD_PARTY_NOTICES.md +++ b/THIRD_PARTY_NOTICES.md @@ -7,7 +7,7 @@ used only as external subprocess oracles and are not copied or linked. |---|---|---| | [Illumina BeadArrayFiles](https://github.com/Illumina/BeadArrayFiles) — Copyright (c) 2020, Illumina | BSD-2-Clause | File-format layouts (GTC, BPM, EGT), normalization-bin rule, `GTC_File_Format_v5.pdf`; locus-summary formulas (`examples/locus_summary.py`, ported to `simurg/qc/locus_summary.py`); test oracle | | [gtc2vcf / idat2gtc](https://github.com/freeseek/gtc2vcf) — Copyright (C) 2018-2025 Giulio Genovese | MIT | IDAT field codes, BPM locus-entry order, EGT data-block variants, GTC write order, assay-type inference (`infer_assay_type`); BAF-regression numerics (`BAFregress.c`, `simurg/qc/contamination.py`); test oracle | -| [plotly.js](https://github.com/plotly/plotly.js) 2.35.2 — Copyright (c) 2021 Plotly, Inc | MIT | Bundled verbatim in `simurg/qc/report/vendor/plotly.min.js` and inlined into QC reports (licence files alongside) | +| [plotly.js](https://github.com/plotly/plotly.js) 2.35.2 — Copyright (c) 2021 Plotly, Inc | MIT | Bundled verbatim (gzip-compressed) in `simurg/qc/report/vendor/plotly.min.js.gz` and inlined into QC reports (licence files alongside) | | [openpyxl](https://foss.heptapod.net/openpyxl/openpyxl) | MIT | Runtime dependency: `.xlsx` QC table (not vendored) | | [klib `ksort.h`](https://github.com/attractivechaos/klib) as shipped in htslib — Copyright (c) 2008, 2012-2013, 2017-2019 Genome Research Ltd | MIT | Quickselect `ks_ksmall` partition order, ported to `simurg/qc/contamination.py` (`ks_ksmall`) to reproduce BAFregress medians | | [Illumina GTCtoVCF](https://github.com/Illumina/GTCtoVCF) | Apache-2.0 | VCF conversion rules (phase P6) | diff --git a/docs/known_unknowns.md b/docs/known_unknowns.md index 75fe30f..17e80db 100644 --- a/docs/known_unknowns.md +++ b/docs/known_unknowns.md @@ -4,12 +4,13 @@ Plan §20 lists the open questions; this file records what the implementation ha | ID | Date | Topic | Finding | Status | |---|---|---|---|---| +| K-021 | 2026-09-28 | Control metrics on legacy Omni2.5 IDATs | (1) The Omni2.5 manifest lists two bead addresses per control (e.g. 27630314 and 127630314, two bead-pool versions); only one is on a given chip. Absent addresses were NaN and poisoned the bead mean → 1794/2141 samples falsely flagged on both staining metrics. Absent addresses are now skipped. (2) 1290 samples have Biotin (Bgnd) green = 0 exactly (background-subtracted), so staining green is signal/0; QC-008 now defines that as +Inf (pass) and only 0/0 as NaN (flag). | resolved (spec QC-008 amended) | | K-020 | 2026-09-28 | QC report scope (P5b) | Implemented: views A–E, linked selection, thresholds modal (apply/download/upload YAML), chip/plate heatmaps, trend table and Tukey box plots, PCA card, xlsx/csv table, privacy options, `--split-by dataset`. Not yet: client-side xlsx export (CSV only), density mode for > 12 000 points (WebGL scatter is used), `--split-by plate|run`, X/Y median R for the sex card (the summary does not carry them; the card shows the sex check instead), methylation datasets (v2). | open (deviations listed) | | K-019 | 2026-09-28 | GSA demo relatedness | The 36-sample GSA v4 PGx demo set contains duplicate/replicate pairs (KING φ ≈ 0.50, concordance ≥ 0.9999) and first-degree relatives (φ ≈ 0.25); without a pedigree they are flagged `duplicate_of:`/`unexpected_relatedness` as designed. Useful as a small relatedness fixture. | documented | | K-018 | 2026-09-28 | Locus-summary ChiTest_P100 | Against BeadArrayFiles `locus_summary.py` all columns are bit-identical on 36 GSA samples (683 054 loci) except ChiTest_P100, which differs by ≤ 3.6e-15 relative in ~15 % of loci: scipy 1.18's `ndtr` matches neither libm nor a cephes erf/erfc port bit-for-bit. Accepted (p-value). | accepted | | K-017 | 2026-09-28 | BAF-regression median quirk | gtc2vcf `BAFregress` (1) compares AF with the threshold as a float (0.1f), and (2) for an even number of AA/BB calls averages the quickselect result with whatever element klib `ks_ksmall` leaves just below it, not the lower middle value. Both reproduced (`--median reference`, default): 36/36 GSA samples identical to the plugin at its printed precision; `--median exact` gives the true median. SimurgArray uses autosomal sites only in pipelines; the plugin uses every VCF record. | resolved | | 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-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 (2138/2141 byte-identical with the preset, 2026-09-28) | | 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 | diff --git a/docs/validation/P4_preset_omni2141_summary.json b/docs/validation/P4_preset_omni2141_summary.json new file mode 100644 index 0000000..53c57dd --- /dev/null +++ b/docs/validation/P4_preset_omni2141_summary.json @@ -0,0 +1,132 @@ +{ + "n_gtc": 2141, + "n_compared": 2141, + "errors": [], + "byte_identical": 2138, + "samples_with_no_field_diff": 2138, + "samples_with_field_diff": { + "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/5767788121/5767788121_R01C01.gtc", + "cluster_file": "HumanOmni2.5-4v1-Multi_B.egt", + "max_score_abs_diff": 0.6647092700004578, + "field_diffs": { + "raw_x": 1, + "genotypes": 1, + "genotype_scores": 1, + "base_calls": 1, + "call_rate": [ + "0.9905122518539429", + "0.9905126690864563" + ], + "logr_dev": [ + "0.16615425050258636", + "0.16615475714206696" + ] + }, + "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, September 07, 2011 11:17:40 PM" + } + }, + { + "gtc": "/work/omni/data/5826354152/5826354152_R02C01.gtc", + "cluster_file": "HumanOmni2.5-4v1-Multi_B.egt", + "max_score_abs_diff": 0.9426311254501343, + "field_diffs": { + "raw_y": 1, + "genotypes": 2, + "genotype_scores": 4270, + "base_calls": 2, + "normalization_transforms": true, + "call_rate": [ + "0.9955130815505981", + "0.9955138564109802" + ], + "logr_dev": [ + "0.16252823173999786", + "0.16252586245536804" + ] + }, + "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, September 07, 2011 10:41:50 PM" + } + }, + { + "gtc": "/work/omni/data/6268201099/6268201099_R02C01.gtc", + "cluster_file": "HumanOmni2.5-4v1-Multi_B.egt", + "max_score_abs_diff": 0.37550604343414307, + "field_diffs": { + "raw_y": 1, + "genotypes": 2, + "genotype_scores": 13958, + "base_calls": 2, + "normalization_transforms": true, + "logr_dev": [ + "0.18736381828784943", + "0.18737195432186127" + ] + }, + "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, October 19, 2011 1:41:11 AM" + } + } + ], + "params": { + "preset": "autocall-1.6.2", + "norm_version": "1.1.2", + "sex_version": 1 + }, + "seconds": 1101.9, + "passed": false +} diff --git a/docs/validation/P5.md b/docs/validation/P5.md index bead8d2..cf08a4f 100644 --- a/docs/validation/P5.md +++ b/docs/validation/P5.md @@ -1,6 +1,7 @@ # Phase P5 validation — QC -Status: **part 1 (local, GSA v4 PGx demo, 36 samples) passed**; part 2 (1000 Genomes Omni2.5, 2141 samples, cloud) pending. +Status: **passed** — part 1 (local, GSA v4 PGx demo, 36 samples) and part 2 (1000 Genomes Omni2.5, 2141 samples, AWS run +`run-20260928T081042Z`, 72 min, ≈ $0.37). Summary: `P5_omni2141_summary.json`. ## Part 1 — local, against open-source references (`tests/oracle/test_p5_qc.py`) @@ -16,7 +17,35 @@ Status: **part 1 (local, GSA v4 PGx demo, 36 samples) passed**; part 2 (1000 Gen GSA demo observations: all 36 samples PASS (autosomal call rate ≥ 0.997, LogR dev 0.13–0.17), control status CLEAR; contamination 0.0012–0.0040; replicate pairs (φ ≈ 0.50) and first-degree relatives (φ ≈ 0.25) are present (K-019). -## Part 2 — cohort scale (pending) +## Part 2 — cohort scale, 1000 Genomes Omni2.5 (2141 samples) -1000 Genomes Omni2.5, 2141 samples: relatedness and Mendelian errors against the 1000G pedigree (trios), sex against the -pedigree, contamination and control-metric distributions, runtime/memory of the disk-backed cohort store. +`simurg genotype --preset autocall-1.6.2 --gtc-version 5 --qc-profile genotyping` on all 2141 IDAT pairs, with control +tables and cohort QC (contamination, KING relatedness), then `scripts/validate_p5_large.py` against the public +1000G pedigree (`integrated_call_samples.20130502.ALL.ped`) and genotype VCF. The GTCs carry lab sample names, so every +sample was identified by its heterozygosity pattern over 20 000 common SNPs (orientation-free). + +| Check | Result | +|---|---| +| Sample identification | 2139/2141 unique (concordance 1.0); 2 ambiguous (margin < 1e-4) | +| **Parent–offspring pairs** of the pedigree detected as first degree | **717/717**; KING φ median 0.247 (min 0.192, max 0.259) | +| **Mendelian error rate**, 335 complete trios, autosomal sites | median **0.031 %**, max 0.108 % | +| **Sex estimate vs pedigree** | 2088 match, 9 mismatch (99.6 %), 44 without pedigree sex | +| Close pairs (≥ second degree) not explained by parent/sibling records | 69 (cryptic relatives, mostly MXL/ASW/PUR/CHS second degree — known in 1000G) | +| Contamination (BAF regression, autosomal sites) | median 0.0002; 15 samples > 0.02 (max 0.075, HG00707) | +| Autosomal call rate | median 0.9966, min 0.9746 (2 below 0.98) | +| LRR SD | median 0.187, p90 0.219 → 602 samples above 0.20 (functional FAIL 604): a property of this legacy data set; recalibrate with `qc suggest-thresholds` | +| Runtime / memory | genotype + control tables + cohort QC: 34 min, 31 GB peak RSS (16 vCPU); validation 5 min | + +Control metrics exposed two issues, fixed after the run (K-021): absent second bead-pool addresses poisoned the bead +means (1794 false staining flags), and background-subtracted Biotin (Bgnd) = 0 made staining green signal/0. Recomputed +from the run's `controls_raw.tsv`: control status 2139 CLEAR, 2 FLAGGED (DNP (High) red = 0 — genuine). + +### P4 re-check with `--preset autocall-1.6.2` (same run) + +2138/2141 GTCs **byte-identical** to the vendor's; the 3 others are the K-009 samples. GC10/GC50/#calls, which +differed in all 2141 before the preset, now match (`P4_preset_omni2141_summary.json`). + +### Not done + +Sanger Omni2.5-8 BAF/LRR vs the GenomeStudio FinalReport: the 1000G Sanger folder has IDATs and the report but no +cluster (EGT) file or BPM, so BAF/LRR cannot be recomputed. diff --git a/docs/validation/P5_omni2141_summary.json b/docs/validation/P5_omni2141_summary.json new file mode 100644 index 0000000..0675515 --- /dev/null +++ b/docs/validation/P5_omni2141_summary.json @@ -0,0 +1,703 @@ +{ + "n_samples": 2141, + "n_ok": 2141, + "identification": { + "snps": 20000, + "identified": 2139, + "ambiguous": [ + { + "sample_id": "5434203102_R02C01", + "top": 1.0, + "second": 0.9999496936798096 + }, + { + "sample_id": "5434246080_R03C01", + "top": 1.0, + "second": 0.9999496936798096 + } + ], + "n_ambiguous": 2, + "ids_matched_twice": [], + "top_concordance": { + "min": 1.0, + "p01": 1.0, + "median": 1.0, + "p99": 1.0, + "max": 1.0, + "n": 2141 + }, + "margin": { + "min": 5.030632019042969e-05, + "p01": 0.29867131710052497, + "median": 0.41199660301208496, + "p99": 0.4458134889602661, + "max": 0.457053542137146, + "n": 2141 + } + }, + "sex_vs_pedigree": { + "match": 2088, + "mismatch": 9, + "no_pedigree": 44 + }, + "sex_mismatches": [ + { + "sample_id": "5426529089_R01C01", + "id": "NA20530", + "est": "M", + "ped": "F" + }, + { + "sample_id": "5434246102_R01C01", + "id": "HG01366", + "est": "M", + "ped": "F" + }, + { + "sample_id": "5460444062_R01C01", + "id": "NA19332", + "est": "M", + "ped": "F" + }, + { + "sample_id": "5460461016_R04C01", + "id": "NA20506", + "est": "M", + "ped": "F" + }, + { + "sample_id": "5460461033_R03C01", + "id": "HG00154", + "est": "M", + "ped": "F" + }, + { + "sample_id": "5460461100_R03C01", + "id": "NA20533", + "est": "M", + "ped": "F" + }, + { + "sample_id": "5460461103_R02C01", + "id": "HG00361", + "est": "M", + "ped": "F" + }, + { + "sample_id": "5885434025_R02C01", + "id": "HG01683", + "est": "F", + "ped": "M" + }, + { + "sample_id": "6298012108_R02C01", + "id": "HG02300", + "est": "F", + "ped": "M" + } + ], + "parent_offspring": { + "pairs_in_pedigree": 717, + "detected_first_degree": 717, + "kinship": { + "min": 0.19192453, + "p01": 0.2261421876, + "median": 0.24747136, + "p99": 0.25384499520000003, + "max": 0.25912801, + "n": 717 + }, + "missed": [] + }, + "close_pairs_total": 808, + "close_pairs_unexplained_by_pedigree": 69, + "unexplained_examples": [ + { + "sample_a": "5426511011_R02C01", + "sample_b": "5434246044_R01C01", + "n_sites": "99409", + "kinship": "0.09722701", + "concordance": "0.70339708", + "degree": "second", + "expected": "", + "id_a": "NA19664", + "id_b": "NA19660" + }, + { + "sample_a": "5426511011_R02C01", + "sample_b": "5434246104_R03C01", + "n_sites": "99304", + "kinship": "0.12977189", + "concordance": "0.71363691", + "degree": "second", + "expected": "", + "id_a": "NA19664", + "id_b": "NA19672" + }, + { + "sample_a": "5426529060_R04C01", + "sample_b": "5550917041_R04C01", + "n_sites": "99458", + "kinship": "0.12910780", + "concordance": "0.67726075", + "degree": "second", + "expected": "", + "id_a": "NA20334", + "id_b": "NA20337" + }, + { + "sample_a": "5426529068_R04C01", + "sample_b": "5557094087_R04C01", + "n_sites": "99667", + "kinship": "0.13348281", + "concordance": "0.67120511", + "degree": "second", + "expected": "", + "id_a": "NA20414", + "id_b": "NA19625" + }, + { + "sample_a": "5426529089_R04C01", + "sample_b": "5557094057_R02C01", + "n_sites": "99650", + "kinship": "0.10322544", + "concordance": "0.71299548", + "degree": "second", + "expected": "", + "id_a": "NA19679", + "id_b": "NA19677" + }, + { + "sample_a": "5426529120_R03C01", + "sample_b": "5434246104_R03C01", + "n_sites": "99603", + "kinship": "0.13687614", + "concordance": "0.71472747", + "degree": "second", + "expected": "", + "id_a": "NA19685", + "id_b": "NA19672" + }, + { + "sample_a": "5434078004_R02C01", + "sample_b": "5550917042_R02C01", + "n_sites": "99521", + "kinship": "0.12556672", + "concordance": "0.67380754", + "degree": "second", + "expected": "", + "id_a": "NA20332", + "id_b": "NA20343" + }, + { + "sample_a": "5434078008_R04C01", + "sample_b": "5460461105_R02C01", + "n_sites": "99521", + "kinship": "0.12256670", + "concordance": "0.73606575", + "degree": "second", + "expected": "", + "id_a": "HG00418", + "id_b": "HG00427" + }, + { + "sample_a": "5434078010_R02C01", + "sample_b": "5550950006_R02C01", + "n_sites": "99399", + "kinship": "0.14619155", + "concordance": "0.68763267", + "degree": "second", + "expected": "", + "id_a": "NA20282", + "id_b": "NA20302" + }, + { + "sample_a": "5434078010_R02C01", + "sample_b": "5826291134_R02C01", + "n_sites": "99144", + "kinship": "0.13752765", + "concordance": "0.68179618", + "degree": "second", + "expected": "", + "id_a": "NA20282", + "id_b": "NA20313" + }, + { + "sample_a": "5434078014_R03C01", + "sample_b": "5443767014_R03C01", + "n_sites": "99547", + "kinship": "0.14121405", + "concordance": "0.69652526", + "degree": "second", + "expected": "", + "id_a": "NA19382", + "id_b": "NA19380" + }, + { + "sample_a": "5434078014_R04C01", + "sample_b": "5557094086_R02C01", + "n_sites": "99525", + "kinship": "0.12724911", + "concordance": "0.67710625", + "degree": "second", + "expected": "", + "id_a": "NA20336", + "id_b": "NA20335" + }, + { + "sample_a": "5434078024_R02C01", + "sample_b": "5557086034_R01C01", + "n_sites": "99590", + "kinship": "0.12186217", + "concordance": "0.73817652", + "degree": "second", + "expected": "", + "id_a": "HG00584", + "id_b": "HG00579" + }, + { + "sample_a": "5434078027_R02C01", + "sample_b": "5550917017_R04C01", + "n_sites": "99408", + "kinship": "0.11196304", + "concordance": "0.73934693", + "degree": "second", + "expected": "", + "id_a": "HG00501", + "id_b": "HG00514" + }, + { + "sample_a": "5434078027_R02C01", + "sample_b": "5557086084_R04C01", + "n_sites": "99533", + "kinship": "0.11746296", + "concordance": "0.73728311", + "degree": "second", + "expected": "", + "id_a": "HG00501", + "id_b": "HG00526" + }, + { + "sample_a": "5434203040_R03C01", + "sample_b": "5568738044_R02C01", + "n_sites": "99242", + "kinship": "0.10730422", + "concordance": "0.73397352", + "degree": "second", + "expected": "", + "id_a": "HG00656", + "id_b": "HG00703" + }, + { + "sample_a": "5434203048_R01C01", + "sample_b": "5557107038_R03C01", + "n_sites": "99616", + "kinship": "0.12640012", + "concordance": "0.71804730", + "degree": "second", + "expected": "", + "id_a": "NA19657", + "id_b": "NA19753" + }, + { + "sample_a": "5434203052_R02C01", + "sample_b": "5528550084_R02C01", + "n_sites": "99487", + "kinship": "0.13330070", + "concordance": "0.69495512", + "degree": "second", + "expected": "", + "id_a": "NA19027", + "id_b": "NA19311" + }, + { + "sample_a": "5434203065_R04C01", + "sample_b": "5443767041_R02C01", + "n_sites": "99564", + "kinship": "0.10208690", + "concordance": "0.67972359", + "degree": "second", + "expected": "", + "id_a": "NA19434", + "id_b": "NA19453" + }, + { + "sample_a": "5434203066_R01C01", + "sample_b": "5568738044_R02C01", + "n_sites": "99209", + "kinship": "0.13288993", + "concordance": "0.74309790", + "degree": "second", + "expected": "", + "id_a": "HG00657", + "id_b": "HG00703" + }, + { + "sample_a": "5434203074_R03C01", + "sample_b": "5568738010_R02C01", + "n_sites": "99489", + "kinship": "0.12890681", + "concordance": "0.74493663", + "degree": "second", + "expected": "", + "id_a": "HG00577", + "id_b": "HG00585" + }, + { + "sample_a": "5434203099_R04C01", + "sample_b": "5557107035_R04C01", + "n_sites": "99515", + "kinship": "0.15078348", + "concordance": "0.72261468", + "degree": "second", + "expected": "", + "id_a": "NA19661", + "id_b": "NA19686" + }, + { + "sample_a": "5434203102_R02C01", + "sample_b": "5434246080_R03C01", + "n_sites": "99550", + "kinship": "0.27029949", + "concordance": "0.82480161", + "degree": "first", + "expected": "", + "id_a": "", + "id_b": "" + }, + { + "sample_a": "5434203102_R02C01", + "sample_b": "5460461110_R04C01", + "n_sites": "99646", + "kinship": "0.24725288", + "concordance": "0.81377075", + "degree": "first", + "expected": "", + "id_a": "", + "id_b": "HG00635" + }, + { + "sample_a": "5434203102_R02C01", + "sample_b": "5550888027_R04C01", + "n_sites": "99305", + "kinship": "0.11918687", + "concordance": "0.73763657", + "degree": "second", + "expected": "", + "id_a": "", + "id_b": "HG00636" + } + ], + "mendel": { + "trios": 335, + "error_rate": { + "min": 0.00013780579359228424, + "p01": 0.0001580795944988077, + "median": 0.00031072140415760394, + "p99": 0.0007011283263471403, + "max": 0.0010779448028232972, + "n": 335 + }, + "worst": [ + { + "child": "NA12329", + "father": "NA06984", + "mother": "NA06989", + "errors": 2559, + "informative": 2373962, + "rate": 0.0010779448028232972 + }, + { + "child": "HG01099", + "father": "HG01097", + "mother": "HG01098", + "errors": 1954, + "informative": 2349923, + "rate": 0.0008315166071398935 + }, + { + "child": "NA19258", + "father": "NA19256", + "mother": "NA19257", + "errors": 1951, + "informative": 2358051, + "rate": 0.0008273782034400444 + }, + { + "child": "NA19199", + "father": "NA19198", + "mother": "NA19197", + "errors": 1656, + "informative": 2357810, + "rate": 0.0007023466691548513 + }, + { + "child": "HG02480", + "father": "HG02478", + "mother": "HG02479", + "errors": 1644, + "informative": 2352728, + "rate": 0.0006987633079557008 + }, + { + "child": "NA19918", + "father": "NA19916", + "mother": "NA19917", + "errors": 1536, + "informative": 2365464, + "rate": 0.0006493440610383417 + }, + { + "child": "NA19778", + "father": "NA19777", + "mother": "NA19776", + "errors": 1439, + "informative": 2367621, + "rate": 0.0006077830869045341 + }, + { + "child": "HG02222", + "father": "HG02221", + "mother": "HG02220", + "errors": 1397, + "informative": 2358117, + "rate": 0.0005924218348792702 + }, + { + "child": "HG01463", + "father": "HG01461", + "mother": "HG01462", + "errors": 1389, + "informative": 2350850, + "rate": 0.0005908501180424102 + }, + { + "child": "NA20279", + "father": "NA20278", + "mother": "NA20282", + "errors": 1378, + "informative": 2368169, + "rate": 0.0005818841476262885 + } + ] + }, + "contamination": { + "min": -0.00131277, + "p01": -0.000799082, + "median": 0.00024026, + "p99": 0.009337409999999999, + "max": 0.07474245, + "n": 2141, + "above_0.02": 15 + }, + "contamination_top": [ + { + "sample_id": "5434246028_R02C01", + "id": "HG00707", + "c": 0.07474245 + }, + { + "sample_id": "5434246028_R04C01", + "id": "HG01072", + "c": 0.06748388 + }, + { + "sample_id": "5460479065_R04C01", + "id": "NA18982", + "c": 0.06604837 + }, + { + "sample_id": "5767788059_R02C01", + "id": "HG01667", + "c": 0.06352254 + }, + { + "sample_id": "5426529051_R04C01", + "id": "NA20760", + "c": 0.06220859 + }, + { + "sample_id": "5885434025_R01C01", + "id": "HG00098", + "c": 0.05900927 + }, + { + "sample_id": "5434246054_R02C01", + "id": "HG00335", + "c": 0.05507963 + }, + { + "sample_id": "5434246100_R04C01", + "id": "HG00637", + "c": 0.03604213 + }, + { + "sample_id": "5434203056_R03C01", + "id": "HG00358", + "c": 0.02990194 + }, + { + "sample_id": "5434246020_R01C01", + "id": "NA19129", + "c": 0.02691854 + }, + { + "sample_id": "5434246020_R03C01", + "id": "NA12234", + "c": 0.0254513 + }, + { + "sample_id": "5826338161_R02C01", + "id": "HG01686", + "c": 0.02254977 + }, + { + "sample_id": "5443767155_R04C01", + "id": "NA20526", + "c": 0.02212635 + }, + { + "sample_id": "5460479059_R04C01", + "id": "NA19058", + "c": 0.02205388 + }, + { + "sample_id": "6317383048_R03C01", + "id": "HG02229", + "c": 0.02015661 + } + ], + "functional_status": { + "FAIL": 604, + "PASS": 1537 + }, + "control_status": { + "FLAGGED": 1795, + "CLEAR": 346 + }, + "control_flag_codes": { + "stainingGreen": 1794, + "stainingRed": 1794 + }, + "flags": { + "high_lrr_sd": 602, + "unexpected_relatedness": 1139, + "contaminated": 15, + "low_call_rate": 2 + }, + "control_metrics": { + "restoration": {}, + "staining_green": { + "min": 39.41527002, + "p01": 44.6584822776, + "median": 95.23627288, + "p99": 715.4660557490015, + "max": 2566.66666667, + "n": 347 + }, + "staining_red": { + "min": 0.0, + "p01": 43.783606462600005, + "median": 106.485488455, + "p99": 202.30416815859996, + "max": 224.93461538, + "n": 348 + }, + "extension_green": { + "min": 14.32592078, + "p01": 17.28083184, + "median": 31.76146789, + "p99": 61.591069746, + "max": 84.0969163, + "n": 2141 + }, + "extension_red": { + "min": 10.56216216, + "p01": 14.684177382, + "median": 23.79825835, + "p99": 44.49232804199999, + "max": 49.67384615, + "n": 2141 + }, + "hyb_high_med": { + "min": 1.63464397, + "p01": 1.831392898, + "median": 2.10967896, + "p99": 2.4262993479999997, + "max": 2.64515327, + "n": 2141 + }, + "hyb_med_low": { + "min": 1.81739684, + "p01": 1.999406742, + "median": 2.32677346, + "p99": 2.697378338, + "max": 3.23596807, + "n": 2141 + }, + "target_removal_1": { + "min": 3.7875, + "p01": 4.624796872, + "median": 9.41374838, + "p99": 23.095572163999982, + "max": 39.10588235, + "n": 2141 + }, + "stringency": { + "min": 4.61498708, + "p01": 5.13971662, + "median": 9.30527069, + "p99": 19.436929915999976, + "max": 23.71089631, + "n": 2141 + }, + "nsb_green": { + "min": 3.72138047, + "p01": 4.183322406, + "median": 9.24747475, + "p99": 23.321340729999996, + "max": 39.625, + "n": 2141 + }, + "nsb_red": { + "min": 5.03819444, + "p01": 5.9731914900000005, + "median": 13.42666667, + "p99": 26.774626863999966, + "max": 38.24742268, + "n": 2141 + }, + "np_green": { + "min": 9.49642346, + "p01": 11.214550714, + "median": 22.8027027, + "p99": 41.890710461999994, + "max": 57.11235955, + "n": 2141 + }, + "np_red": { + "min": 6.56808511, + "p01": 8.242032115999999, + "median": 14.65871122, + "p99": 23.917212037999985, + "max": 26.11714286, + "n": 2141 + } + }, + "het_rate": { + "min": 0.14461898, + "p01": 0.156512408, + "median": 0.18633936, + "p99": 0.215229702, + "max": 0.22765567, + "n": 2141 + }, + "seconds": 291.3, + "checks": { + "identified_ge_99pct": true, + "sex_agreement_ge_99pct": true, + "parent_offspring_detected_ge_99pct": true, + "median_mendel_rate_lt_0.001": true + }, + "passed": true +} diff --git a/src/simurg/pipeline/qc.py b/src/simurg/pipeline/qc.py index f55430b..9ee54ab 100644 --- a/src/simurg/pipeline/qc.py +++ b/src/simurg/pipeline/qc.py @@ -57,6 +57,8 @@ def fmt(v: Any, digits: int = 8) -> str: return "" if isinstance(v, float | np.floating): f = float(v) + if math.isinf(f): + return "Inf" if f > 0 else "-Inf" return "NaN" if math.isnan(f) else f"{f:.{digits}f}" return str(v).replace("\t", " ").replace("\n", " ") diff --git a/src/simurg/pipeline/qc_report.py b/src/simurg/pipeline/qc_report.py index 913ff82..f1b6af9 100644 --- a/src/simurg/pipeline/qc_report.py +++ b/src/simurg/pipeline/qc_report.py @@ -76,6 +76,7 @@ FORMULAS = { CHIP_FORMATS = {4: (4, 1), 8: (8, 1), 12: (12, 1), 16: (8, 2), 24: (12, 2), 48: (12, 4), 96: (12, 8)} POS_RE = re.compile(r"^(\d+)_R(\d{2})C(\d{2})$") LABEL_RE = re.compile(r"[^A-Za-z0-9_.-]") +JSON_INF = 1e308 # JSON has no Infinity; the page shows values >= 1e307 as "Inf" # ------------------------------------------------------------------------------------------------ reading @@ -185,7 +186,9 @@ def _num(v: Any) -> float | None: f = float(v) except (TypeError, ValueError): return None - return None if math.isnan(f) or math.isinf(f) else f + if math.isinf(f): + return JSON_INF if f > 0 else -JSON_INF + return None if math.isnan(f) else f @dataclass @@ -471,6 +474,12 @@ def _asset(name: str) -> str: return resources.files("simurg.qc.report").joinpath(name).read_text(encoding="utf-8") +def _plotly() -> str: + """Bundled plotly.js 2.35.2 (MIT), stored gzip-compressed in the package.""" + data = resources.files("simurg.qc.report").joinpath("vendor", "plotly.min.js.gz").read_bytes() + return gzip.decompress(data).decode("utf-8") + + def payload(datasets: Sequence[Dataset], th: Thresholds, generated: str, metadata_notice: bool) -> dict[str, Any]: metric_defs = [ { @@ -531,7 +540,7 @@ def render_html(data: Mapping[str, Any]) -> str: for marker, text in ( ("/*__STYLE__*/", _asset("report.css")), ("/*__APP__*/", _asset("report.js")), - ("/*__PLOTLY__*/", _asset("vendor/plotly.min.js").replace(" Any: if isinstance(o, float): - return None if math.isnan(o) or math.isinf(o) else o + if math.isinf(o): + return JSON_INF if o > 0 else -JSON_INF + return None if math.isnan(o) else o if isinstance(o, dict): return {k: _sanitize_json(v) for k, v in o.items()} if isinstance(o, list): diff --git a/src/simurg/qc/controls.py b/src/simurg/qc/controls.py index dce681b..e879c15 100644 --- a/src/simurg/qc/controls.py +++ b/src/simurg/qc/controls.py @@ -192,12 +192,17 @@ class _Probes: self.absent: set[str] = set() def vals(self, ids: Sequence[str], ch: str) -> npt.NDArray[np.float64]: - rows: list[int] = [] + """Intensities of every bead of ``ids`` in channel ``ch``. Bead addresses absent from the IDAT (NaN; e.g. + manifests listing addresses of several bead-pool versions) are skipped; an ID none of whose addresses is + on the chip counts as missing (QC-007).""" + out: list[float] = [] for cid in ids: - if cid not in self.idx: + v = self.ch[ch][self.idx.get(cid, [])] + v = v[~np.isnan(v)] + if not len(v): self.absent.add(cid) - rows.extend(self.idx.get(cid, [])) - return self.ch[ch][rows] + out.extend(v.tolist()) + return np.asarray(out, dtype=np.float64) def mean(self, cid: str, ch: str) -> float: v = self.vals([cid], ch) @@ -249,9 +254,9 @@ def compute_control_metrics( if p.absent: res.missing[m.column] = sorted(p.absent) res.values[m.column] = math.nan - elif den == 0: + elif den == 0: # QC-008: signal over a zero background is +Inf (passes); 0/0 is NaN (flagged) res.zero_divisor.append(m.column) - res.values[m.column] = math.nan + res.values[m.column] = math.inf if num > 0 else math.nan else: - res.values[m.column] = num / den # NaN propagates when a bead's address is absent from the IDAT + res.values[m.column] = num / den return res diff --git a/src/simurg/qc/report/report.js b/src/simurg/qc/report/report.js index 1f9684c..e489fb2 100644 --- a/src/simurg/qc/report/report.js +++ b/src/simurg/qc/report/report.js @@ -66,7 +66,7 @@ options.map(([v, t]) => { const o = $("option", { value: v }, t); if (String(v) === String(value)) o.selected = true; return o; })); return $("label", {}, label, s); }; - const fmt = (v, d = 4) => (v === null || v === undefined || Number.isNaN(v) ? "" : typeof v === "number" ? +v.toFixed(d) : v); + const fmt = (v, d = 4) => (v === null || v === undefined || Number.isNaN(v) ? "" : typeof v === "number" ? (Math.abs(v) >= 1e307 ? (v > 0 ? "Inf" : "-Inf") : +v.toFixed(d)) : v); const ds = () => D.datasets[cur]; const val = (s, col) => (col in s.m ? s.m[col] : s[col] !== undefined ? s[col] : s.meta[col]); const trunc = (x) => { x = String(x ?? ""); return x.length > 200 ? x.slice(0, 200) + "…" : x; }; diff --git a/tests/unit/test_qc.py b/tests/unit/test_qc.py index 54226df..aa45e00 100644 --- a/tests/unit/test_qc.py +++ b/tests/unit/test_qc.py @@ -147,7 +147,24 @@ def test_control_missing_and_zero_divisor() -> None: # QC-007, QC-008 r[i] = 0 q = compute_control_metrics(beads, g, r) assert math.isnan(q.values["restoration"]) and q.missing["restoration"] == ["restoration"] - assert q.zero_divisor == ["stringency"] and math.isnan(q.values["stringency"]) + assert q.zero_divisor == ["stringency"] and q.values["stringency"] == math.inf # QC-008: signal / 0 passes + r[i - 1] = 0 + q = compute_control_metrics(beads, g, r) + assert math.isnan(q.values["stringency"]) # 0 / 0 is NaN (flagged) + + +def test_absent_addresses_are_skipped() -> None: # QC-008, K-021: two bead-pool addresses, one on the chip + beads = beads_from_definitions( + [ + {"address": "5:105", "category": "Staining", "color": "Green", "name": "Biotin (High)"}, + {"address": "6:106", "category": "Staining", "color": "Blue", "name": "Biotin (Bgnd)"}, + ] + ) + g = np.array([np.nan, 900.0, 100.0, np.nan]) + q = compute_control_metrics(beads, g, np.zeros(4)) + assert q.values["staining_green"] == pytest.approx(9.0) + q = compute_control_metrics(beads, np.array([np.nan, 900.0, np.nan, np.nan]), np.zeros(4)) + assert q.missing["staining_green"] == ["biotin_bkg"] def test_methylation_is_v2() -> None: