37 lines
2.0 KiB
R
37 lines
2.0 KiB
R
# 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()
|