v2 methylation: sesame oracle scripts and cloud task
Co-Authored-By: Claude Opus 5.5 <[email protected]> Claude-Session: https://claude.ai/code/session_01814GKVzZVL3FuF3DeiZXhv
This commit is contained in:
kemalandClaude Opus 5.5 committed
1 parent
c1c2c30c87
commit
2c503d2742
3 files changed
+144
No files matched your search
@@ -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: <unused data dir> <results dir>. 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) ==="
|
||||
@@ -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())
|
||||
@@ -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()
|
||||
Reference in new issue
Block a user