From da4a67b59b06490a45aa4d4a8dfa36c85a8b97d1 Mon Sep 17 00:00:00 2001 From: Kemal Yaylali Date: Tue, 29 Sep 2026 00:14:18 +0100 Subject: [PATCH] v2 methylation: sesame oracle scripts and cloud task --- scripts/cloud/tasks/meth.sh | 33 +++++++++++++ scripts/v2/methyl_oracle_compare.py | 75 +++++++++++++++++++++++++++++ scripts/v2/sesame_oracle.R | 36 ++++++++++++++ 3 files changed, 144 insertions(+) create mode 100644 scripts/cloud/tasks/meth.sh create mode 100644 scripts/v2/methyl_oracle_compare.py create mode 100644 scripts/v2/sesame_oracle.R diff --git a/scripts/cloud/tasks/meth.sh b/scripts/cloud/tasks/meth.sh new file mode 100644 index 0000000..c5c7fb1 --- /dev/null +++ b/scripts/cloud/tasks/meth.sh @@ -0,0 +1,33 @@ +#!/bin/bash +# v2 methylation QC: numerics vs the sesame oracle (R, MIT) on the public minfiDataEPIC set (3 x NA12878, EPIC v1). +# Args: . Run with OMNI=0. +set -uo pipefail +OUT=$2 +R=$OUT/meth; W=/work/meth; mkdir -p $R $W +ARCH=$(uname -m) +echo "=== data $(date -u +%T) ===" +cd $W +curl -sSfL -O https://bioconductor.org/packages/release/data/experiment/src/contrib/minfiDataEPIC_1.38.0.tar.gz \ + && tar -xzf minfiDataEPIC_1.38.0.tar.gz minfiDataEPIC/inst/extdata +curl -sSfL -O https://webdata.illumina.com/downloads/productfiles/methylationEPIC/infinium-methylationepic-v-1-0-b5-manifest-file-csv.zip \ + && unzip -oq infinium-methylationepic-v-1-0-b5-manifest-file-csv.zip +cd - > /dev/null +MAN=$W/infinium-methylationepic-v-1-0-b5-manifest-file.csv +E=$W/minfiDataEPIC/inst/extdata +PREFIXES="$E/200144450018/200144450018_R04C01 $E/200144450019/200144450019_R07C01 $E/200144450021/200144450021_R05C01" +.venv/bin/simurg methyl convert-manifest --manifest $MAN --out $W/mft --quiet + +echo "=== sesame (conda) $(date -u +%T) ===" +curl -sSfL -o $W/mf.sh "https://github.com/conda-forge/miniforge/releases/latest/download/Miniforge3-Linux-$ARCH.sh" \ + && bash $W/mf.sh -b -p $W/conda > /dev/null +$W/conda/bin/mamba create -y -q -p $W/renv -c conda-forge -c bioconda r-base bioconductor-sesame > $R/conda.log 2>&1 \ + || { echo "!! sesame install failed"; tail -30 $R/conda.log; } +$W/renv/bin/Rscript scripts/v2/sesame_oracle.R $W/mft/probes.tsv $W/mft/controls.tsv $W/sesame $PREFIXES \ + > $R/sesame.log 2>&1 || echo "!! sesame run failed" +tail -40 $R/sesame.log + +echo "=== compare $(date -u +%T) ===" +.venv/bin/python scripts/v2/methyl_oracle_compare.py --manifest $MAN --sesame $W/sesame --out $R/compare.json $PREFIXES +.venv/bin/simurg methyl qc --manifest $MAN --idats $E --sheet $E/SampleSheet.csv --out $W/qc --threads 0 --quiet +cp $W/qc/methyl_qc_summary.tsv $W/qc/controls_qc.tsv $R/ +echo "=== task done $(date -u +%T) ===" diff --git a/scripts/v2/methyl_oracle_compare.py b/scripts/v2/methyl_oracle_compare.py new file mode 100644 index 0000000..f9bdddb --- /dev/null +++ b/scripts/v2/methyl_oracle_compare.py @@ -0,0 +1,75 @@ +#!/usr/bin/env python3 +"""Compare `simurg methyl qc` numerics with the sesame oracle (scripts/v2/sesame_oracle.R) on the same probes. + +python scripts/v2/methyl_oracle_compare.py --manifest M.csv --sesame DIR --out summary.json PREFIX... +""" + +from __future__ import annotations + +import argparse +import csv +import json +from pathlib import Path + +import numpy as np + +from simurg.formats.idat import read_idat +from simurg.formats.methyl_manifest import read_methyl_manifest +from simurg.signal import methyl as sig + + +def _q(a: np.ndarray) -> dict[str, float]: + a = a[np.isfinite(a)] + if not len(a): + return {} + return {k: float(np.quantile(a, q)) for k, q in (("median", 0.5), ("p99", 0.99), ("max", 1.0))} + + +def main() -> int: + ap = argparse.ArgumentParser() + ap.add_argument("--manifest", required=True) + ap.add_argument("--sesame", required=True, type=Path) + ap.add_argument("--out", required=True, type=Path) + ap.add_argument("prefixes", nargs="+") + a = ap.parse_args() + man = read_methyl_manifest(a.manifest) + pos = {p: i for i, p in enumerate(man.probe_ids)} + res = {} + for prefix in a.prefixes: + s = sig.assemble(read_idat(prefix + "_Grn.idat"), read_idat(prefix + "_Red.idat"), man) + p = sig.detection_p(s) + n, fg, fr = sig.noob(s) + b = sig.beta_values(n) + rows = list(csv.DictReader((a.sesame / (Path(prefix).name + ".tsv")).open(), delimiter="\t")) + idx = np.array([pos[r["probe_id"]] for r in rows]) + + def col(name: str, rows: list[dict[str, str]] = rows) -> np.ndarray: + return np.array([float(r[name]) if r[name] not in ("NA", "") else np.nan for r in rows]) + + mask_s = np.array([r["mask"] == "TRUE" for r in rows]) + keep = ~mask_s & ~s.masked[idx] + sp, sm, su, sb = col("pval"), col("m_noob"), col("u_noob"), col("beta_noob") + op, om, ou, ob = p[idx], n.m[idx], n.u[idx], b[idx] + ana = man.analytical()[idx] & keep + res[Path(prefix).name] = { + "probes": len(rows), + "masked_sesame": int(mask_s.sum()), + "masked_simurg": int(s.masked.sum()), + "mask_agreement": float((mask_s == s.masked[idx]).mean()), + "noob_fit": {"green": vars(fg), "red": vars(fr)}, + "pval_abs_diff": _q(np.abs(op - sp)[keep]), + "pval_identical_fraction": float((op[keep] == sp[keep]).mean()), + "detection_call_agreement_0.05": float(((op < 0.05) == (sp < 0.05))[keep].mean()), + "prop_passing": {"simurg": float((op[ana] < 0.05).mean()), "sesame": float((sp[ana] < 0.05).mean())}, + "noob_m_rel_diff": _q(np.abs(om - sm)[keep] / np.maximum(np.abs(sm[keep]), 1)), + "noob_u_rel_diff": _q(np.abs(ou - su)[keep] / np.maximum(np.abs(su[keep]), 1)), + "beta_abs_diff": _q(np.abs(ob - sb)[keep]), + "beta_corr": float(np.corrcoef(ob[keep & np.isfinite(sb)], sb[keep & np.isfinite(sb)])[0, 1]), + } + a.out.write_text(json.dumps(res, indent=1)) + print(json.dumps(res, indent=1)) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/scripts/v2/sesame_oracle.R b/scripts/v2/sesame_oracle.R new file mode 100644 index 0000000..5be5be9 --- /dev/null +++ b/scripts/v2/sesame_oracle.R @@ -0,0 +1,36 @@ +# Oracle for `simurg methyl qc` numerics: sesame (MIT) pOOBAH, noob and beta values on the same probe set. +# Usage: Rscript sesame_oracle.R probes.tsv controls.tsv OUTDIR IDAT_PREFIX... +# probes.tsv / controls.tsv are written by `simurg methyl convert-manifest` (MET-003: sesame manifest layout). +suppressPackageStartupMessages(library(sesame)) +args <- commandArgs(trailingOnly = TRUE) +p <- read.delim(args[1], stringsAsFactors = FALSE, na.strings = "NA", colClasses = c("character", "numeric", + "numeric", "character")) +mft <- data.frame(Probe_ID = p$probe_id, M = p$m, U = p$u, col = ifelse(is.na(p$col), "2", p$col), + stringsAsFactors = FALSE) +mft$col <- factor(mft$col, levels = c("G", "R", "2")) +ctl <- read.delim(args[2], stringsAsFactors = FALSE) +controls <- data.frame(Address = as.character(ctl$address), Type = ctl$type, Color_Channel = ctl$color_channel, + Name = ctl$name, stringsAsFactors = FALSE) +out <- args[3] +dir.create(out, showWarnings = FALSE, recursive = TRUE) +cat("sesame", as.character(packageVersion("sesame")), "\n") +inband <- function(s) { + m <- ifelse(s$col == "R", s$MR, ifelse(s$col == "G", s$MG, s$UG)) + u <- ifelse(s$col == "R", s$UR, ifelse(s$col == "G", s$UG, s$UR)) + 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) }) + pv <- pOOBAH(sdf, return.pval = TRUE) + sn <- noob(sdf) + ib <- inband(sn) + b <- getBetas(sn, mask = FALSE) + res <- data.frame(probe_id = sdf$Probe_ID, mask = sdf$mask, pval = pv[sdf$Probe_ID], + m_noob = ib$m, u_noob = ib$u, beta_noob = b[sdf$Probe_ID]) + f <- file.path(out, paste0(basename(prefix), ".tsv")) + write.table(res, f, sep = "\t", quote = FALSE, row.names = FALSE) + cat(basename(prefix), nrow(res), "probes;", sum(res$mask), "masked\n") +} +sessionInfo()