.Rhistory

mz <<- unlist(mapply(mzl, oidx, FUN = function(y, z) {
return(y[z])
}, SIMPLIFY = FALSE, USE.NAMES = FALSE), use.names = FALSE)
int <<- unlist(mapply(as.list(split(int, f = splitF)), oidx,
FUN=function(y, z) {
return(y[z])
}, SIMPLIFY = FALSE, USE.NAMES = FALSE),
use.names = FALSE)
rm(mzl)
rm(splitF)
tmp <- capture.output(
roiList <<- .Call("findmzROI",
mz, int, scanindex,
as.double(c(0.0, 0.0)),
as.integer(scanrange),
as.integer(length(scantime)),
as.double(ppm * 1e-6),
as.integer(minCentroids),
as.integer(prefilter),
as.integer(noise),
PACKAGE ='xcms' )
)
}
)
message("OK")
## ROI.list <- findmzROI(object,scanrange=scanrange,dev=ppm * 1e-6,minCentroids=minCentroids, prefilter=prefilter, noise=noise)
if (length(roiList) == 0) {
warning("No ROIs found! \n")
if (verboseColumns) {
nopeaks <- matrix(nrow = 0, ncol = length(basenames) +
length(verbosenames))
colnames(nopeaks) <- c(basenames, verbosenames)
} else {
nopeaks <- matrix(nrow = 0, ncol = length(basenames))
colnames(nopeaks) <- c(basenames)
}
return(invisible(nopeaks))
}
}
return(roiList)
}
.centWave_orig(mzs, ints, scant, valsPerSpect)
.centWave_orig <- function(mz, int, scantime, valsPerSpect,
ppm = 25, peakwidth = c(20,50), snthresh = 10,
prefilter = c(3,100), mzCenterFun = "wMean",
integrate = 1, mzdiff = -0.001, fitgauss = FALSE,
noise = 0, ## noise.local=TRUE,
sleep = 0, verboseColumns = FALSE, roiList = list(),
firstBaselineCheck = TRUE, roiScales = NULL,
extendLengthMSW = FALSE) {
## Input argument checking.
if (missing(mz) | missing(int) | missing(scantime) | missing(valsPerSpect))
stop("Arguments 'mz', 'int', 'scantime' and 'valsPerSpect'",
" are required!")
if (length(mz) != length(int) | length(valsPerSpect) != length(scantime)
| length(mz) != sum(valsPerSpect))
stop("Lengths of 'mz', 'int' and of 'scantime','valsPerSpect'",
" have to match. Also, 'length(mz)' should be equal to",
" 'sum(valsPerSpect)'.")
scanindex <- xcms:::valueCount2ScanIndex(valsPerSpect) ## Get index vector for C calls
if (!is.double(mz))
mz <- as.double(mz)
if (!is.double(int))
int <- as.double(int)
## Fix the mzCenterFun
# mzCenterFun <- paste("mzCenter",
#                      gsub(mzCenterFun, pattern = "mzCenter.",
#                           replacement = "", fixed = TRUE), sep=".")
# if (!exists(mzCenterFun, mode="function"))
#     stop("Function '", mzCenterFun, "' not defined !")
if (!is.logical(firstBaselineCheck))
stop("Parameter 'firstBaselineCheck' should be logical!")
if (length(firstBaselineCheck) != 1)
stop("Parameter 'firstBaselineCheck' should be a single logical !")
if (length(roiScales) > 0)
if (length(roiScales) != length(roiList) | !is.numeric(roiScales))
stop("If provided, parameter 'roiScales' has to be a numeric with",
" length equal to the length of 'roiList'!")
## if (!is.null(roiScales)) {
##     if (!is.numeric(roiScales) | length(roiScales) != length(roiList))
##         stop("Parameter 'roiScales' has to be a numeric of length equal to",
##              " parameter 'roiList'!")
##}
basenames <- c("mz", "mzmin", "mzmax", "rt", "rtmin", "rtmax",
"into", "intb", "maxo", "sn")
verbosenames <- c("egauss", "mu", "sigma", "h", "f", "dppm", "scale",
"scpos", "scmin", "scmax", "lmin", "lmax")
## Peak width: seconds to scales
scalerange <- round((peakwidth / mean(diff(scantime))) / 2)
if (length(z <- which(scalerange == 0)))
scalerange <- scalerange[-z]
if (length(scalerange) < 1) {
warning("No scales? Please check peak width!")
if (verboseColumns) {
nopeaks <- matrix(nrow = 0, ncol = length(basenames) +
length(verbosenames))
colnames(nopeaks) <- c(basenames, verbosenames)
} else {
nopeaks <- matrix(nrow = 0, ncol = length(basenames))
colnames(nopeaks) <- c(basenames)
}
return(invisible(nopeaks))
}
if (length(scalerange) > 1)
scales <- seq(from = scalerange[1], to = scalerange[2], by = 2)
else
scales <- scalerange
minPeakWidth <-  scales[1]
noiserange <- c(minPeakWidth * 3, max(scales) * 3)
maxGaussOverlap <- 0.5
minPtsAboveBaseLine <- max(4, minPeakWidth - 2)
minCentroids <- minPtsAboveBaseLine
scRangeTol <-  maxDescOutlier <- floor(minPeakWidth / 2)
scanrange <- c(1, length(scantime))
## If no ROIs are supplied then search for them.
if (length(roiList) == 0) {
message("Detecting mass traces at ", ppm, " ppm ... ", appendLF = FALSE)
## flush.console();
## We're including the findmzROI code in this function to reduce
## the need to copy objects etc.
## We could also sort the data by m/z anyway; wouldn't need that
## much time. Once we're using classes from MSnbase we can be
## sure that values are correctly sorted.
withRestarts(
tryCatch({
tmp <- capture.output(
roiList <- .Call("findmzROI",
mz, int, scanindex,
as.double(c(0.0, 0.0)),
as.integer(scanrange),
as.integer(length(scantime)),
as.double(ppm * 1e-6),
as.integer(minCentroids),
as.integer(prefilter),
as.integer(noise),
PACKAGE ='xcms' )
)
},
error = function(e){
if (grepl("m/z sort assumption violated !", e$message)) {
invokeRestart("fixSort")
} else {
simpleError(e)
}
}),
fixSort = function() {
## Force ordering of values within spectrum by mz:
##  o split values into a list -> mz per spectrum, intensity per
##    spectrum.
##  o define the ordering.
##  o re-order the mz and intensity and unlist again.
## Note: the Rle split is faster than the "conventional" factor split.
splitF <- Rle(1:length(valsPerSpect), valsPerSpect)
mzl <- as.list(S4Vectors::split(mz, f = splitF))
oidx <- lapply(mzl, order)
mz <<- unlist(mapply(mzl, oidx, FUN = function(y, z) {
return(y[z])
}, SIMPLIFY = FALSE, USE.NAMES = FALSE), use.names = FALSE)
int <<- unlist(mapply(as.list(split(int, f = splitF)), oidx,
FUN=function(y, z) {
return(y[z])
}, SIMPLIFY = FALSE, USE.NAMES = FALSE),
use.names = FALSE)
rm(mzl)
rm(splitF)
tmp <- capture.output(
roiList <<- .Call("findmzROI",
mz, int, scanindex,
as.double(c(0.0, 0.0)),
as.integer(scanrange),
as.integer(length(scantime)),
as.double(ppm * 1e-6),
as.integer(minCentroids),
as.integer(prefilter),
as.integer(noise),
PACKAGE ='xcms' )
)
}
)
message("OK")
## ROI.list <- findmzROI(object,scanrange=scanrange,dev=ppm * 1e-6,minCentroids=minCentroids, prefilter=prefilter, noise=noise)
if (length(roiList) == 0) {
warning("No ROIs found! \n")
if (verboseColumns) {
nopeaks <- matrix(nrow = 0, ncol = length(basenames) +
length(verbosenames))
colnames(nopeaks) <- c(basenames, verbosenames)
} else {
nopeaks <- matrix(nrow = 0, ncol = length(basenames))
colnames(nopeaks) <- c(basenames)
}
return(invisible(nopeaks))
}
}
return(roiList)
}
.centWave_orig(mzs, ints, scant, valsPerSpect)
roi <- .centWave_orig(mzs, ints, scant, valsPerSpect)
View(roi)
roi <- .centWave_orig(mzs, ints, scant, valsPerSpect,
ppm = 3, prefilter = c(3,100),
peakwidth = c(10,60),
mzdiff = -1e-2,
snthresh = 6,
noise = 0)
View(roi)
roi <- do.call(rbind, roi)
View(roi)
roi <- .centWave_orig(mzs, ints, scant, valsPerSpect,
ppm = 3, prefilter = c(3,100),
peakwidth = c(10,60),
mzdiff = -1e-2,
snthresh = 6,
noise = 0)
res <- xcmsOlap("./HermesResults/Ecoli/MS1_positive/UNL_Ecoli_MS1_pos_70-1000_01.mzML")
roi <- do.call(rbind, roi)
View(roi)
View(roi)
roi <- as.data.frame(roi)
roi %<>% mutate(., rtmin = scmin * mean(diff(scant))*2) %>%
mutate(., rtmax = scmax * mean(diff(scant))*2)
roi <-  mutate(roi, rtmin = scmin * mean(diff(scant))*2) %>%
mutate(., rtmax = scmax * mean(diff(scant))*2)
View(roi)
View(roi)
roi <- apply(roi, 2, as.numeric)
roi <- as.data.frame(roi)
View(roi)
roi <-  mutate(roi, rtmin = scmin * mean(diff(scant))*2) %>%
mutate(., rtmax = scmax * mean(diff(scant))*2)
View(roi)
RHermesGUI()
View(pks)
pks <- as.data.frame(pks)
filter(pks, between(V1, 323.063, 323.065))
filter(pks, between(V1, 323.063, 323.065)) %>% plot
filter(pks, between(V1, 323.063, 323.065)) %>% plot(.[,2])
filter(pks, between(V1, 323.063, 323.065)) %>% plot(y=.[,2])
filter(pks, between(V1, 323.063, 323.065)) %>% plot(y=.$V2)
filter(pks, between(V1, 323.063, 323.065)) %>% plot(y=.$V2)
filter(pks, between(V1, 323.063, 323.065)) %>% plot(y=.$V2)
filter(pks, between(V1, 323.063, 323.065))$V2
filter(pks, between(V1, 323.063, 323.065))$V2 %>% plot
RHermesGUI()
RHermesGUI()
roi <-  mutate(roi, rtmin = scmin * mean(diff(scant))) %>%
mutate(., rtmax = scmax * mean(diff(scant)))
View(roi)
roi <-  mutate(roi, rtmin = scmin * mean(diff(scant))) %>%
mutate(., rtmax = scmax * mean(diff(scant))) %>%
mutate(., length = rtmax - rtmin)
View(res)
View(roi)
roi <-  mutate(roi, rtmin = scmin * mean(diff(scant))) %>%
mutate(., rtmax = scmax * mean(diff(scant))) %>%
mutate(., length = rtmax - rtmin)
roi <-  mutate(roi, rtmin = scmin * mean(diff(scant))) %>%
mutate(., rtmax = scmax * mean(diff(scant))) %>%
mutate(., rtlength = rtmax - rtmin)
View(roi)
mypeaks <- roi
raw  <- RHermes:::import_and_filter(file)[[3]]
cnt <- 0
pkscans <- list()
for (j in 1:nrow(mypeaks)) {
j <- mypeaks[j,]
rt1 <- j$rtmin; rt2 <- j$rtmax
mz1 <- j$mzmin; mz2 <- j$mzmax
cur_scans <- filter(raw, between(mz, mz1, mz2) &
between(rt, rt1, rt2))
cnt <- cnt + nrow(cur_scans)
pkscans <- c(pkscans, list(cur_scans))
}
pkscans <- rbindlist(pkscans)
View(pkscans)
res <- list(pkscans = pkscans)
total <- nrow(PL(ecoliPos,2)@raw)
anotated <- nrow(PL(ecoliPos,2)@peaklist %>% distinct())
crude <- select(SOI(ecoliPos, 2)@PlotDF, c("rtiv", "rt")) %>% distinct() %>% nrow()
clean <-  select(SOI(ecoliPos, 5)@PlotDF,c("rtiv", "rt")) %>% distinct() %>% nrow()
#Intersection with peaklist, which contains all "annotated data points" but where SOI
#have not been detected yet.
intersection <- dplyr::intersect(res$pkscans,
select(PL(ecoliPos, 2)@peaklist,
c("mz", "rtiv", "rt")))
#Intersection with "crude" SOI list that doesn't have blank substraction, refining, etc.
crudeSOIintersection <- dplyr::intersect(select(res$pkscans, c("rtiv", "rt")),
select(SOI(ecoliPos, 2)@PlotDF,
c("rtiv", "rt")))
#Intersection with "definitive" SOI list that goes to IL
SOIintersection <- dplyr::intersect(select(res$pkscans, c("rtiv", "rt")),
select(SOI(ecoliPos, 5)@PlotDF,
c("rtiv", "rt")))
SL <- SOI(ecoliPos, 5)@SOIList
points <- lapply(1:nrow(SL), function(soi){
cur <- SL[soi,]
filter(SOI(ecoliPos, 5)@PlotDF, form == cur$formula) %>%
filter(., between(rt, cur$start, cur$end)) %>%
mutate(., FrC = cur$FrC) %>%
mutate(., MIRS = cur$MIRS) %>%
mutate(., SOI = soi)
})
points <- do.call(rbind, points)
ggplot(points, aes(x=FrC, y=MIRS)) + geom_density_2d_filled()
xcms_lab <- lapply(1:nrow(res$pkscans), function(scan){
cur <- res$pkscans[scan,]
matching <- filter(points, rt == cur$rt & rtiv == cur$rtiv)
return(matching)
})
xcms_lab <- do.call(rbind, xcms_lab)
xcms_lab$SOI %>% unique() %>% length
table(xcms_lab$FrC > 0.5, xcms_lab$MIRS > 0.5) %>% prop.table()*100
points$SOI %>% unique() %>% length
table(points$FrC > 0.5, points$MIRS > 0.5) %>% prop.table()*100
pdf("./testSOIXCMS_Common_3.pdf")
for(soi in which(unique(points$SOI) %in% unique(xcms_lab$SOI))){
cur <- points[points$SOI == soi, ]
print(ggplot(cur, aes(x=rt, y=rtiv)) + geom_point() + scale_y_log10() +
ggtitle(paste("SOI", soi, "anot=", cur$form[1], "MIRS=",
cur$MIRS[1], "FrC=", cur$FrC[1])) + theme_minimal())
}
dev.off()
pdf("./testSOIXCMS_UniqueHermes_3.pdf")
for(soi in which(!unique(points$SOI) %in% unique(xcms_lab$SOI))){
cur <- points[points$SOI == soi, ]
print(ggplot(cur, aes(x=rt, y=rtiv)) + geom_point() + scale_y_log10() +
ggtitle(paste("SOI", soi, "anot=", cur$form[1], "MIRS=",
cur$MIRS[1], "FrC=", cur$FrC[1])) + theme_minimal())
}
dev.off()
print(paste("A", round(nrow(res$pkscans)/total*100, 3),
"% of all acquired data points have an associated XCMS peak"))
#Intersection with peaklist, which contains all "annotated data points" but where SOI
#have not been detected yet.
intersection <- dplyr::intersect(res$pkscans,
select(PL(ecoliPos, 2)@peaklist,
c("mz", "rtiv", "rt")))
anot_percentage <- nrow(intersection) / nrow(res$pkscans) * 100
print(paste("Only", round(anot_percentage,2),
"% of XCMS peak data points have an",
"ionic formula annotation present in ECMDB+KEGG"))
#Intersection with "crude" SOI list that doesn't have blank substraction, refining, etc.
crudeSOIintersection <- dplyr::intersect(select(res$pkscans, c("rtiv", "rt")),
select(SOI(ecoliPos, 2)@PlotDF,
c("rtiv", "rt")))
crudeSOIpercentage <- nrow(crudeSOIintersection) / nrow(res$pkscans) * 100
print(paste("Only", round(crudeSOIpercentage,2),
"% of XCMS peak data points are",
"represented in the crude SOI list"))
#Intersection with "definitive" SOI list that goes to IL
SOIintersection <- dplyr::intersect(select(res$pkscans, c("rtiv", "rt")),
select(SOI(ecoliPos, 5)@PlotDF,
c("rtiv", "rt")))
SOI_percentage <- nrow(SOIintersection) / nrow(res$pkscans) * 100
print(paste("Only", round(SOI_percentage,2),
"% of XCMS peak data points are",
"represented in the final SOI list"))
xs <- c(0, clean/total, crude/total, anotated/total, 1)
ys <- c(nrow(SOIintersection)/clean,
(nrow(crudeSOIintersection) - nrow(SOIintersection))/(total * (xs[3]-xs[2])),
(nrow(intersection) - nrow(crudeSOIintersection))/(total * (xs[4]-xs[3])),
(nrow(res$pkscans) - nrow(intersection))/(total * (xs[5]-xs[4])))
p <- ggplot()
p <- p + geom_rect(aes(xmin=0,xmax=1,min=0,max=1), alpha = 0.2, fill = "grey90") +
geom_rect(aes(xmin=xs[1], xmax=xs[2], ymin = 0, ymax = 1), fill = "darkred", alpha = 0.5) +
geom_rect(aes(xmin=xs[2], xmax=xs[3], ymin = 0, ymax = 1), fill = "red", alpha = 0.5) +
geom_rect(aes(xmin=xs[3], xmax=xs[4], ymin = 0, ymax = 1), fill = "pink", alpha = 0.5)
p <- p + geom_rect(aes(xmin=xs[1], xmax=xs[2], ymin = 0, ymax = ys[1]), fill = "blue", alpha = 0.3) +
geom_rect(aes(xmin=xs[2], xmax=xs[3], ymin = 0, ymax = ys[2]), fill = "blue", alpha = 0.3) +
geom_rect(aes(xmin=xs[3], xmax=xs[4], ymin = 0, ymax = ys[3]), fill = "blue", alpha = 0.3) +
geom_rect(aes(xmin=xs[4], xmax=xs[5], ymin = 0, ymax = ys[4]), fill = "blue", alpha = 0.3)
p <- p +
theme_minimal() + theme(axis.text = element_blank(),axis.ticks = element_blank(), panel.grid = element_blank())
p
View(roi)
RHermesGUI()
library(ggplot2)
library(plotly)
library(dplyr)
library(RHermes)
rho_chaos <- function(data, nlevels = 20, fillGaps = TRUE){
data$rtiv <- (data$rtiv - min(data$rtiv)) /
(max(data$rtiv) - min(data$rtiv))
ns <- sapply(seq_len(nlevels), function(n){
t_n <- n/nlevels
data$above <- data$rtiv > t_n
nr <- nrow(data)
#Filling 1-gaps
if(fillGaps){
for(i in seq_len(nr-2)){
if(data$above[i] & !data$above[i+1] & data$above[i+2]){
data$above[i+1] <- TRUE
}
}
}
bad <- which(!data$above)
if(length(bad) == 0){return(1)}
else{
add <- 1
if(length(bad) == 1){return(1+add)}
return(length(which(diff(bad) > 1)) + add)
}
})
rho <- 1 - sum(ns) / (nrow(data) * nlevels)
return(rho)
}
rho_spatial <- function(PL, data, f, ppm, ad){
data(isotopes, package = "enviPat")
iso <- enviPat::isopattern(f, isotopes = isotopes, verbose = FALSE)[[1]]
iso <- iso[order(iso[,2], decreasing = TRUE), ]
num_iso <- min(nrow(iso) - 1, 4)
iso <- iso[seq(num_iso) + 1,]
iso[,1] <- iso[,1] + ad
scores <- apply(iso, 1, function(x){
isodata <- filter(PL, between(mz,
x[1] * (1 - ppm * 1e-6),
x[1] * (1 + ppm * 1e-6)))
RHermes:::cosineSim(data, isodata) * x[2]
})
rho <- sum(scores) / sum(iso[, 2])
}
calculateMSM <- function(PL, data, f, ppm, ad){
r_c <- rho_chaos(data, nlevels = 20, fillGaps = TRUE)
if(max(data$rtiv) < 1e5){
r_spatial <- 1
} else {
r_spatial <- rho_spatial(PL, data, f, ppm, ad)
}
# msm <- r_c * r_spatial * r_spectral
msm <- r_c * r_spatial
return(msm)
}
filterFDR <- function(struct, id, fdr=0.05, iter = 20){
SOI <- SOI(struct, id)
SOIList <- SOI@SOIList
DB <- DB(struct)
PL <- PL(struct, which(struct@metadata@filenames == SOI@filename)[1])@raw
ppm <- struct@metadata@ExpParam@ppm
decoyList <- SOIList
message("Calculating MSM on Sample")
SOIList$MSM <- sapply(1:nrow(SOIList), function(i){
cur <- SOIList[i, ]
addelta <- adlist(struct)$mass[adlist(struct)$adduct == cur[["ad"]][[1]][1]]
return(calculateMSM(PL, cur[["peaks"]][[1]], cur[["f"]][[1]][1], ppm, addelta))
})
message("Calculating MSM on Decoy")
data(isotopes,package = "enviPat", envir = environment())
rare_ads <- isotopes[c(3:295), ]
rare_ads <- do.call(rbind, lapply(unique(rare_ads$element), function(atom){
cur <- rare_ads[rare_ads$element == atom,]
return(cur[which.max(cur$abundance),])
}))
decoy_MSM <- lapply(1:nrow(decoyList), function(i){
if(i%%100 == 0){message(paste(i, "out of", nrow(decoyList)))}
lapply(1:iter, function(j){
adID  <- sample(x = seq_len(nrow(rare_ads)), size = 1)
cur <- decoyList[i, ]
cur$mass <- DB$EnviPatMass[which(DB$MolecularFormula == cur$f[[1]][1])[1]] +
rare_ads$mass[adID]
# pks <- filter(PL, between(rt, cur$start, cur$end))
pks <- filter(PL, between(mz, cur$mass*(1 - ppm*1e-6),
cur$mass*(1 + ppm*1e-6)))
if(nrow(pks) < 3){return(0)}
return(calculateMSM(PL, pks, cur$f[[1]][1], ppm, rare_ads$mass[adID]))
})
})
thr <- c()
for(i in 1:iter){
dec <- data.frame(msm = sapply(decoy_MSM, function(x){x[[i]]}) %>% sort,
x=sapply(decoy_MSM, function(x){x[[i]]}) %>% seq_along(), class = "Decoy")
samp <- data.frame(msm = SOIList$MSM %>% sort,
x=1:nrow(SOIList), class = "Sample")
print(ggplot(rbind(dec, samp)) + geom_step(aes(x=x, y=msm, color = class), size=2,
position = "identity", alpha = 0.5)+
theme_minimal() + ylab("MSM") + xlab("Annotation"))
joint <- rbind(dec,samp)
joint <- joint[order(joint$msm, decreasing = TRUE),]
estimated_fdr <- c()
for(i in 1:nrow(samp)){
estimated_fdr <- c(estimated_fdr,
length(which(joint$class[1:i] == "Decoy")) / i)
}
print(ggplot(data.frame(FDR = estimated_fdr,
Annotation = seq_along(estimated_fdr)),
aes(x = Annotation, y=FDR)) +
geom_line(size = 1) +
theme_minimal() + scale_y_continuous(limits = c(0,1))
)
thr <- samp$msm[max(which(estimated_fdr < fdr) %>% max, nrow(samp))]
}
thr <- median(thr)
struct@data@SOI[[id]]@SOIList <- SOIList[SOIList$MSM > thr, ]
return(struct)
}
struct <- readRDS("D:/HermesResults/HumanPlasma/Hermes_objects/PosResults.rds")
id <- 3
pdf("PlasmaPositive_FDR.pdf")
struct <- filterFDR(struct, id, fdr = 0.2, iter = 5)
dev.off()
test <- MSnbase::readMSData(list.files("D:/HermesResults/Ecoli/MS1_positive/",pattern = "mzML", full.names = T))
test <- MSnbase::readMSData(list.files("D:/HermesResults/Ecoli/MS1_positive/", pattern = "mzML", full.names = T)[1])
test <- MSnbase::readMSData(list.files("D:/HermesResults/Ecoli/MS1_positive/", pattern = "mzML", full.names = T)[1], mode = "onDisk")
View(test)
library(MSnbase)
MSnbase::filterFile(test)
?MSnbase::filterFile(test)
test <- MSnbase::readMSData(list.files("D:/HermesResults/Ecoli/MS1_positive/", pattern = "mzML", full.names = T)[1:3], mode = "onDisk")
test2 <- MSnbase::filterFile(test, 1)
View(test2)
MSnbase::extractSpectraData(test2)
sp <- MSnbase::extractSpectraData(test2)
View(sp)