utils::globalVariables(c("error", "name", "out", "intervalRef", "signal", "signal0",
"signal1", "Mz"))
### plot ----
#' Plot a ptrSet object
#'
#' plot a ptrSet object
#' @aliases plot.ptrSet plot,ptrSet-method
#' @param x a ptrSet object
#' @param y not use
#' @param typePlot could be : \code{calibError}, \code{resolution},
#' \code{peakShape}, or
#' a empty character if you want all.
#' @return plot
#' @rdname plot
#' @export
#' @examples
#' dirRaw <- system.file("extdata/exhaledAir", package = "ptairData")
#' exhaledPtrset <- createPtrSet(dir=dirRaw, setName="exhaledPtrset",
#' mzCalibRef = c(21.022, 60.0525), fracMaxTIC = 0.7, saveDir = NULL )
#' plot(exhaledPtrset )
#' plot(exhaledPtrset ,typePlot='calibError')
#' plot(exhaledPtrset ,typePlot='resolution')
#' plot(exhaledPtrset ,typePlot='peakShape')
methods::setMethod(f = "plot", signature = "ptrSet", function(x, y, typePlot = "") {
# if(!is.null(saveFile)) { ext<- tools::file_ext(saveFile) if(ext == 'pdf')
# pdf(saveFile) if(ext == 'png') png(saveFile) if(ext == 'jpeg') jpeg(saveFile) }
if (!typePlot %in% c("", "calibError", "resolution", "peakShape"))
warning("typePolt should be calibError,resolution or peakShape")
if (typePlot == "calibError") {
return(plotCalibError(x))
} else if (typePlot == "resolution") {
return(plotResolution(x)$plot)
} else if (typePlot == "peakShape") {
return(plotPeakShape(x))
} else {
reso <- try(plotResolution(x))
calib <- try(plotCalibError(x))
peakShape <- try(plotPeakShape(x))
reaction <- try(plotPtrReaction(x))
left <- ggpubr::ggarrange(calib, reso$plot, reso$table, nrow = 3,
heights = c(0.5,
0.4, 0.1))
right <- ggpubr::ggarrange(peakShape, reaction, nrow = 2)
return(ggpubr::ggarrange(left, right, ncol = 2, align = "hv"))
}
# if(!is.null(saveFile)) dev.off()
})
# plot resolution boxplot for a ptrSet
plotResolution <- function(set) {
mzCalibRef <- getParameters(set)$mzCalibRef
resolution <- getPTRInfo(set)$resolution
# complete missing masses
index <- which(vapply(resolution, length, 1) != length(mzCalibRef))
resolution[index] <- lapply(resolution[index], function(x) {
missing <- which(!mzCalibRef %in% names(x))
x_new <- c(x, rep(as.numeric(NA), length(missing)))
names(x_new) <- c(names(x), mzCalibRef[missing])
x_new <- x_new[order(as.numeric(names(x_new)))]
x_new
})
# concatenate
resolutiontemp <- lapply(resolution, function(x) {
mat <- cbind(resolution = x, Mz = as.numeric(names(x)))
rownames(mat) <- NULL
mat
})
resolMat <- do.call(rbind, resolutiontemp)
resolMat <- data.table::as.data.table(resolMat)
`:=` <- data.table::`:=`
resolMat[, `:=`(Mz, factor(Mz))]
resolMat <- cbind(resolMat, name = rep(names(resolution),
each = nlevels(resolMat$Mz)))
# identify outliers
resolMat <- resolMat[!is.na(resolMat$resolution), ]
resolMat <- resolMat[, list(resolution, name, out = ifelse(is_outlier(resolution),
resolution, as.numeric(NA))), by = Mz]
# boxplot with outlier labels
g <- ggplot2::ggplot(subset(resolMat, !is.na(resolution)), ggplot2::aes(y = resolution,
x = Mz)) + ggplot2::geom_boxplot() + ggplot2::geom_text(data = resolMat[!is.na(out),
], ggplot2::aes(x = Mz, y = resolution, label = name), vjust = -0.5, size = 3) +
ggplot2::ggtitle(label = expression("Resolution m/" ~ Delta ~ "(m)")) + ggplot2::ylab("m/delta(m)") +
ggplot2::theme_classic()
resolutionEstimated <- Reduce(c, getPTRInfo(set)$resolution)
info <- gridExtra::tableGrob(t(data.frame(resolution = c(min = floor(min(resolutionEstimated)/1000) *
1000, mean = round(mean(resolutionEstimated)/1000) * 1000,
max = ceiling(max(resolutionEstimated)/1000) *
1000))), theme = gridExtra::ttheme_minimal(base_size = 12))
return(list(plot = g, table = info))
}
#' @importFrom rlang .data
plotCalibError <- function(set) {
massCalib <- getCalibrationInfo(set)$mzCalibRef
errorList <- getCalibrationInfo(set)$errorCalibPpm
names(errorList) <- paste(seq(1, length(errorList)), names(errorList),
sep = "-")
# get list files
listFiles <- getParameters(set)$listFile
if (methods::is(listFiles, "expression"))
listFiles <- eval(listFiles)
# number of files
nfiles <- length(massCalib)
# calibration reference
mzCalibRef <- getParameters(set)$mzCalibRef
# count how many time masses where excluded
count <- table(unlist(massCalib))
exclud <- nfiles - count
# TODO: affichier cette info sur le graph
# complete missing masses
index <- which(vapply(errorList, length, 1) != length(mzCalibRef))
errorList[index] <- lapply(errorList[index], function(x) {
missing <- which(!mzCalibRef %in% names(x))
x_new <- c(x, rep(NA, length(missing)))
names(x_new) <- c(names(x), mzCalibRef[missing])
x_new <- x_new[order(as.numeric(names(x_new)))]
x_new
})
# concatenate
errorListemp <- lapply(errorList, function(x) {
mat <- cbind(error = x, Mz = as.numeric(names(x)))
rownames(mat) <- NULL
mat
})
calibErrorMat <- do.call(rbind, errorListemp)
calibErrorMat <- data.table::as.data.table(calibErrorMat)
`:=` <- data.table::`:=`
calibErrorMat[, `:=`(Mz, factor(Mz))]
calibErrorMat <- cbind(calibErrorMat, name = rep(names(errorList),
each = nlevels(calibErrorMat$Mz)))
# identify outliers
calibErrorMat <- calibErrorMat[!is.na(calibErrorMat$error), ]
calibErrorMat <- calibErrorMat[, list(error, name,
out = ifelse(is_outlier(error),
error, as.numeric(NA))), by = Mz]
# boxplot with labels
p <- ggplot2::ggplot(subset(calibErrorMat, !is.na(error)),
ggplot2::aes(y = .data$error, x = .data$Mz)) +
ggplot2::geom_boxplot() +
ggplot2::geom_text(data = calibErrorMat[!is.na(out),],
ggplot2::aes(x = .data$Mz,
y = .data$error,
label = .data$name), vjust = -0.5, size = 3.5) +
ggplot2::ggtitle("Calibration error") + ggplot2::ylab("ppm") +
ggplot2::theme_classic()
return(p)
}
is_outlier <- function(x) {
return(x < stats::quantile(x, 0.25, na.rm = TRUE) - 1.5 * stats::IQR(x, na.rm = TRUE) |
x > stats::quantile(x, 0.75, na.rm = TRUE) + 1.5 * stats::IQR(x, na.rm = TRUE))
}
#' plot the average peak shape of reference calibration masses for a ptrSet
#' @param set ptrSet object
#' @param showAverage boolean
#' @return ggplot object
#' @importFrom scales hue_pal
#' @keywords internal
plotPeakShape <- function(set, showAverage = FALSE) {
mzRef <- getParameters(set)$mzCalibRef
n.files <- length(getCalibrationInfo(set)$mzCalibRef)
npoints <- 50
interval.ref <- seq(-3, 3, length = npoints) # +- 3 * peak width(=1)
# loop over file
peak_ref <- array(NA, dim = c(length(mzRef), n.files, npoints))
for (n.file in seq_len(n.files)) {
massRef <- getCalibrationInfo(set)$mzCalibRef[[n.file]]
interval <- getCalibrationInfo(set)$signalCalibRef[[n.file]]
names(interval) <- massRef
n.mass <- length(interval)
delta <- vapply(interval, function(x) width(tof = x$mz, peak = x$signal)$delta,
FUN.VALUE = 1.1)
t_centre <- vapply(massRef, function(x) {
sp <- interval[[as.character(x)]]$signal
mz <- interval[[as.character(x)]]$mz
localMax <- LocalMaximaSG(sp = sp, minPeakHeight = max(sp) * 0.2)
interpol <- stats::spline(mz[(localMax - 4):(localMax + 4)], sp[(localMax -
4):(localMax + 4)], n = 100)
sg <- signal::sgolayfilt(interpol$y, n = 51, p = 3) #n/2
center <- interpol$x[which.max(sg)]
return(center)
}, FUN.VALUE = 1.1)
# normalization of intreval
interval.n <- list()
for (j in seq_len(n.mass)) {
interval.n[[j]] <- (interval[[j]]$mz - t_centre[j])/delta[j]
}
names(interval.n) <- massRef
peak <- matrix(0, ncol = n.mass, nrow = length(interval.ref))
for (i in seq_len(n.mass)) {
# interpolation
interpolation <- stats::spline(interval.n[[i]], interval[[i]]$signal,
xout = interval.ref)
# baseline correction
interpolation$y <- interpolation$y - snipBase(interpolation$y, widthI = 4)
# normalization
indexPeakRef <- which(massRef == massRef[i])
peak_ref[indexPeakRef, n.file, ] <- interpolation$y/stats::spline(interval.ref,
interpolation$y, xout = 0)$y
# cumsum(interpolation$y)/sum(interpolation$y)
}
}
# aggreate
# mean
peaksMAT <- apply(peak_ref, 1, function(x) {
if (all(is.na(x)))
return(rep(NA, 4 * npoints))
peak <- apply(x, 2, function(x) mean(x, na.rm = TRUE))
peak <- stats::spline(interval.ref, peak, n = 4 * npoints)
return(peak$y)
})
interval.ref2 <- seq(-3, 3, length = 4 * npoints)
peaks <- matrix(peaksMAT, ncol = 1)
peakData <- data.frame(signal = peaks, Mz = as.factor(rep(round(mzRef, 3),
each = 4 *npoints)),
intervalRef = rep(interval.ref2, length(mzRef)))
peakData <- peakData[!is.na(peakData$signal), ]
# plot
if (showAverage) {
p <- ggplot2::ggplot() +
ggplot2::geom_line(data = peakData, ggplot2::aes(x = .data$intervalRef,
y = .data$signal,
group = .data$Mz,
colour = .data$Mz),
size = 1) +
ggplot2::geom_line(data = data.frame(x = interval.ref2,
y = apply(peaksMAT,
which.max(dim(peaksMAT)),
mean)),
ggplot2::aes(x = .data$x,y = .data$y, colour = "average"),
size = 1.2) +
ggplot2::ggtitle("Average normalized peak shape of
calibration peaks") +
ggplot2::xlab("Mz interval normalized") +
ggplot2::ylab("Intenisty normalized") +
ggplot2::scale_color_manual(values = c((scales::hue_pal())(length(levels(peakData$Mz))),
"black"))
} else {
p <- ggplot2::ggplot() +
ggplot2::geom_line(data = peakData,
ggplot2::aes(x = .data$intervalRef,
y = .data$signal, group = .data$Mz,
colour = .data$Mz), size = 1) +
ggplot2::ggtitle("Average normalized peak shape of
calibration peaks") +
ggplot2::xlab("Mz interval normalized") +
ggplot2::ylab("Intenisty normalized") +
ggplot2::scale_color_manual(values = c((scales::hue_pal())(length(levels(peakData$Mz)))))
}
if (n.files > 1) {
# confidence interval
peaksUp <- apply(peak_ref, 1, function(x) {
if (all(is.na(x)))
return(rep(NA, 4 * npoints))
peak <- apply(x, 2, function(y) mean(y, na.rm = TRUE) +
stats::qnorm(0.975) * stats::sd(y, na.rm = TRUE)/
sqrt(length(y)))
peak <- stats::spline(interval.ref, peak, n = 4 * npoints)
return(peak$y)
})
peaksDown <- apply(peak_ref, 1, function(x) {
if (all(is.na(x)))
return(rep(NA, 4 * npoints))
peak <- apply(x, 2, function(y) mean(y, na.rm = TRUE) - stats::qnorm(0.975) *
stats::sd(y, na.rm = TRUE)/sqrt(length(y)))
peak <- stats::spline(interval.ref, peak, n = 4 * npoints)
return(peak$y)
})
peaksUp <- matrix(peaksUp, ncol = 1)
peaksDown <- matrix(peaksDown, ncol = 1)
peakData <- data.frame(signal0 = peaksDown, signal1 = peaksUp,
Mz = as.factor(rep(round(mzRef, 3), each = 4 * npoints)),
intervalRef = rep(interval.ref2, length(mzRef)))
peakData <- peakData[!is.na(peakData$signal0), ]
p <- p + ggplot2::geom_line(data = peakData, ggplot2::aes(x = .data$intervalRef,
y = .data$signal0, group = .data$Mz, colour = .data$Mz), size = 0.6,
linetype = "dashed") + ggplot2::geom_line(data = peakData, mapping = ggplot2::aes(x = .data$intervalRef,
y = .data$signal1, group = .data$Mz, colour = .data$Mz), size = 0.6,
linetype = "dashed")
}
return(p + ggplot2::theme_classic())
}
plotPeakShapeTof <- function(set) {
mzRef <- getParameters(set)$mzCalibRef
n.files <- length(getCalibrationInfo(set)$mzCalibRef)
npoints <- 50
interval.ref <- seq(-3, 3, length = npoints) # +- 3 * peak width(=1)
# loop over file
peak_ref <- array(NA, dim = c(length(mzRef), n.files, npoints))
for (n.file in seq_len(n.files)) {
massRef <- getCalibrationInfo(set)$mzCalibRef[[n.file]]
interval <- alignCalibrationPeak(getCalibrationInfo(set)$signalCalibRef[[n.file]], massRef,
length(getTimeInfo(set)$TIC[[n.files]]))
names(interval) <- massRef
n.mass <- length(interval)
coefs <- getCalibrationInfo(set)$coefCalib[[n.file]]
# delta <- vapply(interval, function(x) width(tof =
# sqrt(x$mz)*coefs['a',]+coefs['b',], peak = x$signal)$delta, FUN.VALUE = 1.1)
delta <- vapply(interval, function(x) width(tof = x$mz, peak = x$signal)$delta,
FUN.VALUE = 1.1)
t_centre <- vapply(massRef, function(x) {
sp <- interval[[as.character(x)]]$signal
mz <- interval[[as.character(x)]]$mz
localMax <- LocalMaximaSG(sp = sp, minPeakHeight = max(sp) * 0.2)
interpol <- stats::spline(mz[(localMax - 4):(localMax + 4)], sp[(localMax -
4):(localMax + 4)], n = 1000)
sg <- signal::sgolayfilt(interpol$y, n = 501, p = 3)
tofMax <- interpol$x[which.max(sg)]
return(tofMax)
}, FUN.VALUE = 1.1)
# t_centre<-sqrt(t_centre)*coefs['a',]+coefs['b',]
# normalization of intreval
interval.n <- list()
for (j in seq_len(n.mass)) {
# interval.n[[j]]<-(sqrt(interval[[j]]$mz)*coefs['a',]+coefs['b',]-
# t_centre[j])/delta[j]
interval.n[[j]] <- (interval[[j]]$mz - t_centre[j])/delta[j]
}
names(interval.n) <- massRef
peak <- matrix(0, ncol = n.mass, nrow = length(interval.ref))
for (i in seq_len(n.mass)) {
# interpolation
interpolation <- stats::spline(interval.n[[i]], interval[[i]]$signal,
xout = interval.ref)
# baseline correction
interpolation$y <- interpolation$y - snipBase(interpolation$y, widthI = 4)
# normalization
indexPeakRef <- which(massRef == massRef[i])
peak_ref[indexPeakRef, n.file, ] <- interpolation$y/stats::spline(interval.ref,
interpolation$y, xout = 0)$y
# cumsum(interpolation$y)/sum(interpolation$y)
}
}
# aggreate
# mean
peaks <- apply(peak_ref, 1, function(x) {
if (all(is.na(x)))
return(rep(NA, 4 * npoints))
peak <- apply(x, 2, function(x) mean(x, na.rm = TRUE))
peak <- stats::spline(interval.ref, peak, n = 4 * npoints)
return(peak$y)
})
interval.ref2 <- seq(-3, 3, length = 4 * npoints)
peaks <- matrix(peaks, ncol = 1)
peakData <- data.frame(signal = peaks, Mz = as.factor(rep(mzRef, each = 4 * npoints)),
intervalRef = rep(interval.ref2, length(mzRef)))
peakData <- peakData[!is.na(peakData$signal), ]
# plot
p <- ggplot2::ggplot(data = peakData) +
ggplot2::geom_line(ggplot2::aes(x = intervalRef, y = signal,
group = Mz, colour = Mz), size = 1) +
ggplot2::ggtitle("Average normalized peak shape") +
ggplot2::xlab("Mz interval normalized") +
ggplot2::ylab("Intenisty normalized")
if (n.files > 1) {
# confidence interval
peaksUp <- apply(peak_ref, 1, function(x) {
if (all(is.na(x)))
return(rep(NA, 4 * npoints))
peak <- apply(x, 2, function(y) mean(y, na.rm = TRUE) +
stats::qnorm(0.975) * stats::sd(y, na.rm = TRUE)/
sqrt(length(y)))
peak <- stats::spline(interval.ref, peak, n = 4 * npoints)
return(peak$y)
})
peaksDown <- apply(peak_ref, 1, function(x) {
if (all(is.na(x)))
return(rep(NA, 4 * npoints))
peak <- apply(x, 2, function(y) mean(y, na.rm = TRUE) -
stats::qnorm(0.975) * stats::sd(y, na.rm = TRUE)/
sqrt(length(y)))
peak <- stats::spline(interval.ref, peak, n = 4 * npoints)
return(peak$y)
})
peaksUp <- matrix(peaksUp, ncol = 1)
peaksDown <- matrix(peaksDown, ncol = 1)
peakData <- data.frame(signal0 = peaksDown, signal1 = peaksUp,
Mz = as.factor(rep(mzRef,each = 4 * npoints)),
intervalRef = rep(interval.ref2, length(mzRef)))
peakData <- peakData[!is.na(peakData$signal0), ]
p <- p + ggplot2::geom_line(data = peakData, ggplot2::aes(x = intervalRef,
y = signal0, group = Mz, colour = Mz), size = 0.6,
linetype = "dashed") +
ggplot2::geom_line(data = peakData, ggplot2::aes(x = intervalRef,
y = signal1,
group = Mz, colour = Mz), size = 0.6, linetype = "dashed")
}
return(p)
}
plotPtrReaction <- function(pSet) {
# U <- Reduce(c, lapply(getPTRInfo(pSet)$prtReaction, function(x) {
# y <- c(x$TwData[1, , ])
# mean(y[y != 0])
# }))
# PD <- Reduce(c, lapply(getPTRInfo(pSet)$prtReaction, function(x) {
# y <- c(x$TwData[2, , ])
# mean(y[y != 0])
# }))
# TD <- Reduce(c, lapply(getPTRInfo(pSet)$prtReaction, function(x) {
# y <- c(x$TwData[3, , ])
# mean(y[y != 0])
# }))
# EN <- Reduce(c, lapply(getPTRInfo(pSet)$prtReaction, function(x) {
# y <- c(x$TwData[4, , ])
# mean(y[y != 0])
# }))
primaryIon <- Reduce(c, lapply(getPTRInfo(pSet)$primaryIon, function(x) x$primaryIon))
date <- Reduce(c, getDate(pSet))
date <- try(vapply(date, function(x) chron::chron(dates. = strsplit(x, " ")[[1]][1],
times. = strsplit(x, " ")[[1]][2], format = c("d/m/y", "h:m:s")),
FUN.VALUE = chron::chron(1)),silent = TRUE)
if(all(is.na(date) )) {
date <- Reduce(c, getDate(pSet))
date <- try(vapply(date, function(x) chron::chron(dates. = strsplit(x, " ")[[1]][1],
times. = strsplit(x, " ")[[1]][2], format = c("m/d/y", "h:m:s")),
FUN.VALUE = chron::chron(1)),silent = TRUE)
}
if(!is.null(attr(date,"condition"))){
date <- Reduce(c, getDate(pSet))
date <- vapply(date, function(x) chron::chron(dates. = strsplit(x, " ")[[1]][1],
times. = substring(gsub("[[:lower:]]",":", paste(strsplit(x, " ")[[1]][2:4],collapse = "")),first = 1,last = 8), format = c("d/m/y", "h:m:s")),
FUN.VALUE = chron::chron(1))
}
# Udrift <- ggplot2::ggplot() +
# ggplot2::geom_point(mapping = ggplot2::aes(x = .data$date,y = .data$U),
# data = data.frame(date = as.Date(chron::as.dates(date)),
# U = U)) +
# ggplot2::ggtitle("Drift voltage") + ggplot2::ylab("V") +
# ggplot2::theme_classic() +
# ggplot2::theme(title = ggplot2::element_text(size = 9))
#
# Tdrift <- ggplot2::ggplot() +
# ggplot2::geom_point(mapping = ggplot2::aes(x = .data$date,y = .data$TD),
# data = data.frame(date = as.Date(chron::as.dates(date)), TD = TD)) +
# ggplot2::ggtitle("Drift temperature") +
# ggplot2::ylab("degree") +
# ggplot2::theme_classic() +
# ggplot2::theme(title = ggplot2::element_text(size = 9))
#
# Pdrift <- ggplot2::ggplot() + ggplot2::geom_point(mapping = ggplot2::aes(x = .data$date,
# y = .data$PD), data = data.frame(date = as.Date(chron::as.dates(date)), PD = PD)) +
# ggplot2::ggtitle("Drift pressure") +
# ggplot2::ylab("mbar") +
# ggplot2::theme_classic() +
# ggplot2::theme(title = ggplot2::element_text(size = 9))
#
primaryIonPlot <- ggplot2::ggplot() +
ggplot2::geom_point(mapping = ggplot2::aes(x = .data$date, y = .data$cps),
data = data.frame(date = as.Date(chron::as.dates(date)),
cps = primaryIon)) +
ggplot2::ggtitle("Primary ion isotope intensity") +
ggplot2::xlab("Date") + ggplot2::theme_classic() +
ggplot2::theme(title = ggplot2::element_text(size = 9))
# blank<-ggplot2::element_blank()
# reaction <- ggpubr::ggarrange(Udrift +
# ggplot2::theme(axis.ticks.x =blank ,
# axis.text.x = blank,
# axis.title.x = blank),
# Tdrift + ggplot2::theme(axis.ticks.x =blank,
# axis.text.x = blank,
# axis.title.x = blank),
# Pdrift + ggplot2::theme(axis.ticks.x = blank,
# axis.text.x = blank,
# axis.title.x = blank),
# primaryIonPlot, ncol = 1, heights = c(0.23, 0.23, 0.23, 0.31), align = "v")
return(primaryIonPlot)
}
### plotFiles----
#' Plot the calibration peaks after calibration
#'
#' @param object a ptrSet or ptrRaw object
#' @param ppm the width of plot windows
#' @param pdfFile is different of \code{NULL}, the file path to save the
#' plots in pdf
#' @param fileNames the name of the files in the ptrSet object to plot.
#' If \code{NULL}, all files will be
#' plotted
#' @param ... not used
#' @return plot
#' @export
#' @rdname plotCalib
#' @examples
#' ## ptrSet
#' dirRaw <- system.file("extdata/exhaledAir", package = "ptairData")
#' exhaledPtrset <- createPtrSet(dir=dirRaw, setName="exhaledPtrset",
#' mzCalibRef = c(21.022, 60.0525), fracMaxTIC = 0.7, saveDir = NULL )
#' plotCalib(exhaledPtrset ,fileNames=getFileNames(exhaledPtrset )[1])
#'
#' ## ptrRaw
#' filePath<-system.file('extdata/exhaledAir/ind1/ind1-1.h5',
#' package = 'ptairData')
#' raw <- readRaw(filePath,mzCalibRef=c(21.022,59.049))
#' plotCalib(raw)
methods::setMethod(f = "plotCalib", signature = "ptrSet", function(object, ppm = 2000,
pdfFile = NULL, fileNames = NULL, ...) {
fileNamesObject <- getParameters(object)$listFile
if (methods::is(fileNamesObject, "expression"))
fileNamesObject <- eval(fileNamesObject)
set <- object
# get list files
if (is.null(fileNames)) {
fileNames <- basename(fileNamesObject)
} else {
# check if fileNames are in the ptrSet object put in basename
fileNames <- basename(fileNames)
if (methods::is(fileNames, "expression"))
fileNames <- eval(fileNames)
test <- fileNames %in% basename(fileNamesObject)
if (!all(test))
stop("This file(s) names are not in
the directory: \n",
paste(!fileNames[test], collapse = "\n"))
}
# get information for plot
massCalib <- getCalibrationInfo(set)$mzCalibRef[fileNames]
# spectrCalib<-lapply(fileNames, function(x)
# alignCalibrationPeak(getCalibrationInfo(set)$signalCalibRef[[x]], length(getTimeInfo(set)$TIC[[x]])))
# names(spectrCalib)<-fileNames
spectrCalib <- getCalibrationInfo(set)$signalCalibRef[fileNames]
errorList <- getCalibrationInfo(set)$errorCalibPpm[fileNames]
# get the extension of the file
if (!is.null(pdfFile)) {
# check is pdf extension
file_ext <- tools::file_ext(pdfFile)
if (file_ext != "pdf") {
warning("the plot will be saved in pdf extension")
pdfFile <- gsub(paste0(".", file_ext, "$"), ".pdf", pdfFile)
}
pdf(file = pdfFile, width = 7 * 1.8, height = 5 * 1.8)
}
# plot all masses use for calibration for each files loop over files
for (i in seq_along(massCalib)) {
error <- errorList[[i]]
# plot in a window of 2000 ppm
nb_plot <- length(massCalib[[i]])
nb_col <- min(3, nb_plot)
nb_row <- ceiling(nb_plot/nb_col)
graphics::par(oma = c(0, 0, 3, 0))
layout(matrix(seq(1, nb_row * nb_col), nrow = nb_row, ncol = nb_col,
byrow = TRUE))
for (j in seq_along(massCalib[[i]])) {
m <- massCalib[[i]][j]
th <- m * (ppm/2)/10^6
mz <- spectrCalib[[i]][[j]]$mz
index <- which(m - th < mz & mz < m + th)
x <- mz[index]
y <- spectrCalib[[i]][[j]]$signal[index]
plot(x, y, type = "l", lwd = 2, ylab = "intenisty",
xlab = "mz", main = c(m, paste("error:", round(error[j], 2),
"ppm")))
abline(v = m, col = "red", lwd = 2)
} #end loop masses
title(main = paste(i, "-", names(spectrCalib)[i]), outer = TRUE,
line = 0.5,
cex.main = 2)
} #end loop files
if (!is.null(pdfFile))
dev.off()
graphics::par(oma = c(0, 0, 0, 0))
})
#' plot the Total Ion sptectrum (TIC) for one or several files.
#' @param object ptrSet or ptrRaw S4 object
#' @param type set 'plotly' to get an interactive plot, and 'ggplot' for
#' classical plot.
#' @param baselineRm logical. If \code{TRUE}, remove the baseline of the TIC
#' @param showLimits logical. If \code{TRUE}, add the time limits to the plot
#' (obtain with the `fracMaxTIC` argument or `createPtrSet` function)
#' @return a plotly of ggplot2 object.
#' @param pdfFile a absolute file path. A pdf will be generated with a plot
#' for each file, caints TIC and
#'time limits.
#' @param fileNames vector of character. The file names of the ptrSer that
#' you want to plot. Could be in
#' basename or fullname.
#' @param colorBy character. A name of the ptrSet's sampleMetaData column,
#' to display with
#' the same color files of same attributes.
#' @param normalizePrimariIon should the TIC be normalized by the primary ion
#' @param ... not used
#' @rdname plotTIC
#' @importFrom grDevices dev.off pdf
#' @export
#' @examples
#' dirRaw <- system.file("extdata/exhaledAir", package = "ptairData")
#' exhaledPtrset <- createPtrSet(dir=dirRaw, setName="exhaledPtrset",
#' mzCalibRef = c(21.022, 60.0525), fracMaxTIC = 0.7, saveDir = NULL )
#' plotTIC(exhaledPtrset ,type='ggplot')
methods::setMethod(f = "plotTIC", signature = "ptrSet",
function(object, type, baselineRm,
showLimits, pdfFile = NULL,
fileNames = NULL, colorBy = "rownames",
normalizePrimariIon = FALSE, ...) {
fileNamesObject <- getParameters(object)$listFile
if (methods::is(fileNamesObject, "expression"))
fileNamesObject <- eval(fileNamesObject)
set <- object
# get list files
if (is.null(fileNames)) {
fileNames <- fileNamesObject
} else {
# check if fileNames are in the ptrset object
test <- basename(fileNames) %in% basename(fileNamesObject)
if (!all(test))
stop("This file(s) are not in the directory: \n",
paste(basename(fileNames)[which(!test)],
collapse = "\n"))
# put in full name
if (any(dirname(fileNames) == ".")) {
fileNames <- fileNamesObject[basename(fileNamesObject) %in%
basename(fileNames)]
}
}
# get the TIC and time limit
TIC <- getTimeInfo(set)$TIC[basename(fileNames)]
indLim <- getTimeInfo(set)$timeLimit[basename(fileNames)]$exp
# getSampleMetadata
SMD <- getSampleMetadata(set)
# test if colorBy is a colnames of samplemetadata
if (colorBy != "rownames" & !(colorBy %in% colnames(SMD))) {
message(colorBy, " is not a column of sample metadata")
colorBy <- "rownames"
}
# remove baseline
if (baselineRm) {
TIC <- lapply(TIC, function(x) {
bl <- try(baselineEstimation(x, d = 1))
if (is.null(attr(bl, "condition")))
TIC.blrm <- x - bl else TIC.blrm <- x - x[1]
TIC.blrm
})
}
if (!is.null(pdfFile)) {
# check is pdf extension
file_ext <- tools::file_ext(pdfFile)
if (file_ext != "pdf") {
warning("the plot will be saved in pdf extension")
pdfFile <- gsub(paste0(".", file_ext, "$"), ".pdf", pdfFile)
}
pdf(pdfFile, width = 7 * 1.8, height = 5 * 1.8)
}
p <- ggplot2::ggplot()
for (j in seq_along(TIC)) {
ticPlot <- TIC[[j]]
time <- as.numeric(names(ticPlot))
if (normalizePrimariIon)
ticPlot <- ticPlot/getPTRInfo(set)$primaryIon[[j]]
if (!is.null(pdfFile)) {
plot <- ggplot2::qplot(x = time, y = ticPlot/(time[2] - time[1]),
xlab = "time", ylab = "Sum cps",
main = paste(j, names(TIC)[j], sep = " - "),
size = 0.8)
if (ncol(indLim[[j]]) > 0)
plot <- plot +
ggplot2::geom_vline(ggplot2::aes(xintercept = as.numeric(names(TIC[[j]]))[c(indLim[[j]])],
color = "time limits")) +
ggplot2::scale_fill_manual("Legend") +
ggplot2::theme(plot.title = ggplot2::element_text(size = 20,
face = "bold"),
axis.title = ggplot2::element_text(size = 16),
axis.text = ggplot2::element_text(size = 14),
legend.text = ggplot2::element_text(size = 14),
legend.title = ggplot2::element_text(size = 16))
print(plot)
}
if (colorBy == "rownames") {
colour <- names(TIC)[j]
} else {
colour <- as.factor(SMD[names(TIC)[j], colorBy])
}
data <- data.frame(time = time,
Sum_cps = ticPlot/(time[2] - time[1]),
Legend = colour)
p <- p + ggplot2::geom_point(mapping = ggplot2::aes(x = .data$time,
y = .data$Sum_cps,
color = .data$Legend),
data = data) +
ggplot2::geom_line(mapping = ggplot2::aes(x = .data$time,
y = .data$Sum_cps,
color = .data$Legend),
data = data, size = 1)
}
if (!is.null(pdfFile))
dev.off()
if (showLimits) {
if (colorBy == "rownames") {
Legend <- names(TIC)
} else {
Legend <- as.factor(SMD[, colorBy])
}
limitdf <- lapply(Legend,
function(j) {
t<-as.numeric(names(getTimeInfo(set)$TIC[[j]]))[c(getTimeInfo(set)$timeLimit[[j]]$exp)]
data.frame(x = t, Legend = rep(j, 2))}
)
limitdf <- Reduce(rbind, limitdf)
p <- p +
ggplot2::geom_vline(mapping = ggplot2::aes(xintercept = .data$x,
color = .data$Legend), data = limitdf, size = 0.9)
}
# set title and legend
p <- p + ggplot2::ggtitle(paste("TIC of", getParameters(set)$name)) +
ggplot2::theme_classic()
switch(type, ggplot = return(p), plotly = return(plotly::ggplotly(p)))
} #end function
)
#' Plot method for 'ptrRaw' objects
#'
#' Displays the image of the matrix of intensities, the TIC and the TIS,
#' for the selected m/z and time ranges
#'
#' @param object An S4 object of class \code{ptrRaw-class} or \code{ptrSet}
#' @param mzRange Either a vector of 2 numerics indicating the m/z limits
#' or an integer giving a nominal m/z
#' @param timeRange Vector of 2 numerics giving the time limits
#' @param type Character: plot type; either 'classical' [default] or 'plotly'
#' @param ppm Integer: Half size of the m/z window when mzRange is set to a
#' nominal mass
#' @param palette Character: Color palette for the 'classical' plot; either 'heat'
#' [default], 'revHeat', 'grey', 'revGrey' or 'ramp'
#' @param showVocDB Logical: Should putative m/z annotations from the internal
#' package database be displayed (default is TRUE)
#' @param figure.pdf Character: Either 'interactive' [default], or the filename
#' of the figure to be saved (with the 'pdf' extension); only available for the
#' 'classical' display
#' @return Invisibly returns a list of the raw (sub)matrix 'rawsubM' and
#' the voc (sub)database 'vocsubDB'
#' @param fileNames vector of character. The file names of the ptrSer that you
#' want to plot. Could be in
#' basename or fullname.
#' @param ... not used
#' @rdname plotRaw
#' @examples
#' dirRaw <- system.file("extdata/exhaledAir", package = "ptairData")
#' exhaledPtrset <- createPtrSet(dir=dirRaw, setName="exhaledPtrset",
#' mzCalibRef = c(21.022, 60.0525), fracMaxTIC = 0.7, saveDir = NULL )
#' ptairMS::plotRaw(exhaledPtrset ,mzRange=59,fileNames='ind1-1.h5')
#'
#' patientRaw <- ptairMS::readRaw(system.file('extdata/exhaledAir/ind1/ind1-1.h5',
#' package = 'ptairData'), mzCalibRef=c(21.022,59.049,75.05))
#' ptairMS::plotRaw(patientRaw, mzRange = 59)
#' ptairMS::plotRaw(patientRaw, mzRange = 59, type = 'plotly')
#'
#' @importFrom grDevices dev.off pdf
#' @importFrom graphics abline layout par title
#' @export
methods::setMethod(f = "plotRaw", signature = "ptrSet",
function(object, mzRange, timeRange = c(NA, NA),
type = c("classical", "plotly")[1], ppm = 2000,
palette = c("heat", "revHeat", "grey", "revGrey",
"ramp")[1], showVocDB = TRUE,
figure.pdf = "interactive",
fileNames = NULL, ...) {
set <- object
fileNamesObject <- getParameters(set)$listFile
if (methods::is(fileNamesObject, "expression"))
fileNamesObject <- eval(fileNamesObject)
# get list files
if (is.null(fileNames)) {
fileNames <- fileNamesObject
} else {
# check if fileNames are in the ptrset object
test <- basename(fileNames) %in% basename(fileNamesObject)
if (!all(test))
stop("This file(s) are not in the directory: \n",
paste(basename(fileNames)[which(!test)],
collapse = "\n"))
# put in full name
if (any(dirname(fileNames) == ".")) {
fileNames <- fileNamesObject[basename(fileNamesObject) %in%
basename(fileNames)]
}
}
if (length(mzRange) == 1)
mzRange <- mzRange + c(-1, 1) * mzRange * ppm * 1e-06
if (figure.pdf != "interactive") {
if (type == "plotly")
stop("'plotly display is only available in the 'interactive'
mode currently.",
call. = FALSE)
filenameSplitVc <- unlist(strsplit(basename(figure.pdf), ".", fixed = TRUE))
extC <- utils::tail(filenameSplitVc, 1)
if (extC == "pdf") {
pdf(figure.pdf)
} else stop("The extension of the 'figure.pdf' filename argument
should be 'pdf'",
call. = FALSE)
} else if (is.null(fileNames))
graphics::par(ask = TRUE)
for (file in fileNames) {
# open a mz range
name<-rhdf5::h5ls(file)
if(name$group[2] == "/AcquisitionLog"){
mz <- rhdf5::h5read(file, name = "FullSpectra/MassAxis")
FirstcalibCoef <- try(rhdf5::h5read(file, "FullSpectra/MassCalibration", index = list(NULL,
1)))
attributCalib <- try(rhdf5::h5readAttributes(file, "/FullSpectra"))
if (!is.null(attr(FirstcalibCoef, "condition")) & is.null(attr(attributCalib, "condition"))) {
FirstcalibCoef <- matrix(c(attributCalib$`MassCalibration a`, attributCalib$`MassCalibration b`),
ncol = 1)
}
} else if (name$group[2] == "/AddTraces"){
CalibInfo <-rhdf5::h5read(file, "/CALdata", index = list(NULL,1))
if(dim(CalibInfo$Spectrum)[1]!=0) FirstcalibCoef<-as.matrix(c(CalibInfo$Spectrum[,1],2),ncol=1) else{
mzRef <- CalibInfo$Mapping[1,]
tofRef <- CalibInfo$Mapping[2,]
a <- (tofRef[2] - tofRef[1]) / (sqrt(mzRef[2])-sqrt(mzRef[1]))
b<- tofRef[1] - sqrt(mzRef[1])*a
FirstcalibCoef<-as.matrix(c(a,b,2),ncol=1)}
rownames(FirstcalibCoef)<-c("a","b","q")
mz <- tofToMz(seq(0,
as.numeric(name[name$name=="AverageSpec" ,"dim"])-1),
calibCoef = FirstcalibCoef)
}
mzRange[1] <- max(min(mz), mzRange[1], na.rm = TRUE)
mzRange[2] <- min(max(mz), mzRange[2], na.rm = TRUE)
# get the index of the aroud mzRange
index <- which(mzRange[1] < mz & mz < mzRange[2])
# object data
if(name$group[2] == "/AcquisitionLog"){object <- rhdf5::h5read(file, name = "/FullSpectra/TofData",
index = list(index, NULL, NULL, NULL))
} else if(name$group[2] == "/AddTraces") {
object <- rhdf5::h5read(file, name = "/SPECdata/Intensities",
index = list(index, NULL))
}
# time axis
timeVn <- as.numeric(names(getTimeInfo(set)$TIC[[basename(file)]]))
timeiNTER <- (timeVn[3] - timeVn[2])
timeRange.file <- rep(0, 0)
timeRange.file[1] <- max(min(timeVn), timeRange[1], na.rm = TRUE)
timeRange.file[2] <- min(max(timeVn), timeRange[2], na.rm = TRUE)
indexTime <- which(timeVn >= timeRange.file[1] &
timeVn <= timeRange.file[2])
timeVn <- timeVn[indexTime]
# calibrate mass axis
tof <- sqrt(mz) * FirstcalibCoef[1, 1] + FirstcalibCoef[2, 1]
# tof<- seq(0,length(mz))
coefCalib <- getCalibrationInfo(set)$coefCalib[[basename(file)]][[1]]
mzVn <- ((tof - coefCalib["b", ])/coefCalib["a", ])^coefCalib["q", ]
mzVn <- mzVn[index]
# formate object matix
rawSubMN <- matrix(object, nrow = dim(object)[1],
ncol = prod(utils::tail(dim(object), 2)))
rawSubMN <- rawSubMN[, indexTime]
dimnames(rawSubMN) <- list(mzVn, timeVn)
if (showVocDB) {
vocdbDF <- .loadVocDB()
vocdb_sel.vl <- vocdbDF[, "ion_mass"] >= mzRange[1] &
vocdbDF[, "ion_mass"] <= mzRange[2]
if (sum(vocdb_sel.vl)) {
vocdbDF <- vocdbDF[vocdb_sel.vl, , drop = FALSE]
} else vocdbDF <- NULL
} else vocdbDF <- NULL
switch(type, classical = {
imageMN <- t(rawSubMN)[, seq_len(nrow(rawSubMN)), drop = FALSE]
rownames(imageMN) <- round(as.numeric(rownames(imageMN)))
colnames(imageMN) <- round(as.numeric(colnames(imageMN)), 4)
paletteVc <- .palette(palette = palette)
currentParLs <- graphics::par()
for (parC in c("cin", "cra", "csi", "cxy", "din", "page")) currentParLs[[parC]] <- NULL
marLs <- list(chr = c(0.6, 4.1, 4, 0), sca = c(0.6, 0.6, 4, 7.1), ima = c(4.1,
4.1, 0, 0), spe = c(4.1, 0.6, 0, 1.1))
graphics::par(font = 2, font.axis = 2, font.lab = 2, pch = 18)
layout(matrix(c(1, 2, 3, 4), byrow = TRUE, nrow = 2), heights = c(2,
5), widths = c(5, 2))
## chr: Chromatogram
graphics::par(mar = marLs[["chr"]])
chromVn <- apply(rawSubMN, 2, function(intVn) mean(intVn, na.rm = TRUE))
plot(as.numeric(names(chromVn)), chromVn, cex = 0.7, pch = 16, xlab = "",
ylab = "", xaxs = "i", xaxt = "n", yaxs = "i")
graphics::mtext("Mean of intensity", cex = 0.8, side = 2, line = 2.5)
## sca: Color scale
graphics::par(mar = marLs[["sca"]])
.drawScale(imageMN = imageMN, paletteVc = paletteVc)
## ima: Image
graphics::par(mar = marLs[["ima"]])
.drawImage(imageMN = imageMN, paletteVc = paletteVc)
if (showVocDB && !is.null(vocdbDF)) {
mzImaVn <- as.numeric(colnames(imageMN))
abline(h = vapply(vocdbDF[, "ion_mass"],
function(mzN) (mzN - min(mzImaVn))/
diff(range(mzImaVn)) * ncol(imageMN) +
par("usr")[1],
FUN.VALUE = 1.1), lty = "dotted")
}
## spe: Spectrum
graphics::par(mar = marLs[["spe"]])
specVn <- apply(rawSubMN, 1, function(intVn) mean(intVn, na.rm = TRUE)/timeiNTER)
plot(specVn, as.numeric(names(specVn)), cex = 0.7, pch = 16, xlab = "",
ylab = "", xaxs = "i", yaxs = "i", yaxt = "n")
graphics::mtext("Count Per Second (cps)", cex = 0.8, side = 1, line = 2.5)
if (showVocDB && !is.null(vocdbDF)) {
abline(h = vocdbDF[, "ion_mass"], lty = "dotted")
}
title(main = basename(file), outer = TRUE, line = -2, cex.main = 2)
graphics::par(currentParLs)
}, plotly = {
plotlyBuild <- function() {
p <- plotly::subplot(plotly::plot_ly(x = as.numeric(colnames(rawSubMN)),
y = colSums(rawSubMN),
type = "scatter",
mode = "lines"),
plotly::plotly_empty(),
plotly::plot_ly(z = rawSubMN,
x = as.numeric(colnames(rawSubMN)),
y = as.numeric(rownames(rawSubMN)),
type = "heatmap"),
plotly::plot_ly(x = rowSums(rawSubMN),
y = as.numeric(rownames(rawSubMN)),
type = "scatter",
mode = "lines"),
nrows = 2, heights = c(0.2, 0.8),
widths = c(0.8, 0.2), shareX = TRUE,
shareY = TRUE, titleX = FALSE,
titleY = FALSE)
`%>%` <- plotly::`%>%`
p <- p %>% plotly::layout(title = basename(file))
return(p)
}
p <- suppressWarnings(suppressMessages(plotlyBuild()))
p <- plotly::layout(p, showlegend = FALSE)
print(plotly::ggplotly(p))
})
} #end loop file
# reset parameter plot
graphics::par(ask = FALSE)
# close pdf
if (figure.pdf != "interactive")
dev.off()
# print voc ref
if (!is.null(vocdbDF)) {
vocdbDF <- vocdbDF[nrow(vocdbDF):1, , drop = FALSE]
print(vocdbDF[, c("ion_mass", "ion_formula", "name_iupac"), drop = FALSE])
}
})
#' @rdname plotFeatures
#' @export
methods::setMethod(f = "plotFeatures", signature = "ptrSet",
function(set, mz,
typePlot,
addFeatureLine, ppm, pdfFile,
fileNames, colorBy) {
fileNamesObject <- getParameters(set)$listFile
if (methods::is(fileNamesObject, "expression"))
fileNamesObject <- eval(fileNamesObject)
# get list files
if (is.null(fileNames)) {
fileNames <- fileNamesObject
} else {
# check if fileNames are in the ptrSet object put in basename
test <- basename(fileNames) %in% basename(fileNamesObject)
if (!all(test))
stop("This file(s) names are not in
the directory: \n",
paste(!basename(fileNames)[test], collapse = "\n"))
# Put in fullNames
if (any(dirname(fileNames) == ".")) {
fileNames <- fileNamesObject[which(basename(fileNamesObject) %in% fileNames)]
}
}
# getSampleMetadata
SMD <- getSampleMetadata(set)
# test if colorBy is a colnames of samplemetadata
if (colorBy != "rownames" & !(colorBy %in% colnames(SMD))) {
message(colorBy, " is not a column of sample metadata")
colorBy <- "rownames"
}
# get the extension file
if (!is.null(pdfFile)) {
# check is pdf extension
file_ext <- tools::file_ext(pdfFile)
if (file_ext != "pdf") {
warning("the plot will be saved in pdf extension")
pdfFile <- gsub(paste0(".", file_ext, "$"), ".pdf", pdfFile)
}
pdf(file = pdfFile, width = 7 * 1.8, height = 5 * 1.8)
}
# loop over file
plotAll <- ggplot2::ggplot()
listPlotFile <- list()
j <- 0
for (file in fileNames) {
j <- j + 1
name<-rhdf5::h5ls(file)
if(name$group[2] == "/AcquisitionLog"){
mzAxis <- rhdf5::h5read(file, name = "FullSpectra/MassAxis")
thLarge <- max(0.4, mz * (ppm/2)/10^6)
indexMz <- which(mz - thLarge < mzAxis & mzAxis < mz + thLarge)
raw <- rhdf5::h5read(file, name = "/FullSpectra/TofData",
index = list(indexMz, NULL, NULL, NULL))
time <- c(rhdf5::h5read(file, name = "/TimingData/BufTimes"))
index_zero <- which(time == 0)[-1]
if (length(index_zero))
time <- time[-index_zero]
rawMn <- matrix(raw, nrow = dim(raw)[1], ncol = prod(utils::tail(dim(raw),
2)))
FirstcalibCoef <- try(rhdf5::h5read(file, "FullSpectra/MassCalibration", index = list(NULL, 1)))
attributCalib <- try(rhdf5::h5readAttributes(file, "/FullSpectra"))
if (!is.null(attr(FirstcalibCoef, "condition")) & is.null(attr(attributCalib, "condition"))) {
FirstcalibCoef <- matrix(c(attributCalib$`MassCalibration a`, attributCalib$`MassCalibration b`),
ncol = 1)
}
} else if (name$group[2] == "/AddTraces"){
CalibInfo <-rhdf5::h5read(file, "/CALdata", index = list(NULL,1))
if(dim(CalibInfo$Spectrum)[1]!=0) FirstcalibCoef<-as.matrix(c(CalibInfo$Spectrum[,1],2),ncol=1) else{
mzRef <- CalibInfo$Mapping[1,]
tofRef <- CalibInfo$Mapping[2,]
a <- (tofRef[2] - tofRef[1]) / (sqrt(mzRef[2])-sqrt(mzRef[1]))
b<- tofRef[1] - sqrt(mzRef[1])*a
FirstcalibCoef<-as.matrix(c(a,b,2),ncol=1)}
rownames(FirstcalibCoef)<-c("a","b","q")
mzAxis <- tofToMz(seq(0,
as.numeric(name[name$name=="AverageSpec" ,"dim"])-1),
calibCoef = FirstcalibCoef)
thLarge <- max(0.4, mz * (ppm/2)/10^6)
indexMz <- which(mz - thLarge < mzAxis & mzAxis < mz + thLarge)
time <- rhdf5::h5read(file, "/SPECdata/Times", bit64conversion = "bit64",
index = list(NULL, NULL))[4,]
rawMn <- rhdf5::h5read(file, "/SPECdata/Intensities", bit64conversion = "bit64",
index = list(indexMz, NULL))
}
# get the mass axis, must be the same for all files
# get the index of the mzAxis
indLim <- getTimeInfo(set)$timeLimit[[basename(file)]]$exp
indLimBg <- getTimeInfo(set)$timeLimit[[basename(file)]]$backGround
n.limit <- dim(indLim)[2]
mzAxis.j <- mzAxis[indexMz]
# * 0.2 ns / 2.9 (single ion signal) if convert to cps
tof <- sqrt(mzAxis.j) * FirstcalibCoef[1, 1] + FirstcalibCoef[2, 1]
coefCalib <- getCalibrationInfo(set)$coefCalib[[basename(file)]][[1]]
mzNew <- ((tof - coefCalib["b", ])/coefCalib["a", ])^coefCalib["q", ]
# get smaller windows
th <- mz * (ppm/2)/10^6
indexSub <- which(mz - th < mzNew & mzNew < mz + th)
# plot mean spectrum for expriration/time periods
plotFile <- ggplot2::ggplot()
indexTimeVec <- NULL
for (i in seq_len(n.limit)) {
# get index of the time period
indexTime <- seq(indLim["start", i], indLim["end", i])
indexTimeVec <- c(indexTimeVec, indexTime)
spectrum <- rowSums(rawMn[, indexTime])/(
ncol(rawMn[, indexTime]) * (time[3] - time[2]))
data <- data.frame(mz = mzNew[indexSub], cps = spectrum[indexSub],
timePeriods = as.character(i))
plotFile <- plotFile + ggplot2::geom_point(mapping = ggplot2::aes(x = .data$mz,
y = .data$cps,
color = .data$timePeriods),
data = data) +
ggplot2::geom_line(mapping = ggplot2::aes(x = .data$mz,
y = .data$cps,
color = .data$timePeriods),
data = data, size = 1)
}
# get the background
if (!is.null(indLimBg)) {
background <- rowSums(rawMn[, indLimBg])/
(length(indLimBg) * (time[3] - time[2]))
data = data.frame(mz = mzNew[indexSub], cps = background[indexSub],
timePeriods = "background")
# plot background to the file plot
plotFile <- plotFile +
ggplot2::geom_point(mapping = ggplot2::aes(x = .data$mz,
y = .data$cps,
color = .data$timePeriods),
data = data) +
ggplot2::geom_line(mapping = ggplot2::aes(x = .data$mz,
y = .data$cps,
color = .data$timePeriods),
data = data, linetype = "dashed")
}
plotFile <- plotFile + ggplot2::ggtitle(basename(file))
listPlotFile[[j]] <- plotFile
## add summary line to plotAll get the average time periods spectrum
spectrum <- rowSums(rawMn[, indexTimeVec])/(length(indexTimeVec) * (time[3] -
time[2]))
# spline intrepolation for spectrum
splineInterpol <- stats::splinefun(mzNew, spectrum)
# plot spectrum
if (colorBy == "rownames") {
colour <- basename(file)
} else {
colour <- SMD[basename(file), colorBy]
}
data = data.frame(mz = mzNew[indexSub], cps = spectrum[indexSub],
Legend = rep(colour,length(indexSub)))
plotAll <- plotAll +
ggplot2::geom_point(ggplot2::aes(x = .data$mz,
y = .data$cps,
color = .data$Legend),
data = data) +
ggplot2::stat_function(mapping = ggplot2::aes(x = .data$mz,
color = .data$Legend),
data = data, fun = splineInterpol,
n = 1000, size = 1)
# spline interpolation for background
if (!is.null(indLimBg)) {
background <- rowSums(rawMn[, indLimBg])/
(length(indLimBg) * (time[3] - time[2]))
splineInterpol <- stats::splinefun(mzNew, background)
# plot background
data = data.frame(mz = mzNew[indexSub],
Legend = rep(colour, length(indexSub)))
plotAll <- plotAll +
ggplot2::stat_function(mapping = ggplot2::aes(x = .data$mz,
color = .data$Legend),
data = data,
fun = splineInterpol, n = 1000,
linetype = "dashed", size = 1)
}
# plot vertical line line
if (addFeatureLine)
plotAll <- plotAll + ggplot2::geom_vline(xintercept = mz)
} #END loop file
plotAll <- plotAll + ggplot2::ggtitle(paste("Features", round(mz, 3),
"of", getParameters(set)$name),
subtitle = paste(annotateVOC(mz)[2:3], collapse = " ")) +
ggplot2::labs(color = colorBy) +
ggplot2::theme_classic()
if (!is.null(pdfFile)) {
for (j in seq_along(listPlotFile)) print(listPlotFile[[j]])
dev.off()
}
switch(typePlot, ggplot = return(plotAll),
plotly = return(plotly::ggplotly(plotAll)))
})
## sampleMetadata ----
#' reset the default sampleMetadata, containing the suborders names and the
#' acquisition dates as columns.
#' @param ptrset a ptrser object
#' @return a data.frame
#' @export
#' @examples
#' dirRaw <- system.file("extdata/exhaledAir", package = "ptairData")
#' exhaledPtrset <- createPtrSet(dir=dirRaw, setName="exhaledPtrset",
#' mzCalibRef = c(21.022, 60.0525), fracMaxTIC = 0.7, saveDir = NULL )
#' SMD<- resetSampleMetadata(exhaledPtrset)
resetSampleMetadata <- function(ptrset) {
dir <- getParameters(ptrset)$dir
fileNamesObject <- getParameters(ptrset)$listFile
if (methods::is(dir, "expression")) {
dir <- eval(dir)
fileNamesObject <- eval(fileNamesObject)
}
filesFullName <- list.files(dir, recursive = TRUE,
pattern = "\\.h5$", full.names = TRUE)
fileDir <- dirname(list.files(dir,
recursive = TRUE, pattern = "\\.h5$"))
filesProcessed <- basename(fileNamesObject)
# delete new files
newFilesIndex <- which(!basename(filesFullName) %in% filesProcessed)
if (length(newFilesIndex > 0)) {
filesFullName <- filesFullName[-newFilesIndex]
fileDir <- fileDir[-newFilesIndex]
}
fileName <- basename(filesFullName)
# test if there is subfolder if there is no subfolder
if (all(fileDir == ".")) {
sampleMetadata <- data.frame(row.names = fileName)
} else {
subfolder <- strsplit(fileDir, "/")
# max depth of subfolder
nSubfolder <- max(vapply(subfolder, length, 0))
# set every element at the same length
subfolder <- lapply(subfolder, function(x) {
length(x) <- nSubfolder
return(x)
})
group <- do.call(rbind, subfolder)
sampleMetadata <- data.frame(subfolder = group,
row.names = fileName,
stringsAsFactors = FALSE)
}
sampleMetadata$date <- Reduce(c, getDate(ptrset))
return(sampleMetadata)
}
#' get sampleMetadata of a ptrSet
#' @param set a ptrSet object
#' @return a data.frame
#' @export
#' @examples
#' dirRaw <- system.file("extdata/exhaledAir", package = "ptairData")
#' exhaledPtrset <- createPtrSet(dir=dirRaw, setName="exhaledPtrset",
#' mzCalibRef = c(21.022, 60.0525), fracMaxTIC = 0.7, saveDir = NULL)
#' SMD<-getSampleMetadata(exhaledPtrset)
getSampleMetadata <- function(set) {
if (!methods::is(set, "ptrSet"))
stop("set is not a ptrSet object")
sampleMetadata <- set@sampleMetadata
return(sampleMetadata)
}
#' set sampleMetadata in a ptrSet
#'
#' Insert a samplemetada data.frame in a ptrSet object. The dataframe must
#' have all file names in rownames.
#' @param set a ptrSet object
#' @param sampleMetadata a data.frame with all file names of the ptrSet in
#' row names
#' @return the ptrSet object in argument with the sampleMetadata modified
#' @export
#' @examples
#' dirRaw <- system.file("extdata/exhaledAir", package = "ptairData")
#' exhaledPtrset <- createPtrSet(dir=dirRaw, setName="exhaledPtrset",
#' mzCalibRef = c(21.022, 60.0525), fracMaxTIC = 0.7, saveDir = NULL )
#' SMD<-getSampleMetadata(exhaledPtrset )
#' colnames(SMD)[1]<-'individual'
#' exhaledPtrset<-setSampleMetadata(exhaledPtrset ,SMD)
setSampleMetadata <- function(set, sampleMetadata) {
# check if set is ptrSet
if (!methods::is(set, "ptrSet"))
stop("set is not a ptrSet object")
# check if row names contains all files
files <- names(set@mzCalibRef)
fileName <- basename(files)
testFilesName <- fileName %in% row.names(sampleMetadata)
if (!all(testFilesName)) {
stop(fileName[!testFilesName], "not in sampleMetadata
row names, please complete theme \n")
}
set@sampleMetadata <- sampleMetadata
## save
if (!is.null(getParameters(set)$saveDir)) {
changeName <- parse(text = paste0(getParameters(set)$name, "<- set "))
eval(changeName)
eval(parse(text = paste0("save(", getParameters(set)$name, ",file=
paste0( getParameters(set)$saveDir,'/', '",
getParameters(set)$name, ".RData '))")))
}
return(set)
}
#' export sampleMetadata
#' @param set a ptrSet object
#' @param saveFile a file path in tsv extension where the data.frame will be
#' exported
#' @return nothing
#' @export
#' @examples
#' dirRaw <- system.file("extdata/exhaledAir", package = "ptairData")
#' exhaledPtrset <- createPtrSet(dir=dirRaw, setName="exhaledPtrset",
#' mzCalibRef = c(21.022, 60.0525), fracMaxTIC = 0.7, saveDir = NULL )
#' saveFile<-file.path(getwd(),'sampleMetadata.tsv')
#' #exportSampleMetada(exhaledPtrset ,saveFile)
exportSampleMetada <- function(set, saveFile) {
if (!methods::is(set, "ptrSet"))
stop("set is not a ptrSet object")
if (tools::file_ext(saveFile) != "tsv")
stop("saveFile must be in
.tsv extension")
# get sampleMetadata
sampleMetadata <- getSampleMetadata(set)
# write the data.frame in a tsv file
utils::write.table(sampleMetadata, file = saveFile, sep = "\t",
col.names = NA)
}
#' import a sampleMetadata from a tsv file to a ptrSet object
#' @param set a ptrSet object
#' @param file a tsv file contains the sample metada to import, with all
#' file names in row name
#' (the fist column on th excel).
#' @return a ptrSet with th enew sample Metadata
#' @export
#' @examples
#' library(ptairData)
#' dirRaw <- system.file("extdata/exhaledAir", package = "ptairData")
#' exhaledPtrset <- createPtrSet(dir=dirRaw, setName="exhaledPtrset",
#' mzCalibRef = c(21.022, 60.0525), fracMaxTIC = 0.7, saveDir = NULL )
#' saveFile<-file.path(getwd(),'sampleMetadata.tsv')
#' #exportSampleMetada(exhaledPtrset ,saveFile)
#' #exhaledPtrset<-importSampleMetadata(exhaledPtrset ,saveFile)
importSampleMetadata <- function(set, file) {
if (!methods::is(set, "ptrSet"))
stop("set is not a ptrSet object")
sampleMetadata <- try(utils::read.table(file = file, sep = "\t",
header = TRUE,
row.names = 1, quote = ""))
# check if row names contains all files
dir <- getParameters(set)$dir
if (methods::is(dir, "expression"))
dir <- eval(dir)
files <- list.files(dir, recursive = TRUE, pattern = "\\.h5$")
fileName <- basename(files)
testFilesName <- fileName %in% row.names(sampleMetadata)
if (!all(testFilesName)) {
stop(fileName[!testFilesName], "not in sampleMetadata
row names, please complete theme \n")
}
set<- setSampleMetadata(set = set,sampleMetadata = sampleMetadata)
return(set)
}
## time limits ----
#' @rdname timeLimits
#' @export
methods::setMethod(f = "timeLimits", signature = "ptrSet",
function(object, fracMaxTIC = 0.5,fracMaxTICBg = 0.2,
derivThresholdExp = 1, derivThresholdBg = 0.05,
minPoints = 2,
degreeBaseline = 1, baseline = TRUE,
redefineKnots = TRUE, plotDel = FALSE) {
fileNames <- getParameters(object)$listFile
if (methods::is(fileNames, "expression"))
fileNames <- eval(fileNames)
fileNames <- basename(fileNames)
for (file in fileNames) {
TIC <- getTimeInfo(object)$breathTracer[[file]]
indLim <- timeLimitFun(TIC, fracMaxTIC = fracMaxTIC,
fracMaxTICBg = fracMaxTICBg,
derivThresholdExp = derivThresholdExp,
derivThresholdBg = derivThresholdBg,
degreeBaseline = degreeBaseline,
baseline = baseline,
minPoints = minPoints,
plotDel = plotDel)
object<-setTimeLimits(object,newTimeLimites =indLim,index = file )
}
paramterTimeLimit <- list(fracMaxTIC = fracMaxTIC,
fracMaxTICBg = fracMaxTICBg,
derivThresholdExp = derivThresholdExp,
derivThresholdBg = derivThresholdBg,
degreeBaseline = degreeBaseline,
minPoints = minPoints)
parameter<-getParameters(object)
parameter$timeLimit <- paramterTimeLimit
object<-setParameters(object,parameter)
if (redefineKnots)
object <- defineKnots(object, knotsPeriod = getParameters(object)$knotsPeriod)
saveDir <- getParameters(object)$saveDir
objName <- getParameters(object)$name
if (!is.null(saveDir)) {
if (!is.null(objName)) {
changeName <- parse(text = paste0(objName, "<- object "))
eval(changeName)
eval(parse(text = paste0("save(", objName, ",file= paste0(
saveDir,'/', '",
objName, ".RData '))")))
} else save(object, file = paste0(saveDir, "/ptrSet.RData"))
}
return(object)
})
## define knots-----
#' Define the knots location
#'
#' \code{defineKnots} function determine the knots location for a ptrSet or
#' ptrRaw object.
#' There is three possibilities :
#' \itemize{
#' \item \code{method = expiration} in the expiration periods, a knots is
#' placed every \code{knotsPeriod} seconds, and 1 knots in the middle of two
#' expiration, one at begin and at the end
#' \item \code{method= uniforme}, the knots are placed uniformly every
#' \code{knotsPeriod} time points
#' \item give in \code{knotsList} a list of knot, with all base name file in
#' name of the list element. All file must be informed. The knots location
#' must be contained in the time axis
#' }
#' @param object ptrSet object
#' @param knotsPeriod the period in second (times scale) between two knots for
#' the two dimensional modelization
#' @param method expiration or uniform
#' @param knotsList a list of knot location for each files, with all base name
#' file in name of the list element
#' @return a list with numeric vector of knots for each file
#' @examples
#' library(ptairData)
#' dirRaw <- system.file("extdata/exhaledAir", package = "ptairData")
#' exhaledPtrset <- createPtrSet(dir=dirRaw, setName="exhaledPtrset",
#' mzCalibRef = c(21.022, 60.0525), fracMaxTIC = 0.7, saveDir = NULL )
#'
#' #### placed knots every 2 times points
#' exhaledPtrset <- defineKnots(exhaledPtrset ,knotsPeriod=2,method='uniform')
#'
#' #### placed knots every 3 times points in the expiration (default)
#' exhaledPtrset <- defineKnots(exhaledPtrset ,knotsPeriod=3,method='expiration')
#' @rdname defineKnots
#' @export
methods::setMethod(f = "defineKnots", signature = "ptrSet",
function(object, knotsPeriod = 3,
method = c("expiration", "uniform")[1],
knotsList = NULL) {
if (knotsPeriod == 0) {
knots <- list(NULL)
} else {
if (!is.null(knotsList)) {
# check if all file are set
if (any(!names(knotsList) %in% names(getTimeInfo(object)$TIC)))
stop(names(knotsList)[!names(knotsList) %in% names(getTimeInfo(object)$TIC)],
" missing")
# check if interior knots are in the times axis
test <- apply(names(knotsList), function(file) {
knots <- knotsList[[file]]
t <- as.numeric(names(getTimeInfo(object)$TIC[[file]]))
return(knots[1] >= t[1] & utils::tail(knots, 1) <= utils::tail(t,
1))
})
if (any(!test))
stop("knots are not contained in
the time axis for file ",
names(knotsList)[which(!test)], "\n")
knots <- knotsList
} else {
knots <- lapply(names(getTimeInfo(object)$TIC), function(file) {
background <- getTimeInfo(object)$timeLimit[[file]]$backGround
t <- as.numeric(names(getTimeInfo(object)$TIC[[file]]))
res <- try(defineKnotsFunc(t, background, knotsPeriod, method, file))
if (!is.null(attr(res, "condition")))
res <- NULL
res
})
}
names(knots) <- names(getTimeInfo(object)$TIC)
}
object <- setKnots(object,knots)
parameter<-getParameters(object)
saveDir <- parameter$saveDir
objName <- parameter$name
parameter$knotsPeriod <- knotsPeriod
object<-setParameters(object,parameter)
if (!is.null(saveDir)) {
if (!is.null(objName)) {
changeName <- parse(text = paste0(objName, "<- object "))
eval(changeName)
eval(parse(text = paste0("save(", objName, ",file=
paste0( saveDir,'/', '",
objName, ".RData '))")))
} else save(object, file = paste0(saveDir, "/ptrSet.RData"))
}
return(object)
})
## calibration----
#' @param fileNames file to recalibrate
#' @param ... " "
#' @rdname calibration
#' @export
methods::setMethod(f = "calibration", signature = "ptrSet",
function(x, mzCalibRef = c(21.022,
29.013424, 41.03858, 75.04406, 203.943, 330.8495),
calibrationPeriod = 60, tol = 70,fileNames=getParameters(x)$listFile) {
if (methods::is(fileNames, "expression"))
fileNames <- eval(fileNames)
for (file in fileNames) {
raw <- readRaw(file, mzCalibRef = mzCalibRef,
calibrationPeriod = calibrationPeriod)
calibSpectr <- alignCalibrationPeak(getCalibrationInfo(raw)$calibSpectr,
getCalibrationInfo(raw)$calibMassRef,
length(getRawInfo(raw)$time))
resolutionEstimated <- estimateResol(calibMassRef = getCalibrationInfo(raw)$calibMassRef,
calibSpectr = calibSpectr)
calibList<-list(coefCalib=getCalibrationInfo(raw)$calibCoef,
mzCalibRef=getCalibrationInfo(raw)$calibMassRef,
errorCalibPpm=getCalibrationInfo(raw)$calibError,
signalCalibRef=calibSpectr,
resolution=resolutionEstimated,
indexTimeCalib=raw@indexTimeCalib)
x<-setCalibration(x,calibList,index=basename(file))
x@peakShape[[basename(file)]]<-raw@peakShape
}
parameter<-getParameters(x)
if(length(fileNames) == length(x@mzCalibRef)){
parameter$mzCalibRef <- mzCalibRef
x<-setParameters(x,parameter)
}
saveDir <- parameter$saveDir
objName <- parameter$name
if (!is.null(saveDir)) {
if (!is.null(objName)) {
changeName <- parse(text = paste0(objName, "<- x "))
eval(changeName)
eval(parse(text = paste0("save(", objName, ",file= paste0(
saveDir,'/', '",
objName, ".RData '))")))
} else save(x, file = paste0(saveDir, "/ptrSet.RData"))
}
return(x)
})
## other ----
#' get the files directory of a ptrSet
#' @param ptrSet ptrSet object
#' @return the directory in absolute path as character
#' @examples
#' library(ptairData)
#' dirRaw <- system.file("extdata/exhaledAir", package = "ptairData")
#' exhaledPtrset <- createPtrSet(dir=dirRaw, setName="exhaledPtrset",
#' mzCalibRef = c(21.022, 60.0525), fracMaxTIC = 0.7, saveDir = NULL )
#' getDirectory(exhaledPtrset )
#' @export
getDirectory <- function(ptrSet) {
dir <- getParameters(ptrSet)$dir
if (methods::is(dir, "expression"))
dir <- eval(dir)
return(dir)
}
#' remove the peakList of an ptrSet object
#'
#' This function is useful when you want to change the parameters of the detect
#' peak function. First delete the peakList with \code{rmPeakList}, and apply
#' \code{detectPeak}with the new parameters.
#' @param object ptrSet object
#' @return a ptrSet
#' @export
#' @examples
#' library(ptairData)
#' dirRaw <- system.file("extdata/exhaledAir", package = "ptairData")
#' exhaledPtrset <- createPtrSet(dir=dirRaw, setName="exhaledPtrset",
#' mzCalibRef = c(21.022, 60.0525), fracMaxTIC = 0.7, saveDir = NULL )
#' exhaledPtrset <-rmPeakList(exhaledPtrset )
rmPeakList <- function(object) {
object@peakList <- list()
parameter<-getParameters(object)
parameter$detectPeakParam <- NULL
object<-setParameters(object,parameter)
return(object)
}
#' show a ptrSet object
#'
#' It indicates the directory, the numbre of files that contain the directory at
#' the moment, and the number of processed files.
#' The two numbers are different, use \code{updatePtrSet} function.
#' @param object a ptrSet object
#' @return nothing
#' @export
methods::setMethod("show", "ptrSet", function(object) {
dir <- getParameters(object)$dir
fileCheck <- getTimeInfo(object)$TIC
if (methods::is(dir, "expression")) {
dir <- eval(dir)
}
nFiles <- length(list.files(dir, recursive = TRUE, pattern = "\\.h5$"))
nFilesCheck <- length(fileCheck)
nFilesProcess <- length( getPeakList(object) )
cat("ptrSet object :", getParameters(object)$name, "\n")
cat("directory:", dir, "\n")
cat(" ", nFiles, "files contained in the directory \n")
cat(" ", nFilesCheck, "files checked", "\n")
cat(" ", nFilesProcess, "files processed by detectPeak", "\n")
if (!is.null(getParameters(object)$saveDir))
cat("object save in ", getParameters(object)$saveDir)
})
#' get the peak list of a ptrSet object
#' @param set ptrSet object
#' @return a list of expressionSet, where:
#'\itemize{
#'\item assay Data contains the matrix with m/z peak center in row names, and the
#'quantification in cps at each time point
#'\item feature Data the matrix with m/z peak center in row names, and the
#'following columns:
#' \itemize{
#' \item quanti_unit: the mean of the quantification over the expiration/headspace
#'time limits defined
#' \item background_unit: the mean of the quantification over the background time
#'limits defined
#' \item diffAbs_unit: the mean of the quantification over the expiration/headspace
#'time limits defined after subtracting the baseline estimated from the background
#'points defined
#' \item pValLess/ pValGreater: the p-value of the unilateral t-test, who test
#'that quantification (in cps) of expiration points are less/greater than the
#'intensity of the background.
#' \item lower/upper: integration boundaries
#' \item parameter peak: estimated peak parameter
#' }
#'}
#' @examples
#' library(ptairData)
#' dirRaw <- system.file("extdata/exhaledAir", package = "ptairData")
#' exhaledPtrset <- createPtrSet(dir=dirRaw, setName="exhaledPtrset",
#' mzCalibRef = c(21.022, 60.0525), fracMaxTIC = 0.7, saveDir = NULL )
#' exhaledPtrset<-detectPeak(exhaledPtrset,mz=c(21,59))
#' peakList<- getPeakList(exhaledPtrset )
#' X<-Biobase::exprs(peakList[[1]])
#' Y<- Biobase::fData(peakList[[1]])
#' head(Y)
#' @export
getPeakList <- function(set) {
return(set@peakList)
}
#' get the file names containing in the directory of a ptrSet or ptrRaw
#' @param object ptrSet object
#' @param fullNames logical: if \code{TRUE}, it return the the directory path is
#' prepended to the file names.
#' @return a vector of character that contains all file names
#' @examples
#' library(ptairData)
#' dirRaw <- system.file("extdata/exhaledAir", package = "ptairData")
#' exhaledPtrset <- createPtrSet(dir=dirRaw, setName="exhaledPtrset",
#' mzCalibRef = c(21.022, 60.0525), fracMaxTIC = 0.7, saveDir = NULL )
#' getFileNames(exhaledPtrset )
#' @rdname getFileNames
#' @export
methods::setMethod("getFileNames", signature = "ptrSet",
function(object, fullNames) {
fileFullNames <- getParameters(object)$listFile
if (methods::is(fileFullNames, "expression"))
fileFullNames <- eval(fileFullNames)
if (fullNames)
return(fileFullNames) else return(basename(fileFullNames))
})
#' get the file names containing in the directory of a ptrSet
#' @rdname getFileNames
#' @export
methods::setMethod("getFileNames", signature = "ptrRaw",
function(object, fullNames) {
return(object@name)
})
methods::setMethod("getParameters", signature = "ptrSet",
function(object) {
return(object@parameter)
})
methods::setMethod("setParameters", signature = "ptrSet",
function(object,newparam) {
object@parameter<-newparam
return(object)
})
methods::setMethod("getDate", signature = "ptrSet",
function(object) {
return(object@date)
})
methods::setMethod("getDate", signature = "ptrRaw",
function(object) {
return(object@date)
})
methods::setMethod("getCalibrationInfo", signature = "ptrSet",
function(object) {
return(list(mzCalibRef=object@mzCalibRef,
signalCalibRef = object@signalCalibRef,
errorCalibPpm =object@errorCalibPpm,
coefCalib=object@coefCalib,
indexTimeCalib=object@indexTimeCalib)
)
})
methods::setMethod("getCalibrationInfo", signature = "ptrRaw",
function(object) {
return(list(
calibMzToTof = object@calibMzToTof,
calibToftoMz = object@calibToftoMz,
calibCoef=object@calibCoef,
calibMassRef = object@calibMassRef,
calibError=object@calibError,
calibSpectr = object@calibSpectr,
indexTimeCalib=object@indexTimeCalib)
)
})
methods::setMethod("setCalibration", signature = "ptrRaw",
function(object,newCalibrationInfo) {
object@indexTimeCalib <- newCalibrationInfo$indexTimeCalib
object@calibSpectr <- newCalibrationInfo$calibSpectr
object@calibError <- newCalibrationInfo$calibError
object@calibCoef <-newCalibrationInfo$calibCoef
return(object)
})
methods::setMethod("setCalibration", signature = "ptrSet",
function(object,newCalibrationInfo,index=NULL) {
if(is.null(index)) index<- seq_along(newCalibrationInfo)
if( length(index == 1) ){
object@coefCalib[[index]] <- newCalibrationInfo$coefCalib
object@mzCalibRef[[index]] <- newCalibrationInfo$mzCalibRef
object@errorCalibPpm[[index]] <- newCalibrationInfo$errorCalibPpm
object@signalCalibRef[[index]] <- newCalibrationInfo$signalCalibRef
object@resolution[[index]] <- newCalibrationInfo$resolution
object@indexTimeCalib[[index]] <-newCalibrationInfo$indexTimeCalib
} else{
object@coefCalib[index] <- newCalibrationInfo$coefCalib
object@mzCalibRef[index] <- newCalibrationInfo$mzCalibRef
object@errorCalibPpm[index] <- newCalibrationInfo$errorCalibPpm
object@signalCalibRef[index] <- newCalibrationInfo$signalCalibRef
object@resolution[index] <- newCalibrationInfo$resolution
object@indexTimeCalib[index] <-newCalibrationInfo$indexTimeCalib
}
return(object)
})
methods::setMethod("getPeaksInfo", signature = "ptrSet",
function(object) {
return(list(knots=object@knots,
fctFit=object@fctFit,
peakShape=object@peakShape))
})
methods::setMethod("getPeaksInfo", signature = "ptrRaw",
function(object) {
return(list(
peakShape=object@peakShape))
})
methods::setMethod("getTimeInfo", signature = "ptrSet",
function(object) {
return(list(TIC =object@TIC,
breathTracer=object@breathTracer,
timeLimit=object@timeLimit))
})
methods::setMethod("setTimeLimits", signature = "ptrSet",
function(object,newTimeLimites,index=NULL) {
if(is.null(index)){
object@timeLimit <- newTimeLimites
} else if(length(index) == 1 ) {
object@timeLimit[[index]] <- newTimeLimites
} else object@timeLimit[index] <- newTimeLimites
return(object)
})
methods::setMethod("setKnots", signature = "ptrSet",
function(object,newKnots,index=NULL) {
if(is.null(index)){
object@knots <- newKnots
} else if(length(index) == 1 ) {
object@knots[[index]] <- newKnots
} else object@knots[index] <- newKnots
return(object)
})
methods::setMethod("getPTRInfo", signature = "ptrSet",
function(object) {
return(list(primaryIon=object@primaryIon,
ptrTransmisison = object@ptrTransmisison,
prtReaction = object@prtReaction,
resolution = object@resolution))
})
methods::setMethod("getPTRInfo", signature = "ptrRaw",
function(object) {
return(list(
ptrTransmisison = object@ptrTransmisison,
prtReaction = object@prtReaction))
})
methods::setMethod("getRawInfo", signature = "ptrRaw",
function(object) {
return(list(
rawM= object@rawM,
mz = object@mz,
time = object@time))
})
methods::setMethod("setPeakShape", signature = "ptrRaw",
function(object,peakShape) {
object@peakShape<-peakShape
object
})
methods::setMethod("deleteFilePtrSet", signature = "ptrSet",
function(object,deletedFiles) {
object@mzCalibRef[deletedFiles] <- NULL
object@signalCalibRef[deletedFiles] <- NULL
object@errorCalibPpm[deletedFiles] <- NULL
object@coefCalib[deletedFiles] <- NULL
object@indexTimeCalib[deletedFiles] <- NULL
object@knots[deletedFiles] <- NULL
object@fctFit[deletedFiles] <- NULL
object@peakShape[deletedFiles] <- NULL
object@primaryIon[deletedFiles] <- NULL
object@resolution[deletedFiles] <- NULL
object@TIC[deletedFiles] <- NULL
object@timeLimit[deletedFiles] <- NULL
object@peakList[deletedFiles] <- NULL
object@breathTracer[deletedFiles] <- NULL
object@ptrTransmisison[deletedFiles] <- NULL
object@prtReaction[deletedFiles] <- NULL
object@date[deletedFiles] <- NULL
object
})
methods::setMethod("mergedPtrSet", signature = "ptrSet",
function(object,ptrSetNewfile,orderedFile) {
object@mzCalibRef <- c(getCalibrationInfo(object)$mzCalibRef, ptrSetNewfile$mzCalibRefList)[basename(orderedFile)]
object@signalCalibRef <- c(getCalibrationInfo(object)$signalCalibRef, ptrSetNewfile$signalCalibRef)[basename(orderedFile)]
object@errorCalibPpm <- c(getCalibrationInfo(object)$errorCalibPpm, ptrSetNewfile$errorCalibPpm)[basename(orderedFile)]
object@coefCalib <- c(getCalibrationInfo(object)$coefCalib, ptrSetNewfile$coefCalibList)[basename(orderedFile)]
object@resolution <- c(getPTRInfo(object)$resolution, ptrSetNewfile$resolution)[basename(orderedFile)]
object@TIC <- c(getTimeInfo(object)$TIC, ptrSetNewfile$TIC)[basename(orderedFile)]
object@timeLimit <- c(getTimeInfo(object)$timeLimit, ptrSetNewfile$timeLimit)[basename(orderedFile)]
object@primaryIon <- c(getPTRInfo(object)$primaryIon, ptrSetNewfile$primaryIon)[basename(orderedFile)]
object@breathTracer <- c(getTimeInfo(object)$breathTracer, ptrSetNewfile$breathTracer)[basename(orderedFile)]
object@ptrTransmisison <- c(getPTRInfo(object)$ptrTransmisison, ptrSetNewfile$ptrTransmisison)[basename(orderedFile)]
object@prtReaction <- c(getPTRInfo(object)$prtReaction, ptrSetNewfile$prtReaction)[basename(orderedFile)]
object@date <- c(getDate(object), ptrSetNewfile$date)[basename(orderedFile)]
object@peakShape <- c(object@peakShape, ptrSetNewfile$peakShape)[basename(orderedFile)]
object@knots<-c(object@knots, ptrSetNewfile$knots)[basename(orderedFile)]
object@indexTimeCalib<-c(object@indexTimeCalib, ptrSetNewfile$indexTimeCalib)[basename(orderedFile)]
object@fctFit<-c(object@fctFit, ptrSetNewfile$fctFit)[basename(orderedFile)]
return(object)
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.