From cb3f5fac5819d76a7b9ee683ae45360c54a99151 Mon Sep 17 00:00:00 2001 From: Kemal Yaylali Date: Tue, 29 Sep 2026 00:41:54 +0100 Subject: [PATCH] v2 methylation QC validated on minfiDataEPIC and against sesame (K-026) beta r 0.99998 vs sesame 1.28.1 (max |dbeta| 0.011), detection calls 99.1% identical, noob within 0.5% median; replicate r 0.997; all 21 control metrics pass. Spec deviations recorded (5-field control lines, pOOBAH formula, minfi sheet keys). User guide section. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01814GKVzZVL3FuF3DeiZXhv --- CHANGELOG.md | 7 + README.md | 4 +- docs/known_unknowns.md | 1 + docs/validation/V2_methylation.md | 53 ++++ docs/validation/V2_methylation_vs_sesame.json | 263 ++++++++++++++++++ 5 files changed, 326 insertions(+), 2 deletions(-) create mode 100644 docs/validation/V2_methylation.md create mode 100644 docs/validation/V2_methylation_vs_sesame.json diff --git a/CHANGELOG.md b/CHANGELOG.md index 890d1dc..7f27758 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,12 @@ # Changelog +## Unreleased (v2) + +**Methylation control QC** (`simurg methyl qc`, `simurg methyl convert-manifest`; MET-001 … MET-018): methylation +CSV manifests, the 21 control metrics, pOOBAH detection p-values, noob, linear/nonlinear dye bias, β and M-values, +summary tsv/xlsx/html with a β PCA, per-sample tables; methylation folders load in `simurg qc report`; +`qc controls --array-type methylation`. Checked against sesame (`docs/validation/V2_methylation.md`, K-026). + ## 1.1.0 — 2026-09-28 (P9, P10) **Packaging (P10).** `scripts/build_dist.sh` (sdist, pure wheel, manylinux_2_28 kernel wheels for x86_64 and diff --git a/README.md b/README.md index c9286f1..ed5bd63 100644 --- a/README.md +++ b/README.md @@ -24,8 +24,8 @@ Commands: `simurg inspect`, `ingest`, `normalize`, `genotype`, `qc controls|samp `gtc-to-vcf`, `export plink|bedgraph|finalreport|locus-summary`, `compare`. **User guide: `docs/user_guide.md`**; worked example: `examples/gsa_demo.sh`. P9 (performance) is done: optional Rust kernels, byte-identical outputs, ~80 GSA samples/min per server core (`docs/validation/P9.md`). P10 (packaging) is done: wheels, Docker image and -conda recipes, verified on fresh x86_64 and aarch64 machines (`docs/validation/P10.md`). Next (plan §18): v2 modules -(cluster training, CN/LOH, PGx). +conda recipes, verified on fresh x86_64 and aarch64 machines (`docs/validation/P10.md`). v2 in progress: methylation control QC (`simurg methyl qc`, +`docs/validation/V2_methylation.md`); next CN/LOH, PGx, cluster training. ## Quick start diff --git a/docs/known_unknowns.md b/docs/known_unknowns.md index c69450b..6f2d917 100644 --- a/docs/known_unknowns.md +++ b/docs/known_unknowns.md @@ -4,6 +4,7 @@ Plan §20 lists the open questions; this file records what the implementation ha | ID | Date | Topic | Finding | Status | |---|---|---|---|---| +| K-026 | 2026-09-29 | Methylation QC (v2) | Spec deviations: vendor EPIC manifests put a 5th field on some control lines (≥ 4 fields accepted); detection p follows sesame's pOOBAH (1 − max(F(M), F(U))) rather than the ambiguous "1 − ECDF(M + U)"; minfi-style `Sentrix_ID`/`Sentrix_Position` sheet keys accepted. Against sesame 1.28.1 on minfiDataEPIC: β r = 0.99998 (max \|Δβ\| 0.011), detection calls 99.1 % identical, noob signals within 0.5 % (median) / 5 % (max). The systematic remainder is unexplained; it is not the negative controls (sesame did not attach the supplied control table). | Open; `docs/validation/V2_methylation.md`. Next step: dump sesame's background vectors and noob parameters to find the difference. | | K-025 | 2026-09-28 | Scaling limits (P9) | With the Rust kernels a Graviton3 core genotypes 78–80 GSA samples/min. At > 8 workers a 125 MB/s gp3 volume caps the run at ~620 samples/min (GTC ≈ 13 MB each); on a RAM disk efficiency is 0.77 at 16 and 0.48 at 32 workers for 360 samples, explained by a fixed start-up cost (≈ 2.5 s + 0.2 s per worker). Linear scaling for large batches is inferred from that model, not measured. | `docs/validation/P9.md`; measure a ≥ 5000-sample batch when one is at hand. | | K-024 | 2026-09-28 | GSE165845 (P8) | GSA-MD v1 chips have no public cluster file; the GSA v1 C1 EGT with the C2 build-38 BPM covers 99.996 % of bead addresses and gives median LogR deviation 0.246 (cluster mismatch, not noise). GSM5051309 (NA20362 rep 2) is a failed array (call rate 0.762, concordance 0.336 with its own replicates). 27 510 array records have no NYGC 30x record and are scored as hom-ref; 1444 of them are systematically called alt. | Recorded in `docs/validation/P8.md`; results reported with and without these. | | K-023 | 2026-09-28 | P6 optional outputs | Implemented and validated the same day: `--info-fields` (with `--egt`), `--extra-format`, `--merge-snv-indel`, `--auxiliary-loci`, PLINK 2 `.pgen`. Two decisions: OUT-035 rule (d) uses the MNV's last plus-strand base for minus-strand designs (spec amended); `--extra-format` writes NORMX/NORMY missing where the raw intensities are both 0, and THETA/R missing where both normalized intensities are 0 (148 values on 3 GSA samples; gtc2vcf writes 0). | resolved | diff --git a/docs/validation/V2_methylation.md b/docs/validation/V2_methylation.md new file mode 100644 index 0000000..d685c6b --- /dev/null +++ b/docs/validation/V2_methylation.md @@ -0,0 +1,53 @@ +# v2 module validation — methylation control QC (`simurg methyl qc`) + +Date: 2026-09-29. Development on the Raspberry Pi; oracle runs on AWS (c7g.xlarge, five short runs, ≈ $0.025). +Spec: MET-001 … MET-018 (and QC-004/QC-006, QC-062). + +## Data + +* **minfiDataEPIC** (Bioconductor, Artistic-2.0): 3 EPIC v1 arrays of NA12878 (IDATs + minfi-style sample sheet). +* **Infinium MethylationEPIC v1.0 B5 manifest (CSV)**: 865 918 probes (142 158 Infinium I: 92 212 red, 49 946 green), + 635 control beads. Read in 8.6 s on the Pi. + +## Delivered + +`simurg methyl qc`, `simurg methyl convert-manifest`, `simurg qc controls --array-type methylation`, methylation +datasets in `simurg qc report`; modules `formats/methyl_manifest.py`, `signal/methyl.py`, `pipeline/methyl.py`; +the 21 methylation metrics in `qc/controls.py`. User guide: "Methylation arrays". + +## Results on minfiDataEPIC + +| Check | Result | +|---|---| +| Run time, 3 samples, 3 workers (Pi) | 25 s (2.9 s signal processing per sample) | +| 21 control metrics | all evaluable, all pass the `methylation_epic` profile (e.g. staining green 57, BS conversion I green 13, BS conversion II 5.3, specificity II 24) | +| Probes passing detection (p < 0.05, analytical CG + SNP probes) | 90.9 %, 91.1 %, 91.3 % | +| β of the three NA12878 replicates | Pearson r 0.997 – 0.9974 between replicates; bimodal distribution | +| Sheet `Sample_Group` → PCA colours; outputs of MET-013 | present; folder loads in `qc report` | + +## Against the sesame oracle (sesame 1.28.1, same probe set via `probes.tsv`/`controls.tsv`) + +`scripts/v2/sesame_oracle.R` + `scripts/v2/methyl_oracle_compare.py`, bead-count masking off on both sides (sesame's +`min_beads` masks every Infinium II probe because their M address is NA), 865 918 probes per sample: + +| Quantity | Agreement (3 samples) | +|---|---| +| β values | Pearson r 0.99998; \|Δβ\| median 0.0009, 99th percentile 0.008, max 0.011 | +| noob-corrected M/U signals | relative difference median 0.5 %, max 5 % | +| Detection calls at p < 0.05 | 99.0 – 99.1 % identical; passing fraction 91.0 – 91.4 % (simurg) vs 92.0 – 92.3 % (sesame) | +| Detection p-values | not identical (median \|Δp\| 0.0008) | + +The remaining difference is systematic (same size in all three samples) and is not the negative controls: sesame did +not attach the supplied control table (`negControls` failed, so its background is out-of-band Infinium I signal only), +and simurg without negative controls is not closer (noob difference 0.56 % instead of 0.53 %). The exact cause is +open (K-026). For QC use — pass/fail per sample, detection rate, β for a PCA — the agreement is sufficient; simurg's +methylation numerics are documented as its own [OURS] implementation of the published methods (MET-010), not as a +bit-exact port. + +## Spec deviations (K-026) + +* MET-002: the vendor EPIC manifest has a 5th field (`AVG`) on some `NEGATIVE` control lines; ≥ 4 fields are + accepted and the first four used. Fewer than 4 still gives `E_METH_CONTROLS_FORMAT`. +* MET-010: the detection p-value follows sesame's pOOBAH (p = 1 − max(F(M), F(U)), each in its own channel) instead + of the ambiguous "1 − ECDF(M + U)", so an open-source oracle exists. +* MET-004: `Sentrix_ID`/`Sentrix_Position` (minfi/sesame sheets) accepted as a third key-column pair. diff --git a/docs/validation/V2_methylation_vs_sesame.json b/docs/validation/V2_methylation_vs_sesame.json new file mode 100644 index 0000000..8aa2cf7 --- /dev/null +++ b/docs/validation/V2_methylation_vs_sesame.json @@ -0,0 +1,263 @@ +{ + "200144450018_R04C01": { + "probes": 865918, + "masked_sesame": 0, + "masked_simurg": 0, + "combine_neg=True": { + "noob_fit": { + "green": { + "mu": 883.7232985218219, + "sigma": 533.736, + "alpha": 6615.19430861235 + }, + "red": { + "mu": 1475.9648221678526, + "sigma": 974.0681999999999, + "alpha": 7235.1788325587695 + } + }, + "pval_abs_diff": { + "median": 0.0007849171113696351, + "p99": 0.012655374352560589, + "max": 0.7810294625483152 + }, + "pval_identical_fraction": 2.7716250268501175e-05, + "detection_call_agreement_0.05": 0.9904228806884716, + "prop_passing": { + "simurg": 0.9104377127786546, + "sesame": 0.9199708917641769 + }, + "noob_m_rel_diff": { + "median": 0.005309072922507809, + "p99": 0.03942181326074589, + "max": 0.04794913709182962 + }, + "noob_u_rel_diff": { + "median": 0.005799391773898118, + "p99": 0.03459884716823513, + "max": 0.04794913709182962 + }, + "beta_abs_diff": { + "median": 0.0008897385559986926, + "p99": 0.008521359578262878, + "max": 0.011219185321570957 + }, + "beta_corr": 0.9999813374332573 + }, + "combine_neg=False": { + "noob_fit": { + "green": { + "mu": 885.3040134669028, + "sigma": 533.736, + "alpha": 6613.613593667269 + }, + "red": { + "mu": 1481.7310124080875, + "sigma": 975.5508, + "alpha": 7229.412642318534 + } + }, + "pval_abs_diff": { + "median": 0.0007957874497818184, + "p99": 0.013082543985806876, + "max": 0.7841804034247452 + }, + "pval_identical_fraction": 2.7716250268501175e-05, + "detection_call_agreement_0.05": 0.9901872925611894, + "prop_passing": { + "simurg": 0.9102024830066768, + "sesame": 0.9199708917641769 + }, + "noob_m_rel_diff": { + "median": 0.005602741250468986, + "p99": 0.038392170387569056, + "max": 0.047002982560608926 + }, + "noob_u_rel_diff": { + "median": 0.0065460727630588, + "p99": 0.03442105845709709, + "max": 0.047002982560608926 + }, + "beta_abs_diff": { + "median": 0.0009450303433635281, + "p99": 0.008636600261016411, + "max": 0.011550874711670933 + }, + "beta_corr": 0.9999799503820251 + } + }, + "200144450019_R07C01": { + "probes": 865918, + "masked_sesame": 0, + "masked_simurg": 0, + "combine_neg=True": { + "noob_fit": { + "green": { + "mu": 1004.7491370858492, + "sigma": 606.3833999999999, + "alpha": 7726.238894681398 + }, + "red": { + "mu": 1645.4509315964758, + "sigma": 1097.124, + "alpha": 8498.773085184836 + } + }, + "pval_abs_diff": { + "median": 0.000796427075933526, + "p99": 0.012417625677686978, + "max": 0.8027931613655367 + }, + "pval_identical_fraction": 5.65873442981899e-05, + "detection_call_agreement_0.05": 0.9906884947535448, + "prop_passing": { + "simurg": 0.912474825779329, + "sesame": 0.9217102015559928 + }, + "noob_m_rel_diff": { + "median": 0.005314312414259365, + "p99": 0.04273539399348046, + "max": 0.051907861845576456 + }, + "noob_u_rel_diff": { + "median": 0.005468488622733155, + "p99": 0.03355605525549436, + "max": 0.051907861845576456 + }, + "beta_abs_diff": { + "median": 0.0008606416392141769, + "p99": 0.008394710099924632, + "max": 0.010860663350041877 + }, + "beta_corr": 0.9999812358896786 + }, + "combine_neg=False": { + "noob_fit": { + "green": { + "mu": 1006.5499088107555, + "sigma": 606.3833999999999, + "alpha": 7724.4381229564915 + }, + "red": { + "mu": 1651.850010672488, + "sigma": 1098.6065999999998, + "alpha": 8492.374006108823 + } + }, + "pval_abs_diff": { + "median": 0.0008019741538315371, + "p99": 0.012866482361703901, + "max": 0.8059396632834097 + }, + "pval_identical_fraction": 5.65873442981899e-05, + "detection_call_agreement_0.05": 0.9904552163137849, + "prop_passing": { + "simurg": 0.9122453898440994, + "sesame": 0.9217102015559928 + }, + "noob_m_rel_diff": { + "median": 0.005610084745351353, + "p99": 0.04163859014373106, + "max": 0.05094825885592487 + }, + "noob_u_rel_diff": { + "median": 0.006166215991548127, + "p99": 0.033209584923356225, + "max": 0.05094825885592487 + }, + "beta_abs_diff": { + "median": 0.0009157086256231864, + "p99": 0.008483978716612176, + "max": 0.011139138206961108 + }, + "beta_corr": 0.999980025075608 + } + }, + "200144450021_R05C01": { + "probes": 865918, + "masked_sesame": 0, + "masked_simurg": 0, + "combine_neg=True": { + "noob_fit": { + "green": { + "mu": 1075.2185580103821, + "sigma": 618.2442, + "alpha": 7644.388814794943 + }, + "red": { + "mu": 1583.358547469605, + "sigma": 1043.7504, + "alpha": 7872.027822306209 + } + }, + "pval_abs_diff": { + "median": 0.0008097160213839151, + "p99": 0.012557480308668073, + "max": 0.7179490998567224 + }, + "pval_identical_fraction": 8.083906328312843e-05, + "detection_call_agreement_0.05": 0.9906931141285895, + "prop_passing": { + "simurg": 0.9142593274977809, + "sesame": 0.9234993383438433 + }, + "noob_m_rel_diff": { + "median": 0.005312767934935867, + "p99": 0.0416027184427129, + "max": 0.05163890421954045 + }, + "noob_u_rel_diff": { + "median": 0.005627520758233918, + "p99": 0.033403350351144646, + "max": 0.05163890421954045 + }, + "beta_abs_diff": { + "median": 0.000886446134372592, + "p99": 0.008067365546443757, + "max": 0.01057770765422894 + }, + "beta_corr": 0.9999820250823297 + }, + "combine_neg=False": { + "noob_fit": { + "green": { + "mu": 1077.1430682973066, + "sigma": 618.2442, + "alpha": 7642.464304508018 + }, + "red": { + "mu": 1589.590396661226, + "sigma": 1045.233, + "alpha": 7865.7959731145875 + } + }, + "pval_abs_diff": { + "median": 0.0008198122345790535, + "p99": 0.013020971906371236, + "max": 0.7195549302039842 + }, + "pval_identical_fraction": 8.083906328312843e-05, + "detection_call_agreement_0.05": 0.9905060294392771, + "prop_passing": { + "simurg": 0.9140716071871386, + "sesame": 0.9234993383438433 + }, + "noob_m_rel_diff": { + "median": 0.005623618294636237, + "p99": 0.04045226512957654, + "max": 0.05065151257916101 + }, + "noob_u_rel_diff": { + "median": 0.006360115250782237, + "p99": 0.03307310255358611, + "max": 0.05065151257916101 + }, + "beta_abs_diff": { + "median": 0.0009431797599571898, + "p99": 0.008174273334787542, + "max": 0.010885211471086453 + }, + "beta_corr": 0.9999807797628828 + } + } +} \ No newline at end of file