# 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)]) {
  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)
  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()