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.
This commit is contained in:
kemal committed 2026-09-29 00:41:54 +01:00
1 parent 6aebb31991
commit 5164124438
5 files changed
+326 -2

No files matched your search

+1
View File
@@ -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 |
+53
View File
@@ -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.
@@ -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
}
}
}