From a9d2dfbac05dab418a2be4743258678b1fc218a7 Mon Sep 17 00:00:00 2001 From: Kemal Yaylali Date: Tue, 29 Sep 2026 00:18:17 +0100 Subject: [PATCH] sesame oracle: explicit platform, cache fallback --- docs/user_guide.md | 28 +++++++++++++++++++++++++++- scripts/v2/sesame_oracle.R | 9 ++++++--- 2 files changed, 33 insertions(+), 4 deletions(-) diff --git a/docs/user_guide.md b/docs/user_guide.md index 0cfec6e..9dbf8a1 100644 --- a/docs/user_guide.md +++ b/docs/user_guide.md @@ -113,6 +113,32 @@ simurg export finalreport --bpm M.bpm --csv M.csv --egt C.egt --gtcs run/gtc --o Global logging options on every command: `--quiet`, `--log-level`, `--log-json` (JSON Lines on stderr), `--debug`. Run `simurg COMMAND --help` for the full list. +## Methylation arrays + +`simurg methyl qc` runs control-probe QC and a signal summary for Infinium methylation arrays (EPIC v1/v2, MSA, +custom designs) from IDATs and the methylation CSV manifest (spec MET-001 … MET-018). It is a QC tool: there is no +differential-methylation analysis. + +```bash +simurg methyl qc --manifest infinium-methylationepic-v-1-0-b5-manifest-file.csv --idats idats/ \ + --sheet SampleSheet.csv --out run/methyl --profile auto +``` + +* 21 control metrics (staining, extension, hybridization, target removal I/II, bisulfite conversion I/II with their + background ratios, specificity I/II, non-polymorphic, restoration); thresholds from the `methylation_epic` / + `methylation_msa` profiles (`auto` picks by manifest name), `--config` or `--threshold KEY=VALUE`. +* Per sample: detection p-values (pOOBAH: out-of-band Infinium I signal and negative controls), noob background + correction, dye-bias correction (`--dye-bias linear|nonlinear`), β = M/(M+U+α) (`--beta-offset`, default 0) and + M-values. +* `pass_qc` = every enabled, evaluable metric passes; `--min-prop-passing X` also requires that fraction of probes + to pass detection (`--pval-threshold`, default 0.05). +* Outputs: `methyl_qc_summary.tsv/.xlsx/.html` (detection histogram and a β PCA coloured by `Sample_Group`), + `samples//_cgs.tsv.gz` (probe_id, beta, m_value, detection_p), `_ctrl.tsv.gz`, `pcs.tsv.gz`, + `probes.tsv`/`controls.tsv` (the manifest in the open-source toolkits' layout; also written by + `simurg methyl convert-manifest`), and `controls_raw.tsv`/`controls_qc.tsv`, so `simurg qc report` accepts the folder. +* Sample sheets may use `Sentrix_ID`/`Sentrix_Position` (minfi/sesame convention) as key columns. +* `simurg qc controls --array-type methylation --csv MANIFEST.csv` computes only the control tables. + ## Presets and reproducibility * `--preset current` (default): normalization 1.2.0, sex rule 2, GTC v5 — the current vendor software behaviour. @@ -185,7 +211,7 @@ A nightly differential run (`scripts/nightly.py`) re-checks these against the or * Speed: with the compiled kernels a GSA sample takes ~1.1 s on a Raspberry Pi 5 core (pure Python ~3.4 s); see `docs/validation/P9.md` for server numbers and thread scaling. Loading a manifest/cluster file is a fixed cost per run (GSA ~5 s, Omni2.5 ~15 s on the Pi). -* Methylation arrays: control-probe QC only. +* Methylation arrays: control-probe QC, detection and β values only (no differential analysis). ## Troubleshooting diff --git a/scripts/v2/sesame_oracle.R b/scripts/v2/sesame_oracle.R index 5be5be9..2693a3b 100644 --- a/scripts/v2/sesame_oracle.R +++ b/scripts/v2/sesame_oracle.R @@ -20,9 +20,12 @@ inband <- function(s) { list(m = m, u = u) } for (prefix in args[-(1:3)]) { - sdf <- tryCatch(readIDATpair(prefix, manifest = mft, controls = controls, min_beads = 3), - error = function(e) { message("min_beads not accepted: ", conditionMessage(e)); - readIDATpair(prefix, manifest = mft, controls = controls) }) + read <- function() readIDATpair(prefix, platform = "EPIC", manifest = mft, controls = controls, min_beads = 3) + sdf <- tryCatch(read(), error = function(e) { + message("first read failed (", conditionMessage(e), "); caching idatSignature and retrying") + sesameDataCache("idatSignature") + read() + }) pv <- pOOBAH(sdf, return.pval = TRUE) sn <- noob(sdf) ib <- inband(sn)