Documentation and Public Functions Update

xdomingoal · Sep 27, 2023 · 3636ce6 · 3636ce6
1 parent cf015e7
commit 3636ce6
Show file tree

Hide file tree

Showing 11 changed files with 263 additions and 141 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: erah
 Type: Package
 Version: 2.0.0
-Author: Xavier Domingo-Almenara [aut, cre, cph], Jasen P. Finch [ctb], Sara Samino [aut], Maria Vinaixa [aut], Alexandre Perera [aut, ths], Jesus Brezmes [aut, ths], Oscar Yanes [aut, ths]
+Author: Xavier Domingo-Almenara [aut, cre, cph], Jasen P. Finch [ctb], Adria Olomi [ctb], Sara Samino [aut], Maria Vinaixa [aut], Alexandre Perera [aut, ths], Jesus Brezmes [aut, ths], Oscar Yanes [aut, ths]
 Title: Automated Spectral Deconvolution, Alignment, and Metabolite Identification in GC/MS-Based Untargeted Metabolomics
 Depends: R (>= 3.4)
 Imports: osd,  
@@ -32,5 +32,5 @@ biocViews: MassSpectrometry, Metabolomics
 LazyData: yes
 NeedsCompilation: yes
 VignetteBuilder: knitr
-RoxygenNote: 7.1.1
+RoxygenNote: 7.2.3
 Encoding: UTF-8
diff --git a/NAMESPACE b/NAMESPACE
@@ -1,6 +1,7 @@
 # Generated by roxygen2: do not edit by hand
 
 export(compInfo)
+export(computeRIerror)
 export(createInstrumentalTable)
 export(createPhenoTable)
 export(createdt)
@@ -15,6 +16,7 @@ export(setAlPar)
 export(setDecPar)
 export(showRTRICurve)
 exportClasses(MetaboSet)
+exportClasses(RawDataParameters)
 exportClasses(eRah_DB)
 exportMethods(alignComp)
 exportMethods(alignList)

diff --git a/R/RI.R b/R/RI.R
@@ -1,4 +1,22 @@
-#' @importFrom stats smooth.spline predict
+#' @name computeRIerror
+#' @aliases computeRIerror
+#' @title computeRIerror
+#' @description This function uses RI of mslib database and RT of the identified compounds to discrimine proper compound identification.
+#' @param Experiment S4 object with experiment Data, Metadata and Results. Results of experiment are used to extract RT and Compound DB Id.
+#' @param id.database Name of the preloaded database, in this case the regular db used by erah mslib
+#' @param reference.list List with the compounds and their attributes (AlignId...)
+#' @param ri.error.type Specify wether absolute or relative RI error is to be computed.
+#' @param plot.results Shows the RI/RT graphic (True by default)
+#' @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{showRTRICurve}} 
+#' @examples \dontrun{
+#' ex <- computeRIerror(ex, mslib, reference.list=list(AlignID = c(45,67,92,120)), ri.error.type = "relative")
+#' }
+#' @export
+#' @importFrom stats smooth.spline predict 
 #' @importFrom graphics points
 
 computeRIerror <- function(Experiment, id.database=mslib, reference.list, ri.error.type=c('relative','absolute'), plot.results=TRUE){
@@ -86,94 +104,97 @@ computeRIerror <- function(Experiment, id.database=mslib, reference.list, ri.err
 
 #' @name showRTRICurve
 #' @aliases showRTRICurve
-#' @title RTRI CURVE
+#' @title Show RT-RI curve
 #' @description This function uses RI of mslib database and RT of the identified compounds to discrimine proper compound identification.
 #' @param Experiment S4 object with experiment Data, Metadata and Results. Results of experiment are used to extract RT and Compound DB Id.
 #' @param reference.list List with the compounds and their attributes (AlignId...)
 #' @param nAnchors The desired equivalent number of degrees of freedom for the smooth.spline function
 #' @param ri.thrs Retention Index treshold given by the user to discrimine bewteen identification results
-#' @param id.database Name of the preloaded database, in this case the regular db used by erah mslib
+#' @param id.database Name of the preloaded database (mslib by default, the regular db used by erah)
 #' @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{ComputeRIerror}} 
 #' @examples \dontrun{
 #' The following set erah to determine which indetified compounds are in RI treshold
-#' RTRICurve <- showRTRICurve(ex, list, nAnchors=4, ri.thrs='1R', id.database = mslib)
+#' RTRICurve <- showRTRICurve(ex, list, nAnchors=4, ri.thrs='1R')
 #' }
 #' @export
 
 showRTRICurve <- function(Experiment, reference.list, nAnchors=4, ri.thrs='1R', id.database = mslib){
-	 if (class(reference.list) != "list") 
-        stop("The parameter reference.list must be a list")
-    if (is.null(reference.list)) 
-        stop("A reference list must be provided")
-    if (!is.null(reference.list$AlignID) & !(is.null(reference.list$RT) | 
-        is.null(reference.list$RI))) 
-        stop("reference.list must contain a) user-defined RT and RI  values, or b) the AlignID of the compounds to be used as a reference")
-    if (is.null(reference.list$AlignID) & (is.null(reference.list$RT) | 
-        is.null(reference.list$RI))) 
-        stop("reference.list must contain a) user-defined RT and RI  values, or b) the AlignID of the compounds to be used as a reference")
-    if (!is.null(reference.list$RT)) 
-        if (length(reference.list$RT) != length(reference.list$RI)) 
-            stop("Both RI and RT vectors in reference list must have the same length! (You provided a different number of RT and RI values, they must have the same length)")
-    if (is.null(id.database)) 
-        stop("A database is needed for spectra comparison. Select a database or set 'compare' parameter to 'False'")
-    if (nrow(Experiment@Results@Identification) == 0) 
-        stop("Factors must be identified first")
-    if (Experiment@Results@Parameters@Identification$database.name != 
-        id.database@name) {
-        error.msg <- paste("This experiment was not processed with the database selected. Please use ", 
-            Experiment@Results@Parameters@Identification$database.name, 
-            sep = "")
-        stop(error.msg)
+  if (class(reference.list) != "list") 
+    stop("The parameter reference.list must be a list")
+  if (is.null(reference.list)) 
+    stop("A reference list must be provided")
+  if (!is.null(reference.list$AlignID) & !(is.null(reference.list$RT) | 
+                                           is.null(reference.list$RI))) 
+    stop("reference.list must contain a) user-defined RT and RI  values, or b) the AlignID of the compounds to be used as a reference")
+  if (is.null(reference.list$AlignID) & (is.null(reference.list$RT) | 
+                                         is.null(reference.list$RI))) 
+    stop("reference.list must contain a) user-defined RT and RI  values, or b) the AlignID of the compounds to be used as a reference")
+  if (!is.null(reference.list$RT)) 
+    if (length(reference.list$RT) != length(reference.list$RI)) 
+      stop("Both RI and RT vectors in reference list must have the same length! (You provided a different number of RT and RI values, they must have the same length)")
+  if (is.null(id.database)) 
+    stop("A database is needed for spectra comparison. Select a database or set 'compare' parameter to 'False'")
+  if (nrow(Experiment@Results@Identification) == 0) 
+    stop("Factors must be identified first")
+  if (Experiment@Results@Parameters@Identification$database.name != 
+      id.database@name) {
+    error.msg <- paste("This experiment was not processed with the database selected. Please use ", 
+                       Experiment@Results@Parameters@Identification$database.name, 
+                       sep = "")
+    stop(error.msg)
+  }
+
+  idlist <- idList(object = Experiment, id.database = id.database)
+
+  if (!is.null(reference.list$AlignID)) {
+    refInd <- which(Experiment@Results@Identification$AlignID %in% 
+                      reference.list$AlignID)
+    riVect.std <- sapply(as.numeric(as.vector(Experiment@Results@Identification$DB.Id.1[refInd])), 
+                         function(i) id.database@database[[i]]$RI.VAR5.ALK)
+    rtVect.std <- as.numeric(as.vector(Experiment@Results@Identification$tmean[refInd]))
+  }
+  else {
+    riVect.std <- reference.list$RI
+    rtVect.std <- reference.list$RT
+  }   
+
+  relativeThr <- NULL
+  if(length(grep('r', tolower(ri.thrs)))==1) {
+    relativeThr <- TRUE
+    riThrs <- as.numeric(as.vector(gsub('r','', tolower(ri.thrs))))
+  }	
+  if(length(grep('a', tolower(ri.thrs)))==1) {
+    relativeThr <- FALSE
+    riThrs <- as.numeric(as.vector(gsub('a','', tolower(ri.thrs))))
+  }	
+  if(is.null(relativeThr)) stop('Incorrect "ri.thrs" value, please check the documentation for assinging "ri.thrs" values.')
+
+  sSp <- smooth.spline(rtVect.std , riVect.std,df=nAnchors)
+  riError <- abs(sSp$y - riVect.std)
+  if(relativeThr) riError <- 100*(riError/sSp$y)
+
+  plot(rtVect.std, riVect.std, type='n', main='RT/RI curve', xlab='RT (min)', ylab='RI')
+  points(rtVect.std[which(riError<riThrs)], riVect.std[which(riError<riThrs)], pch=19)
+  if(length(which(riError>riThrs))!=0) points(rtVect.std[which(riError>riThrs)], riVect.std[which(riError>riThrs)], pch=19, col='red2')
+  lines(sSp)
+
+  if(length(which(riError>riThrs))!=0){
+
+    if (!is.null(reference.list$AlignID)) {
+      cat('The points outside the RI threshold are: ')
+      cat(paste0(Experiment@Results@Identification$AlignID[refInd][which(riError>riThrs)], collapse=', '))
+    }else{
+      cat('The points outside the RI threshold are: \n')
+      #cbind(reference.list$RI[which(riError>riThrs)],reference.list$RT[which(riError>riThrs)])
     }
+  }else{
+    cat('No points outside the selected RI threshold')
+  }	
+}
+
 
-	idlist <- idList(object = Experiment, id.database = id.database)
-
-	if (!is.null(reference.list$AlignID)) {
-        refInd <- which(Experiment@Results@Identification$AlignID %in% 
-            reference.list$AlignID)
-        riVect.std <- sapply(as.numeric(as.vector(Experiment@Results@Identification$DB.Id.1[refInd])), 
-            function(i) id.database@database[[i]]$RI.VAR5.ALK)
-        rtVect.std <- as.numeric(as.vector(Experiment@Results@Identification$tmean[refInd]))
-    }
-    else {
-        riVect.std <- reference.list$RI
-        rtVect.std <- reference.list$RT
-    }   
-
-	relativeThr <- NULL
-	if(length(grep('r', tolower(ri.thrs)))==1) {
-		relativeThr <- TRUE
-		riThrs <- as.numeric(as.vector(gsub('r','', tolower(ri.thrs))))
-	}	
-	if(length(grep('a', tolower(ri.thrs)))==1) {
-		relativeThr <- FALSE
-		riThrs <- as.numeric(as.vector(gsub('a','', tolower(ri.thrs))))
-	}	
-	if(is.null(relativeThr)) stop('Incorrect "ri.thrs" value, please check the documentation for assinging "ri.thrs" values.')
-
-	sSp <- smooth.spline(rtVect.std , riVect.std,df=nAnchors)
-	riError <- abs(sSp$y - riVect.std)
-	if(relativeThr) riError <- 100*(riError/sSp$y)
 
-	plot(rtVect.std, riVect.std, type='n', main='RT/RI curve', xlab='RT (min)', ylab='RI')
-	points(rtVect.std[which(riError<riThrs)], riVect.std[which(riError<riThrs)], pch=19)
-	if(length(which(riError>riThrs))!=0) points(rtVect.std[which(riError>riThrs)], riVect.std[which(riError>riThrs)], pch=19, col='red2')
-	lines(sSp)
-
-	if(length(which(riError>riThrs))!=0){
-
-		if (!is.null(reference.list$AlignID)) {
-			cat('The points outside the RI threshold are: ')
-			cat(paste0(Experiment@Results@Identification$AlignID[refInd][which(riError>riThrs)], collapse=', '))
-      }else{
-      	cat('The points outside the RI threshold are: \n')
-      	#cbind(reference.list$RI[which(riError>riThrs)],reference.list$RT[which(riError>riThrs)])
-      }
-	}else{
-		cat('No points outside the selected RI threshold')
-	}	
-}
diff --git a/R/allClasses.R b/R/allClasses.R
@@ -89,24 +89,24 @@ setClass(Class = "eRah_DB",
                                          database = "list")
 )
 
-#'  @name RawDataParameters-class
-#'  @docType class
-#'  @aliases RawDataParameters-class
-#'  @title Class \code{"RawDataParameters"}
-#'  @description The RawDataParameters class contains the slots for storing and accessing into a MS sample, and the essential parameters for performing its processing (deconvolution).
-#'  @slot data The data matrix of the sample to be processed
-#'  @slot min.mz The minimum adquired mz number
-#'  @slot max.mz The maximum adquired mz number
-#'  @slot start.time Starting time of adquisition
-#'  @slot mz.resolution Mz resolution
-#'  @slot scans.per.second Scans per second
-#'  @slot avoid.processing.mz Which mz do not have to be processed
-#'  @slot min.peak.width Minimum peak width (stored in scans)
-#'  @slot min.peak.height Minimum peak height
-#'  @slot noise.threshold The noise threshold
-#'  @slot compression.coef Compression coefficient (parameter for Orthogonal Signal Deconvolution)
-#'  @author Xavier Domingo-Almenara.
-#'  @export
+#' @name RawDataParameters-class
+#' @docType class
+#' @aliases RawDataParameters-class
+#' @title Class \code{"RawDataParameters"}
+#' @description The RawDataParameters class contains the slots for storing and accessing into a MS sample, and the essential parameters for performing its processing (deconvolution).
+#' @slot data The data matrix of the sample to be processed
+#' @slot min.mz The minimum adquired mz number
+#' @slot max.mz The maximum adquired mz number
+#' @slot start.time Starting time of adquisition
+#' @slot mz.resolution Mz resolution
+#' @slot scans.per.second Scans per second
+#' @slot avoid.processing.mz Which mz do not have to be processed
+#' @slot min.peak.width Minimum peak width (stored in scans)
+#' @slot min.peak.height Minimum peak height
+#' @slot noise.threshold The noise threshold
+#' @slot compression.coef Compression coefficient (parameter for Orthogonal Signal Deconvolution)
+#' @author Xavier Domingo-Almenara.
+#' @export
 
 setClass(Class = "RawDataParameters", 
          representation = representation(data = "matrix", 

diff --git a/R/allGenerics.R b/R/allGenerics.R
@@ -55,7 +55,7 @@ setGeneric('plotSpectra',function(Experiment, AlignId, n.putative=1, compare=T,
 
 #' @rdname plotProfile
 
-setGeneric('plotProfile',function(Experiment,AlignId, per.class=T, xlim=NULL){
+setGeneric('plotProfile',function(Experiment, AlignId, per.class=T, xlim=NULL, cols=NULL){
   standardGeneric('plotProfile')
 })
 

diff --git a/R/eRah_modified.tar.gz b/R/eRah_modified.tar.gz
diff --git a/R/plotting.R b/R/plotting.R
@@ -112,14 +112,16 @@ setMethod(plotSpectra,signature = 'MetaboSet',
 #' @usage plotProfile(Experiment,AlignId, per.class = T, xlim = NULL)
 #' @param Experiment A 'MetaboSet' S4 object containing the experiment after being deconolved, aligned and (optionally) identified.
 #' @param AlignId the Id identificator for the compound to be shown.
-#' @param per.class logical. if TRUE the profiles are shown one color per class, if FALSE one color per sample.
+#' @param per.class logical. if TRUE (by default) the profiles are shown one color per class, if FALSE one color per sample.
 #' @param xlim x axsis (retention time) limits (see \code{\link{plot.default}}).
+#' @param cols vector of colors. Colors are used cyclically.
 #' @author Xavier Domingo-Almenara. xavier.domingo@urv.cat
 #' @seealso \code{\link{plotSpectra}} \code{\link{plotAlign}}
 #' @export
 
 setMethod('plotProfile',signature = 'MetaboSet',
-          function(Experiment,AlignId, per.class=T, xlim=NULL){	
+          function(Experiment, AlignId, per.class=T, xlim=NULL, cols=NULL)
+          {	
             if(!(any(unlist(lapply(Experiment@Data@FactorList,function(x) {is.null(x$AlignID)} ))==FALSE))) stop("Factors must be aligned first")
 
             if(nrow(Experiment@MetaData@Phenotype)==0) 
@@ -136,7 +138,7 @@ setMethod('plotProfile',signature = 'MetaboSet',
             N.groups <- max(unique(unlist(alignId)))	
             N.samples <- length(Experiment@Data@FactorList)
 
-            samples.name <- metaData(Experiment)$sampleID		
+            samples.name <- names(Experiment@Data@FactorList)		
             for(i in 1:N.samples) samples.name[i] <- strsplit(as.character(samples.name[i]), split="\\.")[[1]][1]
 
             profile.list <- lapply(Experiment@Data@FactorList,function(x) {
@@ -145,7 +147,7 @@ setMethod('plotProfile',signature = 'MetaboSet',
               int <- NA
               if(length(outp)!=0)
               {
-                output <- sparse.to.vector(outp)
+                output <- erah:::sparse.to.vector(outp)
                 time <- output$time
                 int <- output$int*as.numeric(as.character(x[which(x$AlignID==AlignId),"Peak Height"]))
               }
@@ -178,8 +180,7 @@ setMethod('plotProfile',signature = 'MetaboSet',
             compound.name <- as.character(Experiment@Results@Identification[which(Experiment@Results@Identification$AlignID==AlignId),"Name"])
 
             if(per.class==F)
-            {
-
+            {	
               matplot(profile.time, profile.int, type="l", lty=1, col=(1:length(samples.name)), main=paste("Profile Comparison \n",compound.name), xlab="time (min)", ylab="Intensity", xlim=xlim)
               par(font=2)
               legend("topright",legend=samples.name, pch=19, col=(1:length(samples.name)), title="Samples")
@@ -188,17 +189,25 @@ setMethod('plotProfile',signature = 'MetaboSet',
               pn <- Experiment@MetaData@Phenotype
               indx <- apply(as.matrix(samples.name),1,function(x) which(pn[,"sampleID"]==x))		
               class.names <- pn[indx,"class"]
+              ssCC <- split(1:length(indx), pn[indx,"class"])
+
+              class.colors <- unlist(sapply(1:length(ssCC), function(x) rep(x,length(ssCC[[x]]))))
+              if(!is.null(cols)){
+                if(!(length(cols)==length(ssCC))) stop('The length of the specified colors does not match the number of classes used in this experiment (',length(ssCC),')')
+                class.colors <- sapply(class.colors, function(x) cols[x])
+              }	
 
-              samples.class.type <- levels(factor(pn$class))
+              samples.class.type <- unique(pn$class)
 
-              matplot(profile.time,profile.int, type="l", lty=1, col=1:length(samples.class.type), main=paste("Profile Comparison \n",compound.name), xlab="time (min)", ylab="Intensity", xlim=xlim)
+              matplot(profile.time,profile.int, type="l", lty=1, col=class.colors, main=paste("Profile Comparison \n",compound.name), xlab="time (min)", ylab="Intensity", xlim=xlim, lwd=2)
               par(font=2)
-              legend("topright",legend=samples.class.type, pch=19, col=1:length(samples.class.type), title="Classes")
+              legend("topright",legend=samples.class.type, pch=19, col=unique(class.colors), title="Classes")
               par(font=1)
 
             }
 
           }
+
 )
 
 #' @rdname plotAlign