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:
1 parent
b8695c835e
commit
579009271e
5 files changed
+326
-2
No files matched your search
@@ -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 |
|
||||
|
||||
@@ -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
|
||||
}
|
||||
}
|
||||
}
|
||||
Reference in new issue
Block a user