Nothing
R/erah.R
In erah: Automated Spectral Deconvolution, Alignment, and Metabolite Identification in GC/MS-Based Untargeted Metabolomics
Defines functions
scansPerSecond
processSample
setAlPar
setDecPar
Documented in
setAlPar
setDecPar
globalVariables("mslib")
#' @useDynLib erah
## Main Software functions:
#' @name setDecPar
#' @aliases setDecPar
#' @title Set Software Parameters
#' @description Sets Software Parameters for eRah.
#' @param min.peak.width Minimum compound peak width (in seconds). This is a critical parameter that conditions the efficiency of eRah. Typically, this should be the half of the mean compound width.
#' @param min.peak.height Minimum compound peak height
#' @param noise.threshold Data above this threshold will be considered as noise
#' @param avoid.processing.mz The masses that do not want to be considered for processing. Typically, in GC-MS those masses are 73,74,75,147,148 and 149, since they are they are ubiquitous mass fragments typically generated from compounds carrying a trimethylsilyl moiety.
#' @param compression.coef Data is compressed when using the orthogonal signal deconvolution (OSD) algorithm according to this value. A level 2 of compression is recomended.
#' @param analysis.time The chromatographic retention time window to process. If 0, all the chromatogram is processed.
#' @details See eRah vignette for more details. To open the vignette, execute the following code in R:
#' vignette("eRahManual", package="erah")
#' @references [1] Xavier Domingo-Almenara, et al., eRah: A Computational Tool Integrating Spectral Deconvolution and Alignment with Quantification and Identification of Metabolites in GC-MS-Based Metabolomics. Analytical Chemistry (2016). DOI: 10.1021/acs.analchem.6b02927
#' @author Xavier Domingo-Almenara. xavier.domingo@urv.cat
#' @seealso \code{\link{newExp}} \code{\link{deconvolveComp}} \code{\link{alignComp}} \code{\link{setAlPar}}
#' @examples \dontrun{
#' # The following will set eRah for analyzing the chromatograms
#' #from minutes 5 to 15, and withouth taking into account the masses
#' #35:69,73:75,147:149, widht a minimum peak widht of 0.7 seconds.
#' ex.dec.par <- setDecPar(min.peak.width = 0.7,
#' min.peak.height = 5000,
#' noise.threshold = 500,
#' avoid.processing.mz = c(35:69,73:75,147:149),
#' analysis.time = c(5,15))
#' }
#' @export
setDecPar <- function(min.peak.width,
min.peak.height=2500,
noise.threshold=500,
avoid.processing.mz=c(73:75,147:149),
compression.coef=2,
analysis.time=0){
softPar <- new("eRahSoftParameters",algorithm="eRah-OSD", min.peak.width = min.peak.width/60, min.peak.height = min.peak.height, noise.threshold = noise.threshold, avoid.processing.mz = avoid.processing.mz, compression.coef = compression.coef, analysis.time=analysis.time)
softPar
}
#' @name setAlPar
#' @aliases setAlPar
#' @title Set Alignment Parameters
#' @description Setting alignment parameters for eRah.
#' @usage setAlPar(min.spectra.cor, max.time.dist,mz.range = c(70:600))
#' @param min.spectra.cor Minimum spectral correlation value. From 0 (non similar) to 1 (very similar). This value sets how similar two or more compounds have be to be considered for alignment between them.
#' @param max.time.dist Maximum retention time distance. This value (in seconds) sets how far two or more compounds can be to be considered for alignment between them.
#' @param mz.range The range of masses that is considered when comparing spectra.
#' @references [1] Xavier Domingo-Almenara, et al., eRah: A Computational Tool Integrating Spectral Deconvolution and Alignment with Quantification and Identification of Metabolites in GC-MS-Based Metabolomics. Analytical Chemistry (2016). DOI: 10.1021/acs.analchem.6b02927
#' @author Xavier Domingo-Almenara. xavier.domingo@urv.cat
#' @seealso \code{\link{newExp}} \code{\link{setDecPar}} \code{\link{alignComp}}
#' @examples\dontrun{
#' # The following will set eRah for aligning compounds which are
#' # at least 90 (per cent) similar, and which peaks are at a
#' # maximum distance of 2 seconds. All the masses are considered when
#' # computing the spectral similarity.
#'
#' ex.al.par <- setAlPar(min.spectra.cor=0.90, max.time.dist=2,
#' mz.range=1:600)
#' }
#' @export
setAlPar <- function(min.spectra.cor, max.time.dist, mz.range=c(70:600))
{
alPar <- new("eRahAlParameters", algorithm="eRah", min.spectra.cor=min.spectra.cor, max.time.dist=max.time.dist/60, mz.range = mz.range, method="eRah")
alPar
}
#' @name newExp
#' @aliases newExp
#' @title New Experiment
#' @description Sets a new experiment for eRah
#' @param instrumental A data.frame containing the sample instrumental information.
#' @param phenotype (optional) A data.frame containing sample phenotype information.
#' @param info Experiment description
#' @details See eRah vignette for more details. To open the vignette, execute the following code in R:
#' vignette("eRahManual", package="erah")
#' @return \code{newExp} returns an S4 object of the class 'MetaboSet'.
#' @references [1] Xavier Domingo-Almenara, et al., eRah: A Computational Tool Integrating Spectral Deconvolution and Alignment with Quantification and Identification of Metabolites in GC-MS-Based Metabolomics. Analytical Chemistry (2016). DOI: 10.1021/acs.analchem.6b02927
#' @author Xavier Domingo-Almenara. xavier.domingo@urv.cat
#' @examples \dontrun{
#' library(gcspikelite)
#' data(targets)
#'
#' files <- list.files(system.file('data',package = 'gcspikelite'),full.names = TRUE)
#' files <- files[sapply(files,grepl,pattern = 'CDF')]
#'
#' instrumental <- createInstrumentalTable(files)
#' phenotype <- createPhenoTable(files,as.character(targets$Group[order(targets$FileName)]))
#'
#' ex <- newExp(instrumental = instrumental,
#' phenotype = phenotype, info = "DEMO Experiment")
#' }
#' @seealso \code{\link{createInstrumentalTable}} \code{\link{createPhenoTable}} \code{\link{setDecPar}} \code{\link{setAlPar}}
#' @importFrom utils read.csv
#' @export
newExp<- function (instrumental, phenotype = NULL, info = character())
{
if (is.null(phenotype)) {
if (is.data.frame(instrumental)){
instrumental <- instrumental
}
else if (is.character(instrumental) && file.exists(instrumental)){
instrumental<-read.csv(instrumental, sep=";")
}
else stop("Instrumental file format not valid.")
phenotype = as.data.frame(NULL)
factors.list <- lapply(1:nrow(instrumental), function(x) {
as.data.frame(NULL)
})
names(factors.list) <- as.vector(instrumental$sampleID)
ident.list <- as.data.frame(matrix(0, ncol = 7, dimnames = list(row = 0,
col = c("AlignID", "tmean", "Name", "MatchFactor", "CAS",
"Formula", "DB.Id"))))
uni.stats <- as.data.frame(matrix(0, ncol = 3, dimnames = list(row = 0,
col = c("Id", "FoldChangue", "pvalue"))))
multi.stats <- as.data.frame(matrix(0, ncol = 3, dimnames = list(row = 0,
col = c("Id", "CompoundsInvolved", "pvalue"))))
al.par <- list()
id.par <- list()
soft.par <- list()
stat.parameters <- new("MSResultsParameters", Alignment = al.par,
Identification = id.par)
statistics <- new("Statistics", Univariate = uni.stats,
Multivariate = multi.stats)
MS.Results <- new("Results", Parameters = stat.parameters,
Identification = ident.list, Statistics = statistics)
MS.Data <- new("Data", FeatureList = list(NULL), FactorList = factors.list,
Parameters = list(NULL))
MS.MetaData <- new("MetaData", Instrumental = instrumental,
Phenotype = phenotype, DataDirectory = "")
col.correct <- c("sampleID", "filename", "date", "time")
for (i in 1:length(col.correct)) if (length(apply(as.matrix(colnames(MS.MetaData@Instrumental)),
1, function(x) grep(col.correct[i], x))) == 0)
stop("Invalid instrumental file. The file must contain at least the following columns: ",
paste(col.correct, collapse = ", "))
if (!is.null(MS.MetaData@Phenotype))
col.correct <- c("sampleID", "class")
for (i in 1:length(col.correct)) if (length(apply(as.matrix(colnames(MS.MetaData@Phenotype)),
1, function(x) grep(col.correct[i], x))) == 0)
sample.container <- new("MetaboSet", Info = info, Data = MS.Data,
MetaData = MS.MetaData, Results = MS.Results)
sample.container
} else if (!is.null(phenotype)){
if (is.data.frame(instrumental)){
instrumental <- instrumental
}
else if (is.character(instrumental) && file.exists(instrumental)){
instrumental<-read.csv(instrumental, sep=";")
}
else stop("Instrumental file format not valid.")
if (is.data.frame(phenotype)){
phenotype <- phenotype
}
else if (is.character(phenotype) && file.exists(phenotype)){
phenotype <- read.csv(phenotype, sep=";")
}
else stop("Phenotype file format not valid.")
factors.list <- lapply(1:nrow(instrumental), function(x) {
as.data.frame(NULL)
})
names(factors.list) <- as.vector(instrumental$sampleID)
ident.list <- as.data.frame(matrix(0, ncol = 7, dimnames = list(row = 0,
col = c("AlignID", "tmean", "Name", "MatchFactor", "CAS",
"Formula", "DB.Id"))))
uni.stats <- as.data.frame(matrix(0, ncol = 3, dimnames = list(row = 0,
col = c("Id", "FoldChangue", "pvalue"))))
multi.stats <- as.data.frame(matrix(0, ncol = 3, dimnames = list(row = 0,
col = c("Id", "CompoundsInvolved", "pvalue"))))
al.par <- list()
id.par <- list()
soft.par <- list()
stat.parameters <- new("MSResultsParameters", Alignment = al.par,
Identification = id.par)
statistics <- new("Statistics", Univariate = uni.stats,
Multivariate = multi.stats)
MS.Results <- new("Results", Parameters = stat.parameters,
Identification = ident.list, Statistics = statistics)
MS.Data <- new("Data", FeatureList = list(NULL), FactorList = factors.list,
Parameters = list(NULL))
MS.MetaData <- new("MetaData", Instrumental = instrumental,
Phenotype = phenotype, DataDirectory = "")
col.correct <- c("sampleID", "filename", "date", "time")
for (i in 1:length(col.correct)) if (length(apply(as.matrix(colnames(MS.MetaData@Instrumental)),
1, function(x) grep(col.correct[i], x))) == 0)
stop("Invalid instrumental file. The file must contain at least the following columns: ",
paste(col.correct, collapse = ", "))
if (!is.null(MS.MetaData@Phenotype)) {
col.correct <- c("sampleID", "class")
for (i in 1:length(col.correct)) if (length(apply(as.matrix(colnames(MS.MetaData@Phenotype)),
1, function(x) grep(col.correct[i], x))) == 0)
stop("Invalid phenotype file. The file must contain at least the following columns: ",
paste(col.correct, collapse = ", "))
}
sample.container <- new("MetaboSet", Info = info, Data = MS.Data,
MetaData = MS.MetaData, Results = MS.Results)
sample.container
}
}
#' @rdname deconvolveComp
#' @title Deconvolution of compounds in samples
#' @description Deconvolution of GC-MS data
#' @param Experiment A 'MetaboSet' S4 object containing the experiment data previously created by newExp.
#' @param decParameters The software deconvolution parameters object previously created by setDecPar
#' @param samples.to.process Vector indicating which samples are to be processed.
#' @param down.sample If TRUE, chromatograms are down sampled to define one peak with 10 scan points (according to the minimum peak width). This is to process longer chromatograms with wider peak widths (more than 20 seconds peak width and small scans per second values). See details.
#' @param virtualScansPerSecond A virtual scans per second. If chromatograms are downsampled (for example, for a 1 mean peak width a 1 scans per second sampling frequency was used), eRah could not perform as expected. In these cases, the BEST solution is to re-acquire the samples. However, by selecting a different (virtual) scans per second frequency, eRah can upsample the data and process it more effectively.
#' @details See eRah vignette for more details. To open the vignette, execute the following code in R:
#' vignette("eRahManual", package="erah")
#'
#' eRah uses multivariate methods which run-time performance depend on the amount of data to be analyzed. When peaks are wider and the #' scans per second is also a small value, the number of points (scans) that define a peak might be too many, leading eRah to a poor run#'-time performance. To solve that, use down.sample=TRUE to allow eRah to define a peak with 10 seconds, and analyze the data more #' efficiently.
#' @return The function returns an updated S4 'MetaboSet' class, where the GC-MS samples have been now deconvolved.
#' @references [1] Xavier Domingo-Almenara, et al., eRah: A Computational Tool Integrating Spectral Deconvolution and Alignment with Quantification and Identification of Metabolites in GC-MS-Based Metabolomics. Analytical Chemistry (2016). DOI: 10.1021/acs.analchem.6b02927
#' @author Xavier Domingo-Almenara. xavier.domingo@urv.cat
#' @seealso \code{\link{newExp}} \code{\link{setAlPar}}
#' @examples \dontrun{
#' # Deconvolve data from a created experiment by \code{\link{newExp}}.
#' # ex <- newExp(instrumental="path")
#'
#' # The following will set eRah for analyzing the chromatograms
#' # from minutes 5 to 15, and withouth taking into account the masses
#' # 35:69,73:75,147:149, with a minimum peak width of 0.7 seconds.
#'
#' ex.dec.par <- setDecPar(min.peak.width=0.7, min.peak.height=5000,
#' noise.threshold=500, avoid.processing.mz=c(35:69,73:75,147:149),
#' analysis.time=c(5,15))
#'
#' # An now deconvolve the compounds in the samples:
#' # ex <- deconvolveComp(ex, decParameters=ex.dec.par)
#' }
#' @export
setGeneric('deconvolveComp',function(Experiment, decParameters, samples.to.process=NULL, down.sample=FALSE, virtualScansPerSecond=NULL){
standardGeneric('deconvolveComp')
})
#' @rdname deconvolveComp
#' @importFrom furrr future_map furrr_options
setMethod('deconvolveComp',signature = 'MetaboSet',
function(Experiment, decParameters, samples.to.process=NULL, down.sample=FALSE, virtualScansPerSecond=NULL){
plotting=FALSE
Number.of.Samples <- nrow(Experiment@MetaData@Instrumental)
if(is.null(samples.to.process)) samples.to.process <- 1:Number.of.Samples
stopifnot(samples.to.process>=1, max(samples.to.process)<=Number.of.Samples, length(samples.to.process)<=Number.of.Samples)
soft.par <- list(min.peak.width = decParameters@min.peak.width, min.peak.height = decParameters@min.peak.height, noise.threshold = decParameters@noise.threshold, avoid.processing.mz = decParameters@avoid.processing.mz, compression.coef = decParameters@compression.coef, analysis.time = decParameters@analysis.time)
Experiment@Data@Parameters <- soft.par
Experiment@Data@FactorList <- future_map(samples.to.process,~{
k <- which(samples.to.process == .x)
message(paste("\n Deconvolving compounds from",
as.character(Experiment@MetaData@Instrumental$filename[k]),
"... Processing",
k,"/",
length(samples.to.process),"\n"))
processSample(Experiment,
.x,
plotting,
down.sample,
virtualScansPerSecond)
},
.options = furrr_options(seed = TRUE))
names(Experiment@Data@FactorList) <- metaData(Experiment)$sampleID
Experiment <- scansPerSecond(Experiment)
message("Compounds deconvolved")
Experiment
}
)
#' @rdname alignComp
#' @title Alignment of compounds
#' @description Alignment of GC-MS deconvolved compounds
#' @usage alignComp(Experiment, alParameters, blocks.size=NULL)
#' @param Experiment A 'MetaboSet' S4 object containing the experiment data previously created by newExp and deconvolved by deconvolveComp.
#' @param alParameters The software alignment parameters object previously created by setAlPar
#' @param blocks.size For experiment of more than 1000 samples, and depending on the computer, alignment can be conducted by block segmentation. See details.
#' @details See eRah vignette for more details. To open the vignette, execute the following code in R:
#' vignette("eRahManual", package="erah")
#'
#' For experiments containing more than 100 (Windows) or 1000 (Mac or Linux) samples (numbers depending on the computer resoures and sample type). In those cases alignment can be conducted by block segmentation. For an experiment of e.g. 1000 samples, the block.size can be set to 100, so the alignment will perform as multiple (ten) 100-samples experiments, to later align them into a single experiment.
#'
#' This parameter is designed to solve the typical problem that appear when aligning under Windows operating system: "Error: cannot allocate vector of size XX Gb". Such a problem will not appear with Mac or Linux, but several hours of computation are expected when aligning a large number of samples. Using block segmentation provides a greatly improved run-time performance.
#' @return The function returns an updated S4 'MetaboSet' class, where the GC-MS samples have been now aligned.
#' @references [1] Xavier Domingo-Almenara, et al., eRah: A Computational Tool Integrating Spectral Deconvolution and Alignment with Quantification and Identification of Metabolites in GC-MS-Based Metabolomics. Analytical Chemistry (2016). DOI: 10.1021/acs.analchem.6b02927
#' @author Xavier Domingo-Almenara. xavier.domingo@urv.cat
#' @seealso \code{\link{newExp}} \code{\link{setDecPar}} \code{\link{deconvolveComp}}
#' @export
setGeneric('alignComp',function(Experiment, alParameters, blocks.size=NULL){
standardGeneric('alignComp')
})
#' @rdname alignComp
setMethod('alignComp',signature = 'MetaboSet',
function(Experiment, alParameters, blocks.size=NULL){
al.par <- list(alignment.algorithm=alParameters@algorithm, min.spectra.cor=alParameters@min.spectra.cor, max.time.dist=alParameters@max.time.dist, mz.range=alParameters@mz.range)
Experiment@Results@Parameters@Alignment <- al.par
min.spectra.cor <- Experiment@Results@Parameters@Alignment$min.spectra.cor
max.time.dist <- Experiment@Results@Parameters@Alignment$max.time.dist
mz.range <- Experiment@Results@Parameters@Alignment$mz.range
maxMZ <- max(mz.range)
# Experiment@Data@FactorList <- align.factors(Experiment@Data@FactorList, min.spectra.cor, max.time.dist, maxMZ, mz.range)
# Experiment@Results@Alignment <- create.factorlist.table(Experiment)
if(is.null(blocks.size))
{
Experiment@Data@FactorList <- align.factors(Experiment@Data@FactorList, min.spectra.cor, max.time.dist, maxMZ, mz.range)
Experiment@Results@Alignment <- create.factorlist.table(Experiment)
#return(Experiment)
}else{
#blocks.size <- 15
max.mz <- maxMZ
Itrt <- length(Experiment@Data@FactorList)/blocks.size
sequs <- trunc(seq(1, length(Experiment@Data@FactorList), length.out=Itrt))
sequs[1] <- 0
corresponding.list <- list()
block.list <- list()
#i <- 1
for(i in 1:(length(sequs)-1))
{
cat("Aligning block ", i, " of ", length(sequs)-1, "... \n", sep="")
ghost.object <- Experiment
ghost.object@Data@FactorList <- Experiment@Data@FactorList[(sequs[i]+1):sequs[(i+1)]]
factors.list <- ghost.object@Data@FactorList
ghost.object@Data@FactorList <- align.factors(factors.list, min.spectra.cor, max.time.dist, max.mz, mz.range)
ghost.factors.list <- create.factorlist.table(ghost.object)
block.list[[i]] <- data.frame(ID=ghost.factors.list$AlignID, RT=ghost.factors.list$tmean, Spectra=ghost.factors.list$Spectra)
corresponding.list <- c(corresponding.list,lapply(ghost.object@Data@FactorList, function(x) x$AlignID))
}
cat("Aligning factors across blocks... \n")
full.factorlist <- align.factors(block.list, min.spectra.cor, max.time.dist, max.mz, mz.range)
#MaxALID <- max(unlist(lapply(full.factorlist, function(x) x$AlignID)))
factors.list <- Experiment@Data@FactorList
if(!(any(unlist(lapply(factors.list,function(x) {is.null(x$AlignID)}))==FALSE)))
{
factors.list <- lapply(factors.list, function(x){
outp <- cbind(x,matrix(0,nrow=length(x$ID)))
colnames(outp)[ncol(outp)] <- "AlignID"
outp
})
}else{
factors.list <- lapply(factors.list, function(x){
x$AlignID <- rep(0,length(x$ID))
x
})
}
Experiment@Data@FactorList <- factors.list
free.aligned.slots <- list()
for(i in 1:length(full.factorlist))
{
for(j in (sequs[i]+1):sequs[(i+1)])
{
ID.vct <- sapply(full.factorlist[[i]]$ID, function(x) {x.num <- which(corresponding.list[[j]]==x)
if(length(x.num)==0) x.num=0
x.num
})
#full.factorlist[[i]]$AlignID[which(ID.vct!=0)]
#ID.vct[which(ID.vct!=0)]
Experiment@Data@FactorList[[j]]$AlignID[ID.vct[which(ID.vct!=0)]] <- full.factorlist[[i]]$AlignID[which(ID.vct!=0)]
free.aligned.slots[[j]] <- which(full.factorlist[[i]]$AlignID[which(ID.vct!=0)]==0)
}
}
MaxALID <- max(unlist(lapply(Experiment@Data@FactorList, function(x) x$AlignID)))
Alid.counter <- MaxALID + 1
for(i in 1:length(free.aligned.slots))
{
Experiment@Data@FactorList[[i]]$AlignID[free.aligned.slots[[i]]] <- seq(Alid.counter, Alid.counter + (length(free.aligned.slots[[i]])-1) )
Alid.counter <- Alid.counter + length(free.aligned.slots[[i]])
}
cat("Constructing Factor List Table... (This may take a while...)\n")
Experiment@Results@Alignment <- create.factorlist.table(Experiment)
}
Experiment
}
)
#' @rdname identifyComp
#' @title Identification of compounds
#' @description Identification of compounds. Each empirical spectrum is compared against a ms library.
#' @usage identifyComp(Experiment, id.database = mslib,mz.range = NULL, n.putative = 3)
#' @param Experiment A 'MetaboSet' S4 object containing the experiment data previously created by newExp, deconvolved by deconvolveComp and optionally aligned by alignComp.
#' @param id.database The mass-spectra library to be compared with the empirical spectra. By default, the MassBank-[2] - Mass Bank of North America (MoNa) database are employed.
#' @param mz.range The same as in alignComp. If specified already in alignComp, then there is no need to especify it again. If not, it has to be specified.
#' @param n.putative The number of hits (compound candidate names) to be returned for each spectrum found.
#' @return The function returns an updated S4 'MetaboSet' class, where the GC-MS samples have been now aligned.
#' @references
#' [1] Xavier Domingo-Almenara, et al., eRah: A Computational Tool Integrating Spectral Deconvolution and Alignment with Quantification and Identification of Metabolites in GC-MS-Based Metabolomics. Analytical Chemistry (2016). DOI: 10.1021/acs.analchem.6b02927
#'
#' [2] MassBank: A public repository for sharing mass spectral data for life sciences, H. Horai, M. Arita, S. Kanaya, Y. Nihei, T. Ikeda, K. Suwa. Y. Ojima, K. Tanaka, S. Tanaka, K. Aoshima, Y. Oda, Y. Kakazu, M. Kusano, T. Tohge, F. Matsuda, Y. Sawada, M. Yokota Hirai, H. Nakanishi, K. Ikeda, N. Akimoto, T. Maoka, H. Takahashi, T. Ara, N. Sakurai, H. Suzuki, D. Shibata, S. Neumann, T. Iida, K. Tanaka, K. Funatsu, F. Matsuura, T. Soga, R. Taguchi, K. Saito and T. Nishioka, J. Mass Spectrom., 45 (2010) 703-714.
#' @author Xavier Domingo-Almenara. xavier.domingo@urv.cat
#' @seealso \code{\link{newExp}} \code{\link{alignComp}} \code{\link{setAlPar}} \code{\link{setDecPar}}
#' @export
setGeneric('identifyComp',function(Experiment, id.database=mslib, mz.range=NULL, n.putative=3){
standardGeneric('identifyComp')
})
#' @rdname identifyComp
setMethod('identifyComp',signature = 'MetaboSet',
function(Experiment, id.database=mslib, mz.range=NULL, n.putative=3){
#if(!(any(unlist(lapply(Experiment@Data@FactorList,function(x) {is.null(x$AlignID)} ))==FALSE))) stop("Factors must be aligned first")
if(is.null(Experiment@Results@Parameters@Alignment$mz.range) && is.null(mz.range)) stop("A mz.range has to be specified")
if(is.null(mz.range)) compare.only.mz <- min(Experiment@Results@Parameters@Alignment$mz.range):max(Experiment@Results@Parameters@Alignment$mz.range)
if(!is.null(mz.range)) compare.only.mz <- mz.range
id.par <- list(database.name = id.database@name, compare.only.mz = compare.only.mz, n.putative = n.putative)
Experiment@Results@Parameters@Identification <- id.par
avoid.processing.mz <- Experiment@Data@Parameters$avoid.processing.mz
maxMZ <- max(compare.only.mz)
Experiment@Results@Identification <- identify.factors(Experiment, maxMZ, compare.only.mz, avoid.processing.mz, id.database@database, n.putative)
Experiment
}
)
processSample <- function(Experiment, index, plotting, down.sample, virtual.scans.ps)
{
if(Experiment@MetaData@DataDirectory=="") {filename <- as.character(Experiment@MetaData@Instrumental$filename[index])
}else{filename <- paste(Experiment@MetaData@DataDirectory,"/",Experiment@MetaData@Instrumental$filename[index], sep="")}
sampleObject <- NULL
sampleObject <- load.file(filename)
# file.extension <- strsplit(as.character(Experiment@MetaData@Instrumental$filename[index]), split="\\.")[[1]]
# file.type <- file.extension[length(file.extension)]
# if(file.type=="cdf") sampleObject <- load.ncdf(filename)
# if(file.type=="mzXML" || file.type=="xml") sampleObject <- load.xml(filename)
# if(file.type=="MetaboSet")
# {
# load(filename)
# sampleObject <- new("RawDataParameters", data = sampleRD@data, min.mz = sampleRD@min.mz, max.mz = sampleRD@max.mz, start.time = sampleRD@start.time, mz.resolution = 1)
# }
Experiment@Data@Parameters$scans.per.second <- sampleObject@scans.per.second
sampleObject@avoid.processing.mz <- Experiment@Data@Parameters$avoid.processing.mz
sampleObject@min.peak.width <- Experiment@Data@Parameters$min.peak.width*Experiment@Data@Parameters$scans.per.second*60
sampleObject@min.peak.height <- Experiment@Data@Parameters$min.peak.height
sampleObject@noise.threshold <- Experiment@Data@Parameters$noise.threshold
#sampleObject@moving.window.length <- Experiment@Data@Parameters$moving.window.length*Experiment@Data@Parameters$scans.per.second*60
#sampleObject@moving.window.overlap <- Experiment@Data@Parameters$moving.window.overlap
sampleObject@compression.coef <- Experiment@Data@Parameters$compression.coef
#sampleObject@factor.minimum.sd <- Experiment@Data@Parameters$factor.minimum.sd
#sampleObject@filter.matrix <- get.filter.matrix(sampleObject)
sampleObject <- avoid.processing(sampleObject)
factor.list <- try(get.factor.list(sampleObject, analysis.window=Experiment@Data@Parameters$analysis.time, plotting, down.sample, virtual.scans.ps), silent=F)
if(inherits(factor.list,"try-error")) {factor.list <- as.data.frame(NULL); warning("Unable to extract factors from ", Experiment@MetaData@Instrumental$filename[index], ". Data may be corrupted.", sep="")}
factor.list
}
scansPerSecond <- function(Experiment){
if (Experiment@MetaData@DataDirectory=="") {
filename <- as.character(Experiment@MetaData@Instrumental$filename[1])
} else {
filename <- paste(Experiment@MetaData@DataDirectory,"/",Experiment@MetaData@Instrumental$filename[1], sep="")
}
sampleObject <- load.file(filename)
Experiment@Data@Parameters$scans.per.second <- sampleObject@scans.per.second
return(Experiment)
}
Try the erah package in your browser
Any scripts or data that you put into this service are public.
erah documentation built on June 22, 2024, 11:01 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions scansPerSecond processSample setAlPar setDecPar
Documented in setAlPar setDecPar
globalVariables("mslib")
#' @useDynLib erah
## Main Software functions:
#' @name setDecPar
#' @aliases setDecPar
#' @title Set Software Parameters
#' @description Sets Software Parameters for eRah.
#' @param min.peak.width Minimum compound peak width (in seconds). This is a critical parameter that conditions the efficiency of eRah. Typically, this should be the half of the mean compound width.
#' @param min.peak.height Minimum compound peak height
#' @param noise.threshold Data above this threshold will be considered as noise
#' @param avoid.processing.mz The masses that do not want to be considered for processing. Typically, in GC-MS those masses are 73,74,75,147,148 and 149, since they are they are ubiquitous mass fragments typically generated from compounds carrying a trimethylsilyl moiety.
#' @param compression.coef Data is compressed when using the orthogonal signal deconvolution (OSD) algorithm according to this value. A level 2 of compression is recomended.
#' @param analysis.time The chromatographic retention time window to process. If 0, all the chromatogram is processed.
#' @details See eRah vignette for more details. To open the vignette, execute the following code in R:
#' vignette("eRahManual", package="erah")
#' @references [1] Xavier Domingo-Almenara, et al., eRah: A Computational Tool Integrating Spectral Deconvolution and Alignment with Quantification and Identification of Metabolites in GC-MS-Based Metabolomics. Analytical Chemistry (2016). DOI: 10.1021/acs.analchem.6b02927
#' @author Xavier Domingo-Almenara. xavier.domingo@urv.cat
#' @seealso \code{\link{newExp}} \code{\link{deconvolveComp}} \code{\link{alignComp}} \code{\link{setAlPar}}
#' @examples \dontrun{
#' # The following will set eRah for analyzing the chromatograms
#' #from minutes 5 to 15, and withouth taking into account the masses
#' #35:69,73:75,147:149, widht a minimum peak widht of 0.7 seconds.
#' ex.dec.par <- setDecPar(min.peak.width = 0.7,
#' min.peak.height = 5000,
#' noise.threshold = 500,
#' avoid.processing.mz = c(35:69,73:75,147:149),
#' analysis.time = c(5,15))
#' }
#' @export
setDecPar <- function(min.peak.width,
min.peak.height=2500,
noise.threshold=500,
avoid.processing.mz=c(73:75,147:149),
compression.coef=2,
analysis.time=0){
softPar <- new("eRahSoftParameters",algorithm="eRah-OSD", min.peak.width = min.peak.width/60, min.peak.height = min.peak.height, noise.threshold = noise.threshold, avoid.processing.mz = avoid.processing.mz, compression.coef = compression.coef, analysis.time=analysis.time)
softPar
}
#' @name setAlPar
#' @aliases setAlPar
#' @title Set Alignment Parameters
#' @description Setting alignment parameters for eRah.
#' @usage setAlPar(min.spectra.cor, max.time.dist,mz.range = c(70:600))
#' @param min.spectra.cor Minimum spectral correlation value. From 0 (non similar) to 1 (very similar). This value sets how similar two or more compounds have be to be considered for alignment between them.
#' @param max.time.dist Maximum retention time distance. This value (in seconds) sets how far two or more compounds can be to be considered for alignment between them.
#' @param mz.range The range of masses that is considered when comparing spectra.
#' @references [1] Xavier Domingo-Almenara, et al., eRah: A Computational Tool Integrating Spectral Deconvolution and Alignment with Quantification and Identification of Metabolites in GC-MS-Based Metabolomics. Analytical Chemistry (2016). DOI: 10.1021/acs.analchem.6b02927
#' @author Xavier Domingo-Almenara. xavier.domingo@urv.cat
#' @seealso \code{\link{newExp}} \code{\link{setDecPar}} \code{\link{alignComp}}
#' @examples\dontrun{
#' # The following will set eRah for aligning compounds which are
#' # at least 90 (per cent) similar, and which peaks are at a
#' # maximum distance of 2 seconds. All the masses are considered when
#' # computing the spectral similarity.
#'
#' ex.al.par <- setAlPar(min.spectra.cor=0.90, max.time.dist=2,
#' mz.range=1:600)
#' }
#' @export
setAlPar <- function(min.spectra.cor, max.time.dist, mz.range=c(70:600))
{
alPar <- new("eRahAlParameters", algorithm="eRah", min.spectra.cor=min.spectra.cor, max.time.dist=max.time.dist/60, mz.range = mz.range, method="eRah")
alPar
}
#' @name newExp
#' @aliases newExp
#' @title New Experiment
#' @description Sets a new experiment for eRah
#' @param instrumental A data.frame containing the sample instrumental information.
#' @param phenotype (optional) A data.frame containing sample phenotype information.
#' @param info Experiment description
#' @details See eRah vignette for more details. To open the vignette, execute the following code in R:
#' vignette("eRahManual", package="erah")
#' @return \code{newExp} returns an S4 object of the class 'MetaboSet'.
#' @references [1] Xavier Domingo-Almenara, et al., eRah: A Computational Tool Integrating Spectral Deconvolution and Alignment with Quantification and Identification of Metabolites in GC-MS-Based Metabolomics. Analytical Chemistry (2016). DOI: 10.1021/acs.analchem.6b02927
#' @author Xavier Domingo-Almenara. xavier.domingo@urv.cat
#' @examples \dontrun{
#' library(gcspikelite)
#' data(targets)
#'
#' files <- list.files(system.file('data',package = 'gcspikelite'),full.names = TRUE)
#' files <- files[sapply(files,grepl,pattern = 'CDF')]
#'
#' instrumental <- createInstrumentalTable(files)
#' phenotype <- createPhenoTable(files,as.character(targets$Group[order(targets$FileName)]))
#'
#' ex <- newExp(instrumental = instrumental,
#' phenotype = phenotype, info = "DEMO Experiment")
#' }
#' @seealso \code{\link{createInstrumentalTable}} \code{\link{createPhenoTable}} \code{\link{setDecPar}} \code{\link{setAlPar}}
#' @importFrom utils read.csv
#' @export
newExp<- function (instrumental, phenotype = NULL, info = character())
{
if (is.null(phenotype)) {
if (is.data.frame(instrumental)){
instrumental <- instrumental
}
else if (is.character(instrumental) && file.exists(instrumental)){
instrumental<-read.csv(instrumental, sep=";")
}
else stop("Instrumental file format not valid.")
phenotype = as.data.frame(NULL)
factors.list <- lapply(1:nrow(instrumental), function(x) {
as.data.frame(NULL)
})
names(factors.list) <- as.vector(instrumental$sampleID)
ident.list <- as.data.frame(matrix(0, ncol = 7, dimnames = list(row = 0,
col = c("AlignID", "tmean", "Name", "MatchFactor", "CAS",
"Formula", "DB.Id"))))
uni.stats <- as.data.frame(matrix(0, ncol = 3, dimnames = list(row = 0,
col = c("Id", "FoldChangue", "pvalue"))))
multi.stats <- as.data.frame(matrix(0, ncol = 3, dimnames = list(row = 0,
col = c("Id", "CompoundsInvolved", "pvalue"))))
al.par <- list()
id.par <- list()
soft.par <- list()
stat.parameters <- new("MSResultsParameters", Alignment = al.par,
Identification = id.par)
statistics <- new("Statistics", Univariate = uni.stats,
Multivariate = multi.stats)
MS.Results <- new("Results", Parameters = stat.parameters,
Identification = ident.list, Statistics = statistics)
MS.Data <- new("Data", FeatureList = list(NULL), FactorList = factors.list,
Parameters = list(NULL))
MS.MetaData <- new("MetaData", Instrumental = instrumental,
Phenotype = phenotype, DataDirectory = "")
col.correct <- c("sampleID", "filename", "date", "time")
for (i in 1:length(col.correct)) if (length(apply(as.matrix(colnames(MS.MetaData@Instrumental)),
1, function(x) grep(col.correct[i], x))) == 0)
stop("Invalid instrumental file. The file must contain at least the following columns: ",
paste(col.correct, collapse = ", "))
if (!is.null(MS.MetaData@Phenotype))
col.correct <- c("sampleID", "class")
for (i in 1:length(col.correct)) if (length(apply(as.matrix(colnames(MS.MetaData@Phenotype)),
1, function(x) grep(col.correct[i], x))) == 0)
sample.container <- new("MetaboSet", Info = info, Data = MS.Data,
MetaData = MS.MetaData, Results = MS.Results)
sample.container
} else if (!is.null(phenotype)){
if (is.data.frame(instrumental)){
instrumental <- instrumental
}
else if (is.character(instrumental) && file.exists(instrumental)){
instrumental<-read.csv(instrumental, sep=";")
}
else stop("Instrumental file format not valid.")
if (is.data.frame(phenotype)){
phenotype <- phenotype
}
else if (is.character(phenotype) && file.exists(phenotype)){
phenotype <- read.csv(phenotype, sep=";")
}
else stop("Phenotype file format not valid.")
factors.list <- lapply(1:nrow(instrumental), function(x) {
as.data.frame(NULL)
})
names(factors.list) <- as.vector(instrumental$sampleID)
ident.list <- as.data.frame(matrix(0, ncol = 7, dimnames = list(row = 0,
col = c("AlignID", "tmean", "Name", "MatchFactor", "CAS",
"Formula", "DB.Id"))))
uni.stats <- as.data.frame(matrix(0, ncol = 3, dimnames = list(row = 0,
col = c("Id", "FoldChangue", "pvalue"))))
multi.stats <- as.data.frame(matrix(0, ncol = 3, dimnames = list(row = 0,
col = c("Id", "CompoundsInvolved", "pvalue"))))
al.par <- list()
id.par <- list()
soft.par <- list()
stat.parameters <- new("MSResultsParameters", Alignment = al.par,
Identification = id.par)
statistics <- new("Statistics", Univariate = uni.stats,
Multivariate = multi.stats)
MS.Results <- new("Results", Parameters = stat.parameters,
Identification = ident.list, Statistics = statistics)
MS.Data <- new("Data", FeatureList = list(NULL), FactorList = factors.list,
Parameters = list(NULL))
MS.MetaData <- new("MetaData", Instrumental = instrumental,
Phenotype = phenotype, DataDirectory = "")
col.correct <- c("sampleID", "filename", "date", "time")
for (i in 1:length(col.correct)) if (length(apply(as.matrix(colnames(MS.MetaData@Instrumental)),
1, function(x) grep(col.correct[i], x))) == 0)
stop("Invalid instrumental file. The file must contain at least the following columns: ",
paste(col.correct, collapse = ", "))
if (!is.null(MS.MetaData@Phenotype)) {
col.correct <- c("sampleID", "class")
for (i in 1:length(col.correct)) if (length(apply(as.matrix(colnames(MS.MetaData@Phenotype)),
1, function(x) grep(col.correct[i], x))) == 0)
stop("Invalid phenotype file. The file must contain at least the following columns: ",
paste(col.correct, collapse = ", "))
}
sample.container <- new("MetaboSet", Info = info, Data = MS.Data,
MetaData = MS.MetaData, Results = MS.Results)
sample.container
}
}
#' @rdname deconvolveComp
#' @title Deconvolution of compounds in samples
#' @description Deconvolution of GC-MS data
#' @param Experiment A 'MetaboSet' S4 object containing the experiment data previously created by newExp.
#' @param decParameters The software deconvolution parameters object previously created by setDecPar
#' @param samples.to.process Vector indicating which samples are to be processed.
#' @param down.sample If TRUE, chromatograms are down sampled to define one peak with 10 scan points (according to the minimum peak width). This is to process longer chromatograms with wider peak widths (more than 20 seconds peak width and small scans per second values). See details.
#' @param virtualScansPerSecond A virtual scans per second. If chromatograms are downsampled (for example, for a 1 mean peak width a 1 scans per second sampling frequency was used), eRah could not perform as expected. In these cases, the BEST solution is to re-acquire the samples. However, by selecting a different (virtual) scans per second frequency, eRah can upsample the data and process it more effectively.
#' @details See eRah vignette for more details. To open the vignette, execute the following code in R:
#' vignette("eRahManual", package="erah")
#'
#' eRah uses multivariate methods which run-time performance depend on the amount of data to be analyzed. When peaks are wider and the #' scans per second is also a small value, the number of points (scans) that define a peak might be too many, leading eRah to a poor run#'-time performance. To solve that, use down.sample=TRUE to allow eRah to define a peak with 10 seconds, and analyze the data more #' efficiently.
#' @return The function returns an updated S4 'MetaboSet' class, where the GC-MS samples have been now deconvolved.
#' @references [1] Xavier Domingo-Almenara, et al., eRah: A Computational Tool Integrating Spectral Deconvolution and Alignment with Quantification and Identification of Metabolites in GC-MS-Based Metabolomics. Analytical Chemistry (2016). DOI: 10.1021/acs.analchem.6b02927
#' @author Xavier Domingo-Almenara. xavier.domingo@urv.cat
#' @seealso \code{\link{newExp}} \code{\link{setAlPar}}
#' @examples \dontrun{
#' # Deconvolve data from a created experiment by \code{\link{newExp}}.
#' # ex <- newExp(instrumental="path")
#'
#' # The following will set eRah for analyzing the chromatograms
#' # from minutes 5 to 15, and withouth taking into account the masses
#' # 35:69,73:75,147:149, with a minimum peak width of 0.7 seconds.
#'
#' ex.dec.par <- setDecPar(min.peak.width=0.7, min.peak.height=5000,
#' noise.threshold=500, avoid.processing.mz=c(35:69,73:75,147:149),
#' analysis.time=c(5,15))
#'
#' # An now deconvolve the compounds in the samples:
#' # ex <- deconvolveComp(ex, decParameters=ex.dec.par)
#' }
#' @export
setGeneric('deconvolveComp',function(Experiment, decParameters, samples.to.process=NULL, down.sample=FALSE, virtualScansPerSecond=NULL){
standardGeneric('deconvolveComp')
})
#' @rdname deconvolveComp
#' @importFrom furrr future_map furrr_options
setMethod('deconvolveComp',signature = 'MetaboSet',
function(Experiment, decParameters, samples.to.process=NULL, down.sample=FALSE, virtualScansPerSecond=NULL){
plotting=FALSE
Number.of.Samples <- nrow(Experiment@MetaData@Instrumental)
if(is.null(samples.to.process)) samples.to.process <- 1:Number.of.Samples
stopifnot(samples.to.process>=1, max(samples.to.process)<=Number.of.Samples, length(samples.to.process)<=Number.of.Samples)
soft.par <- list(min.peak.width = decParameters@min.peak.width, min.peak.height = decParameters@min.peak.height, noise.threshold = decParameters@noise.threshold, avoid.processing.mz = decParameters@avoid.processing.mz, compression.coef = decParameters@compression.coef, analysis.time = decParameters@analysis.time)
Experiment@Data@Parameters <- soft.par
Experiment@Data@FactorList <- future_map(samples.to.process,~{
k <- which(samples.to.process == .x)
message(paste("\n Deconvolving compounds from",
as.character(Experiment@MetaData@Instrumental$filename[k]),
"... Processing",
k,"/",
length(samples.to.process),"\n"))
processSample(Experiment,
.x,
plotting,
down.sample,
virtualScansPerSecond)
},
.options = furrr_options(seed = TRUE))
names(Experiment@Data@FactorList) <- metaData(Experiment)$sampleID
Experiment <- scansPerSecond(Experiment)
message("Compounds deconvolved")
Experiment
}
)
#' @rdname alignComp
#' @title Alignment of compounds
#' @description Alignment of GC-MS deconvolved compounds
#' @usage alignComp(Experiment, alParameters, blocks.size=NULL)
#' @param Experiment A 'MetaboSet' S4 object containing the experiment data previously created by newExp and deconvolved by deconvolveComp.
#' @param alParameters The software alignment parameters object previously created by setAlPar
#' @param blocks.size For experiment of more than 1000 samples, and depending on the computer, alignment can be conducted by block segmentation. See details.
#' @details See eRah vignette for more details. To open the vignette, execute the following code in R:
#' vignette("eRahManual", package="erah")
#'
#' For experiments containing more than 100 (Windows) or 1000 (Mac or Linux) samples (numbers depending on the computer resoures and sample type). In those cases alignment can be conducted by block segmentation. For an experiment of e.g. 1000 samples, the block.size can be set to 100, so the alignment will perform as multiple (ten) 100-samples experiments, to later align them into a single experiment.
#'
#' This parameter is designed to solve the typical problem that appear when aligning under Windows operating system: "Error: cannot allocate vector of size XX Gb". Such a problem will not appear with Mac or Linux, but several hours of computation are expected when aligning a large number of samples. Using block segmentation provides a greatly improved run-time performance.
#' @return The function returns an updated S4 'MetaboSet' class, where the GC-MS samples have been now aligned.
#' @references [1] Xavier Domingo-Almenara, et al., eRah: A Computational Tool Integrating Spectral Deconvolution and Alignment with Quantification and Identification of Metabolites in GC-MS-Based Metabolomics. Analytical Chemistry (2016). DOI: 10.1021/acs.analchem.6b02927
#' @author Xavier Domingo-Almenara. xavier.domingo@urv.cat
#' @seealso \code{\link{newExp}} \code{\link{setDecPar}} \code{\link{deconvolveComp}}
#' @export
setGeneric('alignComp',function(Experiment, alParameters, blocks.size=NULL){
standardGeneric('alignComp')
})
#' @rdname alignComp
setMethod('alignComp',signature = 'MetaboSet',
function(Experiment, alParameters, blocks.size=NULL){
al.par <- list(alignment.algorithm=alParameters@algorithm, min.spectra.cor=alParameters@min.spectra.cor, max.time.dist=alParameters@max.time.dist, mz.range=alParameters@mz.range)
Experiment@Results@Parameters@Alignment <- al.par
min.spectra.cor <- Experiment@Results@Parameters@Alignment$min.spectra.cor
max.time.dist <- Experiment@Results@Parameters@Alignment$max.time.dist
mz.range <- Experiment@Results@Parameters@Alignment$mz.range
maxMZ <- max(mz.range)
# Experiment@Data@FactorList <- align.factors(Experiment@Data@FactorList, min.spectra.cor, max.time.dist, maxMZ, mz.range)
# Experiment@Results@Alignment <- create.factorlist.table(Experiment)
if(is.null(blocks.size))
{
Experiment@Data@FactorList <- align.factors(Experiment@Data@FactorList, min.spectra.cor, max.time.dist, maxMZ, mz.range)
Experiment@Results@Alignment <- create.factorlist.table(Experiment)
#return(Experiment)
}else{
#blocks.size <- 15
max.mz <- maxMZ
Itrt <- length(Experiment@Data@FactorList)/blocks.size
sequs <- trunc(seq(1, length(Experiment@Data@FactorList), length.out=Itrt))
sequs[1] <- 0
corresponding.list <- list()
block.list <- list()
#i <- 1
for(i in 1:(length(sequs)-1))
{
cat("Aligning block ", i, " of ", length(sequs)-1, "... \n", sep="")
ghost.object <- Experiment
ghost.object@Data@FactorList <- Experiment@Data@FactorList[(sequs[i]+1):sequs[(i+1)]]
factors.list <- ghost.object@Data@FactorList
ghost.object@Data@FactorList <- align.factors(factors.list, min.spectra.cor, max.time.dist, max.mz, mz.range)
ghost.factors.list <- create.factorlist.table(ghost.object)
block.list[[i]] <- data.frame(ID=ghost.factors.list$AlignID, RT=ghost.factors.list$tmean, Spectra=ghost.factors.list$Spectra)
corresponding.list <- c(corresponding.list,lapply(ghost.object@Data@FactorList, function(x) x$AlignID))
}
cat("Aligning factors across blocks... \n")
full.factorlist <- align.factors(block.list, min.spectra.cor, max.time.dist, max.mz, mz.range)
#MaxALID <- max(unlist(lapply(full.factorlist, function(x) x$AlignID)))
factors.list <- Experiment@Data@FactorList
if(!(any(unlist(lapply(factors.list,function(x) {is.null(x$AlignID)}))==FALSE)))
{
factors.list <- lapply(factors.list, function(x){
outp <- cbind(x,matrix(0,nrow=length(x$ID)))
colnames(outp)[ncol(outp)] <- "AlignID"
outp
})
}else{
factors.list <- lapply(factors.list, function(x){
x$AlignID <- rep(0,length(x$ID))
x
})
}
Experiment@Data@FactorList <- factors.list
free.aligned.slots <- list()
for(i in 1:length(full.factorlist))
{
for(j in (sequs[i]+1):sequs[(i+1)])
{
ID.vct <- sapply(full.factorlist[[i]]$ID, function(x) {x.num <- which(corresponding.list[[j]]==x)
if(length(x.num)==0) x.num=0
x.num
})
#full.factorlist[[i]]$AlignID[which(ID.vct!=0)]
#ID.vct[which(ID.vct!=0)]
Experiment@Data@FactorList[[j]]$AlignID[ID.vct[which(ID.vct!=0)]] <- full.factorlist[[i]]$AlignID[which(ID.vct!=0)]
free.aligned.slots[[j]] <- which(full.factorlist[[i]]$AlignID[which(ID.vct!=0)]==0)
}
}
MaxALID <- max(unlist(lapply(Experiment@Data@FactorList, function(x) x$AlignID)))
Alid.counter <- MaxALID + 1
for(i in 1:length(free.aligned.slots))
{
Experiment@Data@FactorList[[i]]$AlignID[free.aligned.slots[[i]]] <- seq(Alid.counter, Alid.counter + (length(free.aligned.slots[[i]])-1) )
Alid.counter <- Alid.counter + length(free.aligned.slots[[i]])
}
cat("Constructing Factor List Table... (This may take a while...)\n")
Experiment@Results@Alignment <- create.factorlist.table(Experiment)
}
Experiment
}
)
#' @rdname identifyComp
#' @title Identification of compounds
#' @description Identification of compounds. Each empirical spectrum is compared against a ms library.
#' @usage identifyComp(Experiment, id.database = mslib,mz.range = NULL, n.putative = 3)
#' @param Experiment A 'MetaboSet' S4 object containing the experiment data previously created by newExp, deconvolved by deconvolveComp and optionally aligned by alignComp.
#' @param id.database The mass-spectra library to be compared with the empirical spectra. By default, the MassBank-[2] - Mass Bank of North America (MoNa) database are employed.
#' @param mz.range The same as in alignComp. If specified already in alignComp, then there is no need to especify it again. If not, it has to be specified.
#' @param n.putative The number of hits (compound candidate names) to be returned for each spectrum found.
#' @return The function returns an updated S4 'MetaboSet' class, where the GC-MS samples have been now aligned.
#' @references
#' [1] Xavier Domingo-Almenara, et al., eRah: A Computational Tool Integrating Spectral Deconvolution and Alignment with Quantification and Identification of Metabolites in GC-MS-Based Metabolomics. Analytical Chemistry (2016). DOI: 10.1021/acs.analchem.6b02927
#'
#' [2] MassBank: A public repository for sharing mass spectral data for life sciences, H. Horai, M. Arita, S. Kanaya, Y. Nihei, T. Ikeda, K. Suwa. Y. Ojima, K. Tanaka, S. Tanaka, K. Aoshima, Y. Oda, Y. Kakazu, M. Kusano, T. Tohge, F. Matsuda, Y. Sawada, M. Yokota Hirai, H. Nakanishi, K. Ikeda, N. Akimoto, T. Maoka, H. Takahashi, T. Ara, N. Sakurai, H. Suzuki, D. Shibata, S. Neumann, T. Iida, K. Tanaka, K. Funatsu, F. Matsuura, T. Soga, R. Taguchi, K. Saito and T. Nishioka, J. Mass Spectrom., 45 (2010) 703-714.
#' @author Xavier Domingo-Almenara. xavier.domingo@urv.cat
#' @seealso \code{\link{newExp}} \code{\link{alignComp}} \code{\link{setAlPar}} \code{\link{setDecPar}}
#' @export
setGeneric('identifyComp',function(Experiment, id.database=mslib, mz.range=NULL, n.putative=3){
standardGeneric('identifyComp')
})
#' @rdname identifyComp
setMethod('identifyComp',signature = 'MetaboSet',
function(Experiment, id.database=mslib, mz.range=NULL, n.putative=3){
#if(!(any(unlist(lapply(Experiment@Data@FactorList,function(x) {is.null(x$AlignID)} ))==FALSE))) stop("Factors must be aligned first")
if(is.null(Experiment@Results@Parameters@Alignment$mz.range) && is.null(mz.range)) stop("A mz.range has to be specified")
if(is.null(mz.range)) compare.only.mz <- min(Experiment@Results@Parameters@Alignment$mz.range):max(Experiment@Results@Parameters@Alignment$mz.range)
if(!is.null(mz.range)) compare.only.mz <- mz.range
id.par <- list(database.name = id.database@name, compare.only.mz = compare.only.mz, n.putative = n.putative)
Experiment@Results@Parameters@Identification <- id.par
avoid.processing.mz <- Experiment@Data@Parameters$avoid.processing.mz
maxMZ <- max(compare.only.mz)
Experiment@Results@Identification <- identify.factors(Experiment, maxMZ, compare.only.mz, avoid.processing.mz, id.database@database, n.putative)
Experiment
}
)
processSample <- function(Experiment, index, plotting, down.sample, virtual.scans.ps)
{
if(Experiment@MetaData@DataDirectory=="") {filename <- as.character(Experiment@MetaData@Instrumental$filename[index])
}else{filename <- paste(Experiment@MetaData@DataDirectory,"/",Experiment@MetaData@Instrumental$filename[index], sep="")}
sampleObject <- NULL
sampleObject <- load.file(filename)
# file.extension <- strsplit(as.character(Experiment@MetaData@Instrumental$filename[index]), split="\\.")[[1]]
# file.type <- file.extension[length(file.extension)]
# if(file.type=="cdf") sampleObject <- load.ncdf(filename)
# if(file.type=="mzXML" || file.type=="xml") sampleObject <- load.xml(filename)
# if(file.type=="MetaboSet")
# {
# load(filename)
# sampleObject <- new("RawDataParameters", data = sampleRD@data, min.mz = sampleRD@min.mz, max.mz = sampleRD@max.mz, start.time = sampleRD@start.time, mz.resolution = 1)
# }
Experiment@Data@Parameters$scans.per.second <- sampleObject@scans.per.second
sampleObject@avoid.processing.mz <- Experiment@Data@Parameters$avoid.processing.mz
sampleObject@min.peak.width <- Experiment@Data@Parameters$min.peak.width*Experiment@Data@Parameters$scans.per.second*60
sampleObject@min.peak.height <- Experiment@Data@Parameters$min.peak.height
sampleObject@noise.threshold <- Experiment@Data@Parameters$noise.threshold
#sampleObject@moving.window.length <- Experiment@Data@Parameters$moving.window.length*Experiment@Data@Parameters$scans.per.second*60
#sampleObject@moving.window.overlap <- Experiment@Data@Parameters$moving.window.overlap
sampleObject@compression.coef <- Experiment@Data@Parameters$compression.coef
#sampleObject@factor.minimum.sd <- Experiment@Data@Parameters$factor.minimum.sd
#sampleObject@filter.matrix <- get.filter.matrix(sampleObject)
sampleObject <- avoid.processing(sampleObject)
factor.list <- try(get.factor.list(sampleObject, analysis.window=Experiment@Data@Parameters$analysis.time, plotting, down.sample, virtual.scans.ps), silent=F)
if(inherits(factor.list,"try-error")) {factor.list <- as.data.frame(NULL); warning("Unable to extract factors from ", Experiment@MetaData@Instrumental$filename[index], ". Data may be corrupted.", sep="")}
factor.list
}
scansPerSecond <- function(Experiment){
if (Experiment@MetaData@DataDirectory=="") {
filename <- as.character(Experiment@MetaData@Instrumental$filename[1])
} else {
filename <- paste(Experiment@MetaData@DataDirectory,"/",Experiment@MetaData@Instrumental$filename[1], sep="")
}
sampleObject <- load.file(filename)
Experiment@Data@Parameters$scans.per.second <- sampleObject@scans.per.second
return(Experiment)
}
Try the erah package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.