diff --git a/tools/metams/metaMS_plot/metaMS_plot.r b/tools/metams/metaMS_plot/metaMS_plot.r index 11f580093..21456fb78 100644 --- a/tools/metams/metaMS_plot/metaMS_plot.r +++ b/tools/metams/metaMS_plot/metaMS_plot.r @@ -1,39 +1,39 @@ #!/usr/bin/env Rscript # metaMS_plot.r version="1.0.0" -#created by Yann GUITTON and updated by Julien SAINT-VANNE +# created by Yann GUITTON and updated by Julien SAINT-VANNE # ----- LOG FILE ----- -#log_file <- file("log.txt", open = "wt") -#sink(log_file) -#sink(log_file, type = "output") +# log_file <- file("log.txt", open = "wt") +# sink(log_file) +# sink(log_file, type = "output") # ----- PACKAGE ----- cat("\tSESSION INFO\n") -#Import the different functions +# Import the different functions source_local <- function(fname) { - argv <- commandArgs(trailingOnly = FALSE) - base_dir <- dirname(substring(argv[grep("--file=", argv)], 8)) - source(paste(base_dir, fname, sep="/")) + argv <- commandArgs(trailingOnly = FALSE) + base_dir <- dirname(substring(argv[grep("--file=", argv)], 8)) + source(paste(base_dir, fname, sep = "/")) } source_local("lib_metams.r") -pkgs <- c("metaMS","batch") #"batch" necessary for parseCommandArgs function +pkgs <- c("metaMS", "batch") # "batch" necessary for parseCommandArgs function loadAndDisplayPackages(pkgs) cat("\n") modNamC <- "plot_metaMS" ## module name -cat("\nStart of the '", modNamC, "' Galaxy module call: ", format(Sys.time(), "%a %d %b %Y %X"), "\n\n", sep="") +cat("\nStart of the '", modNamC, "' Galaxy module call: ", format(Sys.time(), "%a %d %b %Y %X"), "\n\n", sep = "") # ----- ARGUMENTS ----- cat("\tARGUMENTS INFO\n\n") -args = parseCommandArgs(evaluate=FALSE) #interpretation of arguments given in command line as an R list of objects -#write.table(as.matrix(args), col.names=F, quote=F, sep='\t\t') +args <- parseCommandArgs(evaluate = FALSE) # interpretation of arguments given in command line as an R list of objects +# write.table(as.matrix(args), col.names=F, quote=F, sep='\t\t') print(cbind(value = unlist(args))) @@ -44,31 +44,31 @@ cat("\n\n\tARGUMENTS PROCESSING INFO\n\n") load(args$metaMS) if (!exists("resGC")) stop("\n\nERROR: The RData doesn't contain any object called 'resGC' which is provided by the tool: new_metaMS.runGC") -if(args$selecteic) { - #Unknown EIC parameter - if (args$unkn[1] != "NULL") { - #When unkn = 0 user want to process all unknowns - if(args$unkn[1] == 0) { - args$unkn <- c(1:nrow(resGC$PeakTable)) - print("User want to process on all unknown(s) found in metaMS process") - } - #TODO find the biggest number of unkn ask by user cause it can write "1,15,9,8" with a max of 11 unkn. With this code it finds the 8 and it will pass - #Verify that there is not more user's unkn than metaMS unkn (find in resGC$PeakTable) - cat("Number of unknown after metaMS process :",nrow(resGC$PeakTable),"\n") - cat("Number of the last unknown ask by user :",args$unkn[length(args$unkn)],"\n") - cat("Number of unknown ask by user :",length(args$unkn),"\n") - if(args$unkn[length(args$unkn)] <= nrow(resGC$PeakTable)) { - unknarg <- args$unkn - } else { - error_message="Too much unkn compare metaMS results" - print(error_message) - stop(error_message) - } +if (args$selecteic) { + # Unknown EIC parameter + if (args$unkn[1] != "NULL") { + # When unkn = 0 user want to process all unknowns + if (args$unkn[1] == 0) { + args$unkn <- c(1:nrow(resGC$PeakTable)) + print("User want to process on all unknown(s) found in metaMS process") + } + # TODO find the biggest number of unkn ask by user cause it can write "1,15,9,8" with a max of 11 unkn. With this code it finds the 8 and it will pass + # Verify that there is not more user's unkn than metaMS unkn (find in resGC$PeakTable) + cat("Number of unknown after metaMS process :", nrow(resGC$PeakTable), "\n") + cat("Number of the last unknown ask by user :", args$unkn[length(args$unkn)], "\n") + cat("Number of unknown ask by user :", length(args$unkn), "\n") + if (args$unkn[length(args$unkn)] <= nrow(resGC$PeakTable)) { + unknarg <- args$unkn } else { - error_message <- "No EIC selected !" - print(error_message) - stop(error_message) + error_message <- "Too much unkn compare metaMS results" + print(error_message) + stop(error_message) } + } else { + error_message <- "No EIC selected !" + print(error_message) + stop(error_message) + } } cat("\n\n") @@ -92,67 +92,67 @@ cat("\n\n\tMAIN PROCESSING INFO\n") cat("\t\tCOMPUTE\n") -#Use getTIC2s and getBPC2s because getTICs and getBPCs can exists due to transfert of function in Rdata - -if(!is.null(singlefile)) { - files <- paste("./",names(singlefile),sep="") - #WARNING if user has CDF files (not yet good for plotting) - if(MSnbase:::isCdfFile(files)){ - warning_message <- "You have CDF files, for the moment you can't obtain plot after runGC! A new update will follow with the good correction\n" - warning(warning_message) - cat(paste("\n","/!\\Warning/!\\",warning_message,sep="\n")) +# Use getTIC2s and getBPC2s because getTICs and getBPCs can exists due to transfert of function in Rdata + +if (!is.null(singlefile)) { + files <- paste("./", names(singlefile), sep = "") + # WARNING if user has CDF files (not yet good for plotting) + if (MSnbase:::isCdfFile(files)) { + warning_message <- "You have CDF files, for the moment you can't obtain plot after runGC! A new update will follow with the good correction\n" + warning(warning_message) + cat(paste("\n", "/!\\Warning/!\\", warning_message, sep = "\n")) + } + if (!is.null(files)) { + if (args$selectbpc) { + cat("\n\tProcessing BPC(s) from XCMS files...\n") + c <- getBPC2s(files = files, xset = xset, rt = "raw", pdfname = "BPCs_raw.pdf") + cat("BPC(s) created...\n") } - if(!is.null(files)){ - if(args$selectbpc){ - cat("\n\tProcessing BPC(s) from XCMS files...\n") - c <- getBPC2s(files = files, xset = xset, rt="raw", pdfname="BPCs_raw.pdf") - cat("BPC(s) created...\n") - } - if(args$selecttic){ - cat("\n\tProcessing TIC(s) from XCMS files...\n") - b <- getTIC2s(files = files, xset = xset, rt="raw", pdfname="TICs_raw.pdf") - cat("TIC(s) created...\n") - } - if(args$selecteic){ - cat("\n\tProcessing EIC(s) from XCMS files...\n") - cat(length(unknarg),"unknown(s) will be process !\n") - plotUnknowns(resGC=resGC, unkn=unknarg, DB=DBgc, fileFrom="singlefile") - cat("EIC(s) created...\n") - #Merge pdf - system(paste('gs -o GCMS_EIC.pdf -sDEVICE=pdfwrite -dPDFSETTINGS=/prepress Unknown_*.pdf')) - } - } else { - error_message <- paste("Can't find",files,"!") - print(error_message) - stop(error_message) + if (args$selecttic) { + cat("\n\tProcessing TIC(s) from XCMS files...\n") + b <- getTIC2s(files = files, xset = xset, rt = "raw", pdfname = "TICs_raw.pdf") + cat("TIC(s) created...\n") + } + if (args$selecteic) { + cat("\n\tProcessing EIC(s) from XCMS files...\n") + cat(length(unknarg), "unknown(s) will be process !\n") + plotUnknowns(resGC = resGC, unkn = unknarg, DB = DBgc, fileFrom = "singlefile") + cat("EIC(s) created...\n") + # Merge pdf + system(paste("gs -o GCMS_EIC.pdf -sDEVICE=pdfwrite -dPDFSETTINGS=/prepress Unknown_*.pdf")) } + } else { + error_message <- paste("Can't find", files, "!") + print(error_message) + stop(error_message) + } } -if(!is.null(zipfile)) { - files <- getMSFiles(directory) - if(!is.null(files)) { - if(args$selectbpc) { - cat("\n\tProcessing BPC(s) from raw files...\n") - c <- getBPC2s(files = files, rt="raw", pdfname="BPCs_raw.pdf") - cat("BPC(s) created...\n") - } - if(args$selecttic) { - cat("\n\tProcessing TIC(s) from raw files...\n") - b <- getTIC2s(files = files, rt="raw", pdfname="TICs_raw.pdf") - cat("TIC(s) created...\n") - } - if(args$selecteic) { - cat("\n\tProcessing EIC(s) from XCMS files...\n") - cat(length(unknarg),"unknown(s) will be process !\n") - plotUnknowns(resGC=resGC, unkn=unknarg, DB=DBgc, fileFrom="zipfile") - cat("EIC(s) created...\n") - #Merge pdf - system(paste('gs -sOutputFile=GCMS_EIC.pdf -sDEVICE=pdfwrite -dPDFSETTINGS=/prepress Unknown_*.pdf')) - } - } else { - error_message <- paste("Can't find",files,"!") - print(error_message) - stop(error_message) - } +if (!is.null(zipfile)) { + files <- getMSFiles(directory) + if (!is.null(files)) { + if (args$selectbpc) { + cat("\n\tProcessing BPC(s) from raw files...\n") + c <- getBPC2s(files = files, rt = "raw", pdfname = "BPCs_raw.pdf") + cat("BPC(s) created...\n") + } + if (args$selecttic) { + cat("\n\tProcessing TIC(s) from raw files...\n") + b <- getTIC2s(files = files, rt = "raw", pdfname = "TICs_raw.pdf") + cat("TIC(s) created...\n") + } + if (args$selecteic) { + cat("\n\tProcessing EIC(s) from XCMS files...\n") + cat(length(unknarg), "unknown(s) will be process !\n") + plotUnknowns(resGC = resGC, unkn = unknarg, DB = DBgc, fileFrom = "zipfile") + cat("EIC(s) created...\n") + # Merge pdf + system(paste("gs -sOutputFile=GCMS_EIC.pdf -sDEVICE=pdfwrite -dPDFSETTINGS=/prepress Unknown_*.pdf")) + } + } else { + error_message <- paste("Can't find", files, "!") + print(error_message) + stop(error_message) + } } -cat("\nEnd of '", modNamC, "' Galaxy module call: ", as.character(Sys.time()), "\n", sep = "") \ No newline at end of file +cat("\nEnd of '", modNamC, "' Galaxy module call: ", as.character(Sys.time()), "\n", sep = "") diff --git a/tools/metams/metaMS_runGC/metaMS_runGC.r b/tools/metams/metaMS_runGC/metaMS_runGC.r index 7fbdfba3d..070c560a4 100644 --- a/tools/metams/metaMS_runGC/metaMS_runGC.r +++ b/tools/metams/metaMS_runGC/metaMS_runGC.r @@ -1,175 +1,181 @@ #!/usr/bin/env Rscript # metaMS_runGC.r version="3.0.0" -#created by Yann GUITTON and updated by Julien SAINT-VANNE -#use RI options + add try on plotUnknown add session Info -#use make.names in sampleMetadata to avoid issues with files names +# created by Yann GUITTON and updated by Julien SAINT-VANNE +# use RI options + add try on plotUnknown add session Info +# use make.names in sampleMetadata to avoid issues with files names # ----- LOG FILE ----- -#log_file=file("log.txt", open = "wt") -#sink(log_file) -#sink(log_file, type = "output") +# log_file=file("log.txt", open = "wt") +# sink(log_file) +# sink(log_file, type = "output") # ----- PACKAGE ----- cat("\tSESSION INFO\n\n") -#Import the different functions and packages +# Import the different functions and packages source_local <- function(fname) { - argv <- commandArgs(trailingOnly = FALSE) - base_dir <- dirname(substring(argv[grep("--file=", argv)], 8)) - source(paste(base_dir, fname, sep="/")) + argv <- commandArgs(trailingOnly = FALSE) + base_dir <- dirname(substring(argv[grep("--file=", argv)], 8)) + source(paste(base_dir, fname, sep = "/")) } source_local("lib_metams.r") -pkgs <- c("metaMS","stringr","batch","CAMERA") #"batch" necessary for parseCommandArgs function +pkgs <- c("metaMS", "stringr", "batch", "CAMERA") # "batch" necessary for parseCommandArgs function loadAndDisplayPackages(pkgs) cat("\n\n") modNamC <- "metaMS:runGC" ## module name -cat("\nStart of the '", modNamC, "' Galaxy module call: ", format(Sys.time(), "%a %d %b %Y %X"), "\n", sep="") +cat("\nStart of the '", modNamC, "' Galaxy module call: ", format(Sys.time(), "%a %d %b %Y %X"), "\n", sep = "") # ----- PROCESSING INFILE ----- cat("\n\n\tARGUMENTS PROCESSING INFO\n\n") -args = parseCommandArgs(evaluate=FALSE) #interpretation of arguments given in command line as an R list of objects -#write.table(as.matrix(args), col.names=F, quote=F, sep='\t\t') +args <- parseCommandArgs(evaluate = FALSE) # interpretation of arguments given in command line as an R list of objects +# write.table(as.matrix(args), col.names=F, quote=F, sep='\t\t') print(cbind(value = unlist(args))) # ----- PROCESSING INFILE ----- cat("\n\n\tARGUMENTS PROCESSING INFO\n\n") # Saving the specific parameters -#RI parameter -if (args$ri!="NULL"){ - RIarg=read.table(args$ri) - if (ncol(RIarg) < 2) RIarg=read.table(args$ri, h=T, sep=";") - if (ncol(RIarg) < 2) RIarg=read.table(args$ri, h=T, sep="\t") - if (ncol(RIarg) < 2) RIarg=read.table(args$ri, h=T, sep=",") - if (ncol(RIarg) < 2) { - error_message="Your RI file seems not well formatted. The column separators accepted are ; , and tabulation" - print(error_message) - stop(error_message) - } - #to do check real column names - colnames(RIarg)<-c("rt","RI") +# RI parameter +if (args$ri != "NULL") { + RIarg <- read.table(args$ri) + if (ncol(RIarg) < 2) RIarg <- read.table(args$ri, h = T, sep = ";") + if (ncol(RIarg) < 2) RIarg <- read.table(args$ri, h = T, sep = "\t") + if (ncol(RIarg) < 2) RIarg <- read.table(args$ri, h = T, sep = ",") + if (ncol(RIarg) < 2) { + error_message <- "Your RI file seems not well formatted. The column separators accepted are ; , and tabulation" + print(error_message) + stop(error_message) + } + # to do check real column names + colnames(RIarg) <- c("rt", "RI") } else { - RIarg = NULL + RIarg <- NULL } -#RIshift parameter -if (args$rishift!="none"){ - RIshift=args$rishift +# RIshift parameter +if (args$rishift != "none") { + RIshift <- args$rishift } else { - RIshift = "none" + RIshift <- "none" } -#Personal database parameter -if (args$db!="NULL"){ - DBgc=args$db +# Personal database parameter +if (args$db != "NULL") { + DBgc <- args$db } else { - DBgc = NULL -} - -#settings process -if (args$settings=="default") { - cat("Using default parameters") - data(FEMsettings) - if (args$rtrange[1]!="NULL") { - rtrange=args$rtrange - } else { - rtrange=NULL - } - - if (!is.null(DBgc)){ - manual <- read.msp(DBgc) - DBgc <- createSTDdbGC(stdInfo = NULL, settings = TSQXLS.GC, manualDB = manual) - } - - #use RI instead of rt for time comparison vs DB - if (RIshift!="none"){ - TSQXLS.GC@match2DB.timeComparison<-"RI" - TSQXLS.GC@match2DB.RIdiff<-as.numeric(RIshift) - TSQXLS.GC@betweenSamples.timeComparison<-"RI" - TSQXLS.GC@betweenSamples.RIdiff<-as.numeric(RIshift) - } - nSlaves=args$nSlaves + DBgc <- NULL } -if (args$settings=="User_defined") { - cat("Using user's parameters\n") - rtdiffparam=args$rtdiff - minfeatparam=args$minfeat - simthreshparam=args$simthreshold - minclassfractionparam=args$minclassfraction - minclasssizeparam=args$minclasssize - - if (args$rtrange!="NULL") { - rtrange=args$rtrange - cat("rtrange= ",rtrange) - } else { - rtrange=NULL - cat("rtrange= ",rtrange) - } +# settings process +if (args$settings == "default") { + cat("Using default parameters") + data(FEMsettings) + if (args$rtrange[1] != "NULL") { + rtrange <- args$rtrange + } else { + rtrange <- NULL + } + + if (!is.null(DBgc)) { + manual <- read.msp(DBgc) + DBgc <- createSTDdbGC(stdInfo = NULL, settings = TSQXLS.GC, manualDB = manual) + } + + # use RI instead of rt for time comparison vs DB + if (RIshift != "none") { + TSQXLS.GC@match2DB.timeComparison <- "RI" + TSQXLS.GC@match2DB.RIdiff <- as.numeric(RIshift) + TSQXLS.GC@betweenSamples.timeComparison <- "RI" + TSQXLS.GC@betweenSamples.RIdiff <- as.numeric(RIshift) + } + nSlaves <- args$nSlaves +} - nSlaves=args$nSlaves -# peakpicking parameters are not needed anymore xcms does it before runGC but keep it for protocolName and chrom type (2020-06-18 Yann Guitton) - GALAXY.GC <- - metaMSsettings("protocolName" = "GALAXY.GC", - "chrom" = "GC", - PeakPicking = list( - method = "matchedFilter", - step = 0.5, - steps = 2, - mzdiff = .5, - fwhm = 5, - snthresh = 2, - max = 500), - CAMERA = list(perfwhm = 1)) - - metaSetting(GALAXY.GC, "DBconstruction") <- list( - minintens = 0.0, - rttol = rtdiffparam, - intensityMeasure = "maxo", - DBthreshold = .80, - minfeat = minfeatparam) - metaSetting(GALAXY.GC, "match2DB") <- list( - simthresh = simthreshparam, - timeComparison = "rt", - rtdiff = rtdiffparam, - RIdiff = 5, - minfeat = minfeatparam) - - #to used if contaminant filter - - # metaSetting(GALAXY.GC, "matchIrrelevants") <- list( - # irrelevantClasses = c("Bleeding", "Plasticizers"), - # timeComparison = "RI", - # RIdiff = RIdiffparam, - # rtdiff = rtdiffparam, - # simthresh = simthreshparam) - - metaSetting(GALAXY.GC, "betweenSamples") <- list( - min.class.fraction = minclassfractionparam, - min.class.size = minclasssizeparam, - timeComparison = "rt", - rtdiff = rtdiffparam, - RIdiff = 2, - simthresh = simthreshparam) - - #ONLY use RI instead of rt for time comparison vs DB or samples - if (RIshift!="none"){ - GALAXY.GC@match2DB.timeComparison<-"RI" - GALAXY.GC@match2DB.RIdiff<-as.numeric(RIshift) - GALAXY.GC@betweenSamples.timeComparison<-"RI" - GALAXY.GC@betweenSamples.RIdiff<-as.numeric(RIshift) - } - - if (!is.null(DBgc)){ - manual <- read.msp(DBgc) - DBgc <- createSTDdbGC(stdInfo = NULL, settings = GALAXY.GC, manualDB = manual) - } +if (args$settings == "User_defined") { + cat("Using user's parameters\n") + rtdiffparam <- args$rtdiff + minfeatparam <- args$minfeat + simthreshparam <- args$simthreshold + minclassfractionparam <- args$minclassfraction + minclasssizeparam <- args$minclasssize + + if (args$rtrange != "NULL") { + rtrange <- args$rtrange + cat("rtrange= ", rtrange) + } else { + rtrange <- NULL + cat("rtrange= ", rtrange) + } + + nSlaves <- args$nSlaves + # peakpicking parameters are not needed anymore xcms does it before runGC but keep it for protocolName and chrom type (2020-06-18 Yann Guitton) + GALAXY.GC <- + metaMSsettings( + "protocolName" = "GALAXY.GC", + "chrom" = "GC", + PeakPicking = list( + method = "matchedFilter", + step = 0.5, + steps = 2, + mzdiff = .5, + fwhm = 5, + snthresh = 2, + max = 500 + ), + CAMERA = list(perfwhm = 1) + ) + + metaSetting(GALAXY.GC, "DBconstruction") <- list( + minintens = 0.0, + rttol = rtdiffparam, + intensityMeasure = "maxo", + DBthreshold = .80, + minfeat = minfeatparam + ) + metaSetting(GALAXY.GC, "match2DB") <- list( + simthresh = simthreshparam, + timeComparison = "rt", + rtdiff = rtdiffparam, + RIdiff = 5, + minfeat = minfeatparam + ) + + # to used if contaminant filter + + # metaSetting(GALAXY.GC, "matchIrrelevants") <- list( + # irrelevantClasses = c("Bleeding", "Plasticizers"), + # timeComparison = "RI", + # RIdiff = RIdiffparam, + # rtdiff = rtdiffparam, + # simthresh = simthreshparam) + + metaSetting(GALAXY.GC, "betweenSamples") <- list( + min.class.fraction = minclassfractionparam, + min.class.size = minclasssizeparam, + timeComparison = "rt", + rtdiff = rtdiffparam, + RIdiff = 2, + simthresh = simthreshparam + ) + + # ONLY use RI instead of rt for time comparison vs DB or samples + if (RIshift != "none") { + GALAXY.GC@match2DB.timeComparison <- "RI" + GALAXY.GC@match2DB.RIdiff <- as.numeric(RIshift) + GALAXY.GC@betweenSamples.timeComparison <- "RI" + GALAXY.GC@betweenSamples.RIdiff <- as.numeric(RIshift) + } + + if (!is.null(DBgc)) { + manual <- read.msp(DBgc) + DBgc <- createSTDdbGC(stdInfo = NULL, settings = GALAXY.GC, manualDB = manual) + } } @@ -190,155 +196,167 @@ cat("\n\tMAIN PROCESSING INFO\n") cat("\t\tCOMPUTE\n\n") -#runGC accept either a list of files a zip folder or an xset object from xcms.xcmsSet tool -#From xset is an .RData file necessary to use the xcmsSet object generated by xcms.xcmsSet given by previous tools -if (!is.null(args$singlefile_galaxyPath)){ - cat("Loading datas from XCMS files...\n") - load(args$singlefile_galaxyPath) - - #Transform XCMS object if needed - if(!exists("xset")) { - if(exists("xdata")) { - xset<-getxcmsSetObject(xdata) - } else { - error_message="no xset and no xdata... Probably a problem" - print(error_message) - stop(error_message) - } +# runGC accept either a list of files a zip folder or an xset object from xcms.xcmsSet tool +# From xset is an .RData file necessary to use the xcmsSet object generated by xcms.xcmsSet given by previous tools +if (!is.null(args$singlefile_galaxyPath)) { + cat("Loading datas from XCMS files...\n") + load(args$singlefile_galaxyPath) + + # Transform XCMS object if needed + if (!exists("xset")) { + if (exists("xdata")) { + xset <- getxcmsSetObject(xdata) + } else { + error_message <- "no xset and no xdata... Probably a problem" + print(error_message) + stop(error_message) } - #Verify that there are more than 1 file (can't run metaMS on only 1 file) - if(length(rownames(xdata@phenoData)) < 2){ - error_message="You need more than 1 file to be able to run metaMS" - print(error_message) - stop(error_message) + } + # Verify that there are more than 1 file (can't run metaMS on only 1 file) + if (length(rownames(xdata@phenoData)) < 2) { + error_message <- "You need more than 1 file to be able to run metaMS" + print(error_message) + stop(error_message) + } + + # xset from xcms.xcmsSet is not well formatted for metaMS this function do the formatting + if (class(xset) == "xcmsSet") { + if (length(xset@rt$raw) > 1) { + # create an exceptable list of xset for metaMS + xset.l <- vector("list", length(xset@rt$raw)) + for (i in 1:length(xset@rt$raw)) { + xset.l[[i]] <- new("xcmsSet") + xset.l[[i]]@peaks <- xset@peaks[which(xset@peaks[, "sample"] == i), ] + df <- data.frame(class = xset@phenoData[i, ]) + rownames(df) <- rownames(xset@phenoData)[i] + xset.l[[i]]@phenoData <- df + xset.l[[i]]@rt$raw <- xset@rt$raw[[i]] + xset.l[[i]]@rt$corrected <- xset@rt$corrected[[i]] + xset.l[[i]]@filepaths <- xset@filepaths[i] + xset.l[[i]]@profinfo <- xset@profinfo + } + } else { + xset.l <- xset } - #xset from xcms.xcmsSet is not well formatted for metaMS this function do the formatting - if (class(xset)=="xcmsSet"){ - if (length(xset@rt$raw)>1){ - #create an exceptable list of xset for metaMS - xset.l<-vector("list",length(xset@rt$raw)) - for (i in 1:length(xset@rt$raw)){ - xset.l[[i]]<-new("xcmsSet") - xset.l[[i]]@peaks<-xset@peaks[which(xset@peaks[,"sample"]==i),] - df<-data.frame(class=xset@phenoData[i,]) - rownames(df)<-rownames(xset@phenoData)[i] - xset.l[[i]]@phenoData<-df - xset.l[[i]]@rt$raw<-xset@rt$raw[[i]] - xset.l[[i]]@rt$corrected<-xset@rt$corrected[[i]] - xset.l[[i]]@filepaths<-xset@filepaths[i] - xset.l[[i]]@profinfo<-xset@profinfo - } - } else { - xset.l<-xset + # create sampleMetadata, get sampleMetadata and class + sampleMetadata <- xset@phenoData + colnames(sampleMetadata) <- c("sampleMetadata", "sample_group", "class") + sampleMetadata <- sampleMetadata[, -2] + row.names(sampleMetadata) <- NULL + samples <- xset@filepaths + } else { + xset <- NULL + } + if (args$settings == "default") { + settingslist <- TSQXLS.GC + if (class(xset.l[[1]]) != "xsAnnotate") { + cat("Process xsAnnotate with CAMERA package...\n") + xsetCAM <- lapply( + xset.l, + function(x) { + y <- xsAnnotate(x, sample = 1) + capture.output(z <- groupFWHM(y, perfwhm = settingslist@CAMERA$perfwhm), + file = NULL + ) + z } - - #create sampleMetadata, get sampleMetadata and class - sampleMetadata<-xset@phenoData - colnames(sampleMetadata) <- c("sampleMetadata","sample_group","class") - sampleMetadata<-sampleMetadata[,-2] - row.names(sampleMetadata)<-NULL - samples<-xset@filepaths + ) } else { - xset<-NULL + xsetCAM <- xset.l } - if(args$settings == "default") { - settingslist=TSQXLS.GC - if (class(xset.l[[1]])!="xsAnnotate"){ - cat("Process xsAnnotate with CAMERA package...\n") - xsetCAM<-lapply(xset.l, - function(x) { - y <- xsAnnotate(x, sample = 1) - capture.output(z <- groupFWHM(y, perfwhm = settingslist@CAMERA$perfwhm), - file = NULL) - z}) - } else { - xsetCAM <- xset.l - } - - #default settings for GC from Wehrens et al - cat("Process runGC with metaMS package...\n\n") - print(str(TSQXLS.GC)) - resGC<-runGC(xset=xsetCAM,settings=TSQXLS.GC, rtrange=rtrange, DB= DBgc, removeArtefacts = TRUE, - findUnknowns = TRUE, returnXset = TRUE, RIstandards = RIarg, nSlaves = nSlaves) + + # default settings for GC from Wehrens et al + cat("Process runGC with metaMS package...\n\n") + print(str(TSQXLS.GC)) + resGC <- runGC( + xset = xsetCAM, settings = TSQXLS.GC, rtrange = rtrange, DB = DBgc, removeArtefacts = TRUE, + findUnknowns = TRUE, returnXset = TRUE, RIstandards = RIarg, nSlaves = nSlaves + ) + } else { + if (args$settings == "User_defined") { + settingslist <- GALAXY.GC + if (class(xset.l[[1]]) != "xsAnnotate") { + cat("Process xsAnnotate with CAMERA package...\n") + xsetCAM <- lapply( + xset.l, + function(x) { + y <- xsAnnotate(x, sample = 1) + capture.output(z <- groupFWHM(y, perfwhm = settingslist@CAMERA$perfwhm), + file = NULL + ) + z + } + ) + } else { + xsetCAM <- xset.l + } + + # user settings for GC + cat("Process runGC with metaMS package...\n\n") + print(str(GALAXY.GC)) + resGC <- runGC( + xset = xsetCAM, settings = GALAXY.GC, rtrange = rtrange, DB = DBgc, removeArtefacts = TRUE, + findUnknowns = TRUE, returnXset = TRUE, RIstandards = RIarg, nSlaves = nSlaves + ) } else { - if(args$settings == "User_defined") { - settingslist=GALAXY.GC - if (class(xset.l[[1]])!="xsAnnotate") { - cat("Process xsAnnotate with CAMERA package...\n") - xsetCAM<-lapply(xset.l, - function(x) { - y <- xsAnnotate(x, sample = 1) - capture.output(z <- groupFWHM(y, perfwhm = settingslist@CAMERA$perfwhm), - file = NULL) - z}) - } else { - xsetCAM <- xset.l - } - - #user settings for GC - cat("Process runGC with metaMS package...\n\n") - print(str(GALAXY.GC)) - resGC<-runGC(xset=xsetCAM,settings=GALAXY.GC,rtrange = rtrange, DB= DBgc, removeArtefacts = TRUE, - findUnknowns = TRUE, returnXset = TRUE, RIstandards = RIarg, nSlaves = nSlaves) - } else { - error_message <- "There is no xset" - print(error_message) - stop(error_message) - } + error_message <- "There is no xset" + print(error_message) + stop(error_message) } + } } else { - error_message <- "No galaxy path entered" - print(error_message) - stop(error_message) + error_message <- "No galaxy path entered" + print(error_message) + stop(error_message) } # ----- EXPORT ----- -#peakTable ordered by rt +# peakTable ordered by rt cat("\nGenerating peakTable file") -peaktable<-getCorrectFileName(resGC$PeakTable,sampleMetadata) +peaktable <- getCorrectFileName(resGC$PeakTable, sampleMetadata) cat("\t.\t.") -write.table(peaktable, file="peaktable.tsv", sep="\t", row.names=FALSE) +write.table(peaktable, file = "peaktable.tsv", sep = "\t", row.names = FALSE) cat("\t.\tOK") -#peakTable for PCA -#dataMatrix +# peakTable for PCA +# dataMatrix cat("\nGenerating dataMatrix file") -dataMatrix<-cbind(Name=peaktable[,"Name"],peaktable[,(colnames(peaktable) %in% sampleMetadata[,1])]) -rownames(dataMatrix)<-NULL +dataMatrix <- cbind(Name = peaktable[, "Name"], peaktable[, (colnames(peaktable) %in% sampleMetadata[, 1])]) +rownames(dataMatrix) <- NULL cat("\t.\t.") -write.table(dataMatrix, file="dataMatrix.tsv", sep="\t", row.names=FALSE, quote=FALSE) +write.table(dataMatrix, file = "dataMatrix.tsv", sep = "\t", row.names = FALSE, quote = FALSE) cat("\t.\tOK") -#variableMetadata +# variableMetadata cat("\nGenerating variableMetadata file") -variableMetadata<-peaktable[,!(colnames(peaktable) %in% sampleMetadata[,1])] -rownames(variableMetadata)<-NULL +variableMetadata <- peaktable[, !(colnames(peaktable) %in% sampleMetadata[, 1])] +rownames(variableMetadata) <- NULL cat("\t.") -write.table(variableMetadata, file="variableMetadata.tsv", sep="\t", row.names=FALSE, quote=FALSE) +write.table(variableMetadata, file = "variableMetadata.tsv", sep = "\t", row.names = FALSE, quote = FALSE) cat("\t.\tOK") -#sampleMetadata +# sampleMetadata cat("\nGenerating sampleMetadata file") cat("\t.\t.") -write.table(sampleMetadata, file="sampleMetadata.tsv", sep="\t", row.names=FALSE, quote=FALSE) +write.table(sampleMetadata, file = "sampleMetadata.tsv", sep = "\t", row.names = FALSE, quote = FALSE) cat("\t.\tOK") -#peak spectrum as MSP for DB search -cat("\nGenerating",length(resGC$PseudoSpectra),"peakspectra in peakspectra.msp file\n") -write.msp(resGC$PseudoSpectra, file="peakspectra.msp", newFile = TRUE) +# peak spectrum as MSP for DB search +cat("\nGenerating", length(resGC$PseudoSpectra), "peakspectra in peakspectra.msp file\n") +write.msp(resGC$PseudoSpectra, file = "peakspectra.msp", newFile = TRUE) -#saving R data in .Rdata file to save the variables used in the present tool +# saving R data in .Rdata file to save the variables used in the present tool objects2save <- c("resGC", "xset", "singlefile", "zipfile", "DBgc") save(list = objects2save[objects2save %in% ls()], file = "runGC.RData") cat("\nEnd of '", modNamC, "' Galaxy module call: ", as.character(Sys.time()), "\n", sep = "") -#WARNING if user has CDF files (not yet good for plotting) -files <- paste("./",names(singlefile),sep="") -if(MSnbase:::isCdfFile(files)){ - warning_message <- "You have CDF files, for the moment you can't obtain plot after runGC! A new update will follow with the good correction" - warning(warning_message) - cat(paste("\n","/!\\Warning/!\\",warning_message,sep="\n")) +# WARNING if user has CDF files (not yet good for plotting) +files <- paste("./", names(singlefile), sep = "") +if (MSnbase:::isCdfFile(files)) { + warning_message <- "You have CDF files, for the moment you can't obtain plot after runGC! A new update will follow with the good correction" + warning(warning_message) + cat(paste("\n", "/!\\Warning/!\\", warning_message, sep = "\n")) } diff --git a/tools/metams/scripts/lib_metams.r b/tools/metams/scripts/lib_metams.r index 2911add9e..10c356365 100644 --- a/tools/metams/scripts/lib_metams.r +++ b/tools/metams/scripts/lib_metams.r @@ -4,503 +4,517 @@ # author Yann GUITTON Oniris Laberca 2015-2017 -#@author G. Le Corguille +# @author G. Le Corguille # This function will # - load the packages # - display the sessionInfo loadAndDisplayPackages <- function(pkgs) { - for(pkg in pkgs) suppressPackageStartupMessages( stopifnot( library(pkg, quietly=TRUE, logical.return=TRUE, character.only=TRUE))) - - sessioninfo = sessionInfo() - cat(sessioninfo$R.version$version.string,"\n") - cat("Main packages:\n") - for (pkg in names(sessioninfo$otherPkgs)) { cat(paste(pkg,packageVersion(pkg)),"\t") }; cat("\n") - cat("Other loaded packages:\n") - for (pkg in names(sessioninfo$loadedOnly)) { cat(paste(pkg,packageVersion(pkg)),"\t") }; cat("\n") + for (pkg in pkgs) suppressPackageStartupMessages(stopifnot(library(pkg, quietly = TRUE, logical.return = TRUE, character.only = TRUE))) + + sessioninfo <- sessionInfo() + cat(sessioninfo$R.version$version.string, "\n") + cat("Main packages:\n") + for (pkg in names(sessioninfo$otherPkgs)) { + cat(paste(pkg, packageVersion(pkg)), "\t") + } + cat("\n") + cat("Other loaded packages:\n") + for (pkg in names(sessioninfo$loadedOnly)) { + cat(paste(pkg, packageVersion(pkg)), "\t") + } + cat("\n") } -#This function list the compatible files within the directory as xcms did -#@author Gildas Le Corguille lecorguille@sb-roscoff.fr ABiMS TEAM -getMSFiles <- function (directory) { - filepattern <- c("[Cc][Dd][Ff]", "[Nn][Cc]", "([Mm][Zz])?[Xx][Mm][Ll]","[Mm][Zz][Dd][Aa][Tt][Aa]", "[Mm][Zz][Mm][Ll]") - filepattern <- paste(paste("\\.", filepattern, "$", sep=""),collapse="|") - info <- file.info(directory) - listed <- list.files(directory[info$isdir], pattern=filepattern,recursive=TRUE, full.names=TRUE) - files <- c(directory[!info$isdir], listed) - exists <- file.exists(files) - files <- files[exists] - return(files) +# This function list the compatible files within the directory as xcms did +# @author Gildas Le Corguille lecorguille@sb-roscoff.fr ABiMS TEAM +getMSFiles <- function(directory) { + filepattern <- c("[Cc][Dd][Ff]", "[Nn][Cc]", "([Mm][Zz])?[Xx][Mm][Ll]", "[Mm][Zz][Dd][Aa][Tt][Aa]", "[Mm][Zz][Mm][Ll]") + filepattern <- paste(paste("\\.", filepattern, "$", sep = ""), collapse = "|") + info <- file.info(directory) + listed <- list.files(directory[info$isdir], pattern = filepattern, recursive = TRUE, full.names = TRUE) + files <- c(directory[!info$isdir], listed) + exists <- file.exists(files) + files <- files[exists] + return(files) } # This function retrieve a xset like object -#@author Gildas Le Corguille lecorguille@sb-roscoff.fr +# @author Gildas Le Corguille lecorguille@sb-roscoff.fr getxcmsSetObject <- function(xobject) { - # XCMS 1.x - if (class(xobject) == "xcmsSet") - return (xobject) - # XCMS 3.x - if (class(xobject) == "XCMSnExp") { - # Get the legacy xcmsSet object - suppressWarnings(xset <- as(xobject, 'xcmsSet')) - if (!is.null(xset@phenoData$sample_group)) - sampclass(xset) <- xset@phenoData$sample_group - else - sampclass(xset) <- "." - return (xset) + # XCMS 1.x + if (class(xobject) == "xcmsSet") { + return(xobject) + } + # XCMS 3.x + if (class(xobject) == "XCMSnExp") { + # Get the legacy xcmsSet object + suppressWarnings(xset <- as(xobject, "xcmsSet")) + if (!is.null(xset@phenoData$sample_group)) { + sampclass(xset) <- xset@phenoData$sample_group + } else { + sampclass(xset) <- "." } + return(xset) + } } -#J.Saint-Vanne -#Function to correct the file names which can be written like "..alg8.mzData" and we just want "alg8" -getCorrectFileName <- function(peaktable,sampleMetadata){ - - #Try to start for the first sample, avoid description of line with colnamesdontwant - i <- 1 - while(!(sampleMetadata[1,1] %in% strsplit(colnames(peaktable)[i],"\\.")[[1]])) { - if(i < length(peaktable)) { - i <- i + 1 - } else { - break - } +# J.Saint-Vanne +# Function to correct the file names which can be written like "..alg8.mzData" and we just want "alg8" +getCorrectFileName <- function(peaktable, sampleMetadata) { + # Try to start for the first sample, avoid description of line with colnamesdontwant + i <- 1 + while (!(sampleMetadata[1, 1] %in% strsplit(colnames(peaktable)[i], "\\.")[[1]])) { + if (i < length(peaktable)) { + i <- i + 1 + } else { + break } - #i now correspond to the first column with a sample - for(j in 1:(nrow(sampleMetadata))) { - col <- j + i-1 #minus 1 cause i is the good column to start and j start at 1 - if(col <= length(colnames(peaktable))) { - newname <- gsub("(^.*)(\\..*$)","\\1",colnames(peaktable)[col]) - if(newname != sampleMetadata[j,1]){ - #Correction for 2 points starting the name (I don't know why they are here...) - if(".." == gsub("(^\\.+)(.*)","\\1",newname)){ - newname <- sub("(^\\.+)(.*)","\\2",newname) - } - } - colnames(peaktable)[col] <- newname + } + # i now correspond to the first column with a sample + for (j in 1:(nrow(sampleMetadata))) { + col <- j + i - 1 # minus 1 cause i is the good column to start and j start at 1 + if (col <= length(colnames(peaktable))) { + newname <- gsub("(^.*)(\\..*$)", "\\1", colnames(peaktable)[col]) + if (newname != sampleMetadata[j, 1]) { + # Correction for 2 points starting the name (I don't know why they are here...) + if (".." == gsub("(^\\.+)(.*)", "\\1", newname)) { + newname <- sub("(^\\.+)(.*)", "\\2", newname) } + } + colnames(peaktable)[col] <- newname } - return(peaktable) + } + return(peaktable) } -#@author Gildas Le Corguille lecorguille@sb-roscoff.fr +# @author Gildas Le Corguille lecorguille@sb-roscoff.fr # This function get the raw file path from the arguments getRawfilePathFromArguments <- function(singlefile, zipfile, listArguments) { - if (!is.null(listArguments[["zipfile"]])) zipfile = listArguments[["zipfile"]] - if (!is.null(listArguments[["zipfilePositive"]])) zipfile = listArguments[["zipfilePositive"]] - if (!is.null(listArguments[["zipfileNegative"]])) zipfile = listArguments[["zipfileNegative"]] - - if (!is.null(listArguments[["singlefile_galaxyPath"]])) { - singlefile_galaxyPaths = listArguments[["singlefile_galaxyPath"]]; - singlefile_sampleNames = listArguments[["singlefile_sampleName"]] + if (!is.null(listArguments[["zipfile"]])) zipfile <- listArguments[["zipfile"]] + if (!is.null(listArguments[["zipfilePositive"]])) zipfile <- listArguments[["zipfilePositive"]] + if (!is.null(listArguments[["zipfileNegative"]])) zipfile <- listArguments[["zipfileNegative"]] + + if (!is.null(listArguments[["singlefile_galaxyPath"]])) { + singlefile_galaxyPaths <- listArguments[["singlefile_galaxyPath"]] + singlefile_sampleNames <- listArguments[["singlefile_sampleName"]] + } + if (!is.null(listArguments[["singlefile_galaxyPathPositive"]])) { + singlefile_galaxyPaths <- listArguments[["singlefile_galaxyPathPositive"]] + singlefile_sampleNames <- listArguments[["singlefile_sampleNamePositive"]] + } + if (!is.null(listArguments[["singlefile_galaxyPathNegative"]])) { + singlefile_galaxyPaths <- listArguments[["singlefile_galaxyPathNegative"]] + singlefile_sampleNames <- listArguments[["singlefile_sampleNameNegative"]] + } + if (exists("singlefile_galaxyPaths")) { + singlefile_galaxyPaths <- unlist(strsplit(singlefile_galaxyPaths, ",")) + singlefile_sampleNames <- unlist(strsplit(singlefile_sampleNames, ",")) + + singlefile <- NULL + for (singlefile_galaxyPath_i in seq(1:length(singlefile_galaxyPaths))) { + singlefile_galaxyPath <- singlefile_galaxyPaths[singlefile_galaxyPath_i] + singlefile_sampleName <- singlefile_sampleNames[singlefile_galaxyPath_i] + singlefile[[singlefile_sampleName]] <- singlefile_galaxyPath } - if (!is.null(listArguments[["singlefile_galaxyPathPositive"]])) { - singlefile_galaxyPaths = listArguments[["singlefile_galaxyPathPositive"]]; - singlefile_sampleNames = listArguments[["singlefile_sampleNamePositive"]] - } - if (!is.null(listArguments[["singlefile_galaxyPathNegative"]])) { - singlefile_galaxyPaths = listArguments[["singlefile_galaxyPathNegative"]]; - singlefile_sampleNames = listArguments[["singlefile_sampleNameNegative"]] - } - if (exists("singlefile_galaxyPaths")){ - singlefile_galaxyPaths = unlist(strsplit(singlefile_galaxyPaths,",")) - singlefile_sampleNames = unlist(strsplit(singlefile_sampleNames,",")) - - singlefile=NULL - for (singlefile_galaxyPath_i in seq(1:length(singlefile_galaxyPaths))) { - singlefile_galaxyPath=singlefile_galaxyPaths[singlefile_galaxyPath_i] - singlefile_sampleName=singlefile_sampleNames[singlefile_galaxyPath_i] - singlefile[[singlefile_sampleName]] = singlefile_galaxyPath - } - } - for (argument in c("zipfile","zipfilePositive","zipfileNegative","singlefile_galaxyPath","singlefile_sampleName","singlefile_galaxyPathPositive","singlefile_sampleNamePositive","singlefile_galaxyPathNegative","singlefile_sampleNameNegative")) { - listArguments[[argument]]=NULL - } - return(list(zipfile=zipfile, singlefile=singlefile, listArguments=listArguments)) + } + for (argument in c("zipfile", "zipfilePositive", "zipfileNegative", "singlefile_galaxyPath", "singlefile_sampleName", "singlefile_galaxyPathPositive", "singlefile_sampleNamePositive", "singlefile_galaxyPathNegative", "singlefile_sampleNameNegative")) { + listArguments[[argument]] <- NULL + } + return(list(zipfile = zipfile, singlefile = singlefile, listArguments = listArguments)) } -#@author Gildas Le Corguille lecorguille@sb-roscoff.fr +# @author Gildas Le Corguille lecorguille@sb-roscoff.fr # This function retrieve the raw file in the working directory # - if zipfile: unzip the file with its directory tree # - if singlefiles: set symlink with the good filename retrieveRawfileInTheWorkingDirectory <- function(singlefile, zipfile) { - if(!is.null(singlefile) && (length("singlefile")>0)) { - for (singlefile_sampleName in names(singlefile)) { - singlefile_galaxyPath = singlefile[[singlefile_sampleName]] - if(!file.exists(singlefile_galaxyPath)){ - error_message=paste("Cannot access the sample:",singlefile_sampleName,"located:",singlefile_galaxyPath,". Please, contact your administrator ... if you have one!") - print(error_message); stop(error_message) - } - file.symlink(singlefile_galaxyPath,singlefile_sampleName) - } - directory = "." + if (!is.null(singlefile) && (length("singlefile") > 0)) { + for (singlefile_sampleName in names(singlefile)) { + singlefile_galaxyPath <- singlefile[[singlefile_sampleName]] + if (!file.exists(singlefile_galaxyPath)) { + error_message <- paste("Cannot access the sample:", singlefile_sampleName, "located:", singlefile_galaxyPath, ". Please, contact your administrator ... if you have one!") + print(error_message) + stop(error_message) + } + file.symlink(singlefile_galaxyPath, singlefile_sampleName) + } + directory <- "." + } + if (!is.null(zipfile) && (zipfile != "")) { + if (!file.exists(zipfile)) { + error_message <- paste("Cannot access the Zip file:", zipfile, ". Please, contact your administrator ... if you have one!") + print(error_message) + stop(error_message) } - if(!is.null(zipfile) && (zipfile!="")) { - if(!file.exists(zipfile)){ - error_message=paste("Cannot access the Zip file:",zipfile,". Please, contact your administrator ... if you have one!") - print(error_message) - stop(error_message) - } - #list all file in the zip file - #zip_files=unzip(zipfile,list=T)[,"Name"] + # list all file in the zip file + # zip_files=unzip(zipfile,list=T)[,"Name"] - #unzip - suppressWarnings(unzip(zipfile, unzip="unzip")) + # unzip + suppressWarnings(unzip(zipfile, unzip = "unzip")) - #get the directory name - filesInZip=unzip(zipfile, list=T); - directories=unique(unlist(lapply(strsplit(filesInZip$Name,"/"), function(x) x[1]))); - directories=directories[!(directories %in% c("__MACOSX")) & file.info(directories)$isdir] - directory = "." - if (length(directories) == 1) directory = directories + # get the directory name + filesInZip <- unzip(zipfile, list = T) + directories <- unique(unlist(lapply(strsplit(filesInZip$Name, "/"), function(x) x[1]))) + directories <- directories[!(directories %in% c("__MACOSX")) & file.info(directories)$isdir] + directory <- "." + if (length(directories) == 1) directory <- directories - cat("files_root_directory\t",directory,"\n") - } - return (directory) + cat("files_root_directory\t", directory, "\n") + } + return(directory) } -##ADDITIONS FROM Y. Guitton -getBPC <- function(file,rtcor=NULL, ...) { - object <- xcmsRaw(file) - sel <- profRange(object, ...) - cbind(if (is.null(rtcor)) object@scantime[sel$scanidx] else rtcor ,xcms:::colMax(object@env$profile[sel$massidx,sel$scanidx,drop=FALSE])) +## ADDITIONS FROM Y. Guitton +getBPC <- function(file, rtcor = NULL, ...) { + object <- xcmsRaw(file) + sel <- profRange(object, ...) + cbind(if (is.null(rtcor)) object@scantime[sel$scanidx] else rtcor, xcms:::colMax(object@env$profile[sel$massidx, sel$scanidx, drop = FALSE])) } -getBPC2s <- function (files, xset = NULL, pdfname="BPCs.pdf", rt = c("raw","corrected"), scanrange=NULL) { - require(xcms) - - #create sampleMetadata, get sampleMetadata and class - if(!is.null(xset)) { - #When files come from XCMS3 directly before metaMS - sampleMetadata <- xset@phenoData +getBPC2s <- function(files, xset = NULL, pdfname = "BPCs.pdf", rt = c("raw", "corrected"), scanrange = NULL) { + require(xcms) + + # create sampleMetadata, get sampleMetadata and class + if (!is.null(xset)) { + # When files come from XCMS3 directly before metaMS + sampleMetadata <- xset@phenoData + } else { + # When files come from a zip directory with raw files before metaMS + sampleMetadata <- xcms:::phenoDataFromPaths(files) + } + class <- unique(sampleMetadata[, "class"]) # create phenoData like table + classnames <- vector("list", length(class)) + for (i in 1:length(class)) { + classnames[[i]] <- which(sampleMetadata[, "class"] == class[i]) + } + + N <- dim(sampleMetadata)[1] + BPC <- vector("list", N) + + for (j in 1:N) { + BPC[[j]] <- getBPC(files[j]) + # good for raw + # seems strange for corrected + # errors if scanrange used in xcmsSetgeneration + if (!is.null(xcmsSet) && rt == "corrected") { + rtcor <- xcmsSet@rt$corrected[[j]] } else { - #When files come from a zip directory with raw files before metaMS - sampleMetadata<-xcms:::phenoDataFromPaths(files) - } - class<-unique(sampleMetadata[,"class"]) #create phenoData like table - classnames<-vector("list",length(class)) - for (i in 1:length(class)){ - classnames[[i]]<-which( sampleMetadata[,"class"]==class[i]) + rtcor <- NULL } - - N <- dim(sampleMetadata)[1] - BPC <- vector("list",N) - - for (j in 1:N) { - BPC[[j]] <- getBPC(files[j]) - #good for raw - # seems strange for corrected - #errors if scanrange used in xcmsSetgeneration - if (!is.null(xcmsSet) && rt == "corrected"){ - rtcor <- xcmsSet@rt$corrected[[j]] - }else{ - rtcor <- NULL - } - BPC[[j]] <- getBPC(files[j],rtcor=rtcor) - } - - pdf(pdfname,w=16,h=10) - cols <- rainbow(N) - lty = 1:N - pch = 1:N - #search for max x and max y in BPCs - - xlim = range(sapply(BPC, function(x) range(x[,1]))) - ylim = range(sapply(BPC, function(x) range(x[,2]))) - - ylim = c(-ylim[2], ylim[2]) - - ##plot start - if (length(class)>2){ - for (k in 1:(length(class)-1)){ - for (l in (k+1):length(class)){ - cat(paste(class[k],"vs",class[l],sep=" ","\n")) - plot(0, 0, type="n", xlim = xlim/60, ylim = ylim, main = paste("Base Peak Chromatograms \n","BPCs_",class[k]," vs ",class[l], sep=""), xlab = "Retention Time (min)", ylab = "BPC") - colvect<-NULL - for (j in 1:length(classnames[[k]])) { - bpc <- BPC[[classnames[[k]][j]]] - # points(bpc[,1]/60, bpc[,2], col = cols[i], pch = pch[i], type="l") - points(bpc[,1]/60, bpc[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l") - colvect<-append(colvect,cols[classnames[[k]][j]]) - } - for (j in 1:length(classnames[[l]])) { - # i=class2names[j] - bpc <- BPC[[classnames[[l]][j]]] - points(bpc[,1]/60, -bpc[,2], col = cols[classnames[[l]][j]], pch = pch[classnames[[l]][j]], type="l") - colvect<-append(colvect,cols[classnames[[l]][j]]) - } - legend("topright",paste(gsub("(^.+)\\..*$","\\1",basename(files[c(classnames[[k]],classnames[[l]])]))), col = colvect, lty = lty, pch = pch) - } - } - }#end if length >2 - - if (length(class)==2){ - k=1 - l=2 - colvect<-NULL - plot(0, 0, type="n", xlim = xlim/60, ylim = ylim, main = paste("Base Peak Chromatograms \n","BPCs_",class[k],"vs",class[l], sep=""), xlab = "Retention Time (min)", ylab = "BPC") - + BPC[[j]] <- getBPC(files[j], rtcor = rtcor) + } + + pdf(pdfname, w = 16, h = 10) + cols <- rainbow(N) + lty <- 1:N + pch <- 1:N + # search for max x and max y in BPCs + + xlim <- range(sapply(BPC, function(x) range(x[, 1]))) + ylim <- range(sapply(BPC, function(x) range(x[, 2]))) + + ylim <- c(-ylim[2], ylim[2]) + + ## plot start + if (length(class) > 2) { + for (k in 1:(length(class) - 1)) { + for (l in (k + 1):length(class)) { + cat(paste(class[k], "vs", class[l], sep = " ", "\n")) + plot(0, 0, type = "n", xlim = xlim / 60, ylim = ylim, main = paste("Base Peak Chromatograms \n", "BPCs_", class[k], " vs ", class[l], sep = ""), xlab = "Retention Time (min)", ylab = "BPC") + colvect <- NULL for (j in 1:length(classnames[[k]])) { - bpc <- BPC[[classnames[[k]][j]]] - # points(bpc[,1]/60, bpc[,2], col = cols[i], pch = pch[i], type="l") - points(bpc[,1]/60, bpc[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l") - colvect<-append(colvect,cols[classnames[[k]][j]]) + bpc <- BPC[[classnames[[k]][j]]] + # points(bpc[,1]/60, bpc[,2], col = cols[i], pch = pch[i], type="l") + points(bpc[, 1] / 60, bpc[, 2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type = "l") + colvect <- append(colvect, cols[classnames[[k]][j]]) } for (j in 1:length(classnames[[l]])) { - # i=class2names[j] - bpc <- BPC[[classnames[[l]][j]]] - points(bpc[,1]/60, -bpc[,2], col = cols[classnames[[l]][j]], pch = pch[classnames[[l]][j]], type="l") - colvect<-append(colvect,cols[classnames[[l]][j]]) + # i=class2names[j] + bpc <- BPC[[classnames[[l]][j]]] + points(bpc[, 1] / 60, -bpc[, 2], col = cols[classnames[[l]][j]], pch = pch[classnames[[l]][j]], type = "l") + colvect <- append(colvect, cols[classnames[[l]][j]]) } - legend("topright",paste(gsub("(^.+)\\..*$","\\1",basename(files[c(classnames[[k]],classnames[[l]])]))), col = colvect, lty = lty, pch = pch) - }#end length ==2 - - if (length(class)==1){ - k=1 + legend("topright", paste(gsub("(^.+)\\..*$", "\\1", basename(files[c(classnames[[k]], classnames[[l]])]))), col = colvect, lty = lty, pch = pch) + } + } + } # end if length >2 + + if (length(class) == 2) { + k <- 1 + l <- 2 + colvect <- NULL + plot(0, 0, type = "n", xlim = xlim / 60, ylim = ylim, main = paste("Base Peak Chromatograms \n", "BPCs_", class[k], "vs", class[l], sep = ""), xlab = "Retention Time (min)", ylab = "BPC") + + for (j in 1:length(classnames[[k]])) { + bpc <- BPC[[classnames[[k]][j]]] + # points(bpc[,1]/60, bpc[,2], col = cols[i], pch = pch[i], type="l") + points(bpc[, 1] / 60, bpc[, 2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type = "l") + colvect <- append(colvect, cols[classnames[[k]][j]]) + } + for (j in 1:length(classnames[[l]])) { + # i=class2names[j] + bpc <- BPC[[classnames[[l]][j]]] + points(bpc[, 1] / 60, -bpc[, 2], col = cols[classnames[[l]][j]], pch = pch[classnames[[l]][j]], type = "l") + colvect <- append(colvect, cols[classnames[[l]][j]]) + } + legend("topright", paste(gsub("(^.+)\\..*$", "\\1", basename(files[c(classnames[[k]], classnames[[l]])]))), col = colvect, lty = lty, pch = pch) + } # end length ==2 - ylim = range(sapply(BPC, function(x) range(x[,2]))) + if (length(class) == 1) { + k <- 1 - colvect<-NULL - plot(0, 0, type="n", xlim = xlim/60, ylim = ylim, main = paste("Base Peak Chromatograms \n","BPCs_",class[k], sep=""), xlab = "Retention Time (min)", ylab = "BPC") + ylim <- range(sapply(BPC, function(x) range(x[, 2]))) - for (j in 1:length(classnames[[k]])) { - bpc <- BPC[[classnames[[k]][j]]] - # points(bpc[,1]/60, bpc[,2], col = cols[i], pch = pch[i], type="l") - points(bpc[,1]/60, bpc[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l") - colvect<-append(colvect,cols[classnames[[k]][j]]) - } - legend("topright",paste(gsub("(^.+)\\..*$","\\1",basename(files[c(classnames[[k]])]))), col = colvect, lty = lty, pch = pch) - }#end length ==1 - dev.off() + colvect <- NULL + plot(0, 0, type = "n", xlim = xlim / 60, ylim = ylim, main = paste("Base Peak Chromatograms \n", "BPCs_", class[k], sep = ""), xlab = "Retention Time (min)", ylab = "BPC") + + for (j in 1:length(classnames[[k]])) { + bpc <- BPC[[classnames[[k]][j]]] + # points(bpc[,1]/60, bpc[,2], col = cols[i], pch = pch[i], type="l") + points(bpc[, 1] / 60, bpc[, 2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type = "l") + colvect <- append(colvect, cols[classnames[[k]][j]]) + } + legend("topright", paste(gsub("(^.+)\\..*$", "\\1", basename(files[c(classnames[[k]])]))), col = colvect, lty = lty, pch = pch) + } # end length ==1 + dev.off() } -getTIC <- function(file,rtcor=NULL) { - object <- xcmsRaw(file) - cbind(if (is.null(rtcor)) object@scantime else rtcor, rawEIC(object,mzrange=range(object@env$mz))$intensity) +getTIC <- function(file, rtcor = NULL) { + object <- xcmsRaw(file) + cbind(if (is.null(rtcor)) object@scantime else rtcor, rawEIC(object, mzrange = range(object@env$mz))$intensity) } ## overlay TIC from all files in current folder or from xcmsSet, create pdf -getTIC2s <- function(files, xset=NULL, pdfname="TICs.pdf", rt=c("raw","corrected")) { - require(xcms) - - #create sampleMetadata, get sampleMetadata and class - if(!is.null(xset)){ - #When files come from XCMS3 before metaMS treatment - sampleMetadata<-xset@phenoData +getTIC2s <- function(files, xset = NULL, pdfname = "TICs.pdf", rt = c("raw", "corrected")) { + require(xcms) + + # create sampleMetadata, get sampleMetadata and class + if (!is.null(xset)) { + # When files come from XCMS3 before metaMS treatment + sampleMetadata <- xset@phenoData + } else { + # When files come from a zip directory with raw files before metaMS + sampleMetadata <- xcms:::phenoDataFromPaths(files) + } + class <- as.vector(levels(sampleMetadata[, "class"])) # create phenoData like table + classnames <- vector("list", length(class)) + for (i in 1:length(class)) { + classnames[[i]] <- which(sampleMetadata[, "class"] == class[i]) + } + + N <- dim(sampleMetadata)[1] + TIC <- vector("list", N) + + for (i in 1:N) { + if (!is.null(xcmsSet) && rt == "corrected") { + rtcor <- xcmsSet@rt$corrected[[i]] } else { - #When files come from a zip directory with raw files before metaMS - sampleMetadata<-xcms:::phenoDataFromPaths(files) + rtcor <- NULL } - class<-as.vector(levels(sampleMetadata[,"class"])) #create phenoData like table - classnames<-vector("list",length(class)) - for (i in 1:length(class)){ - classnames[[i]]<-which( sampleMetadata[,"class"]==class[i]) - } - - N <- dim(sampleMetadata)[1] - TIC <- vector("list",N) - - for (i in 1:N) { - if (!is.null(xcmsSet) && rt == "corrected") - rtcor <- xcmsSet@rt$corrected[[i]] - else - rtcor <- NULL - TIC[[i]] <- getTIC(files[i],rtcor=rtcor) - } - - pdf(pdfname,w=16,h=10) - cols <- rainbow(N) - lty = 1:N - pch = 1:N - #search for max x and max y in TICs - xlim = range(sapply(TIC, function(x) range(x[,1]))) - ylim = range(sapply(TIC, function(x) range(x[,2]))) - ylim = c(-ylim[2], ylim[2]) - - ##plot start - if (length(class)>2){ - for (k in 1:(length(class)-1)){ - for (l in (k+1):length(class)){ - cat(paste(class[k],"vs",class[l],"\n",sep=" ")) - plot(0, 0, type="n", xlim = xlim/60, ylim = ylim, main = paste("Total Ion Chromatograms \n","TICs_",class[k]," vs ",class[l], sep=""), xlab = "Retention Time (min)", ylab = "TIC") - colvect<-NULL - for (j in 1:length(classnames[[k]])) { - tic <- TIC[[classnames[[k]][j]]] - # points(tic[,1]/60, tic[,2], col = cols[i], pch = pch[i], type="l") - points(tic[,1]/60, tic[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l") - colvect<-append(colvect,cols[classnames[[k]][j]]) - } - for (j in 1:length(classnames[[l]])) { - # i=class2names[j] - tic <- TIC[[classnames[[l]][j]]] - points(tic[,1]/60, -tic[,2], col = cols[classnames[[l]][j]], pch = pch[classnames[[l]][j]], type="l") - colvect<-append(colvect,cols[classnames[[l]][j]]) - } - legend("topright",paste(gsub("(^.+)\\..*$","\\1",basename(files[c(classnames[[k]],classnames[[l]])]))), col = colvect, lty = lty, pch = pch) - } - } - }#end if length >2 - - if (length(class)==2){ - k=1 - l=2 - plot(0, 0, type="n", xlim = xlim/60, ylim = ylim, main = paste("Total Ion Chromatograms \n","TICs_",class[k],"vs",class[l], sep=""), xlab = "Retention Time (min)", ylab = "TIC") - colvect<-NULL + TIC[[i]] <- getTIC(files[i], rtcor = rtcor) + } + + pdf(pdfname, w = 16, h = 10) + cols <- rainbow(N) + lty <- 1:N + pch <- 1:N + # search for max x and max y in TICs + xlim <- range(sapply(TIC, function(x) range(x[, 1]))) + ylim <- range(sapply(TIC, function(x) range(x[, 2]))) + ylim <- c(-ylim[2], ylim[2]) + + ## plot start + if (length(class) > 2) { + for (k in 1:(length(class) - 1)) { + for (l in (k + 1):length(class)) { + cat(paste(class[k], "vs", class[l], "\n", sep = " ")) + plot(0, 0, type = "n", xlim = xlim / 60, ylim = ylim, main = paste("Total Ion Chromatograms \n", "TICs_", class[k], " vs ", class[l], sep = ""), xlab = "Retention Time (min)", ylab = "TIC") + colvect <- NULL for (j in 1:length(classnames[[k]])) { - tic <- TIC[[classnames[[k]][j]]] - # points(tic[,1]/60, tic[,2], col = cols[i], pch = pch[i], type="l") - points(tic[,1]/60, tic[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l") - colvect<-append(colvect,cols[classnames[[k]][j]]) + tic <- TIC[[classnames[[k]][j]]] + # points(tic[,1]/60, tic[,2], col = cols[i], pch = pch[i], type="l") + points(tic[, 1] / 60, tic[, 2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type = "l") + colvect <- append(colvect, cols[classnames[[k]][j]]) } for (j in 1:length(classnames[[l]])) { - # i=class2names[j] - tic <- TIC[[classnames[[l]][j]]] - points(tic[,1]/60, -tic[,2], col = cols[classnames[[l]][j]], pch = pch[classnames[[l]][j]], type="l") - colvect<-append(colvect,cols[classnames[[l]][j]]) - } - legend("topright",paste(gsub("(^.+)\\..*$","\\1",basename(files[c(classnames[[k]],classnames[[l]])]))), col = colvect, lty = lty, pch = pch) - }#end length ==2 - - if (length(class)==1){ - k=1 - ylim = range(sapply(TIC, function(x) range(x[,2]))) - plot(0, 0, type="n", xlim = xlim/60, ylim = ylim, main = paste("Total Ion Chromatograms \n","TICs_",class[k], sep=""), xlab = "Retention Time (min)", ylab = "TIC") - colvect<-NULL - for (j in 1:length(classnames[[k]])) { - tic <- TIC[[classnames[[k]][j]]] - # points(tic[,1]/60, tic[,2], col = cols[i], pch = pch[i], type="l") - points(tic[,1]/60, tic[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l") - colvect<-append(colvect,cols[classnames[[k]][j]]) + # i=class2names[j] + tic <- TIC[[classnames[[l]][j]]] + points(tic[, 1] / 60, -tic[, 2], col = cols[classnames[[l]][j]], pch = pch[classnames[[l]][j]], type = "l") + colvect <- append(colvect, cols[classnames[[l]][j]]) } - legend("topright",paste(gsub("(^.+)\\..*$","\\1",basename(files[c(classnames[[k]])]))), col = colvect, lty = lty, pch = pch) - }#end length ==1 - dev.off() + legend("topright", paste(gsub("(^.+)\\..*$", "\\1", basename(files[c(classnames[[k]], classnames[[l]])]))), col = colvect, lty = lty, pch = pch) + } + } + } # end if length >2 + + if (length(class) == 2) { + k <- 1 + l <- 2 + plot(0, 0, type = "n", xlim = xlim / 60, ylim = ylim, main = paste("Total Ion Chromatograms \n", "TICs_", class[k], "vs", class[l], sep = ""), xlab = "Retention Time (min)", ylab = "TIC") + colvect <- NULL + for (j in 1:length(classnames[[k]])) { + tic <- TIC[[classnames[[k]][j]]] + # points(tic[,1]/60, tic[,2], col = cols[i], pch = pch[i], type="l") + points(tic[, 1] / 60, tic[, 2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type = "l") + colvect <- append(colvect, cols[classnames[[k]][j]]) + } + for (j in 1:length(classnames[[l]])) { + # i=class2names[j] + tic <- TIC[[classnames[[l]][j]]] + points(tic[, 1] / 60, -tic[, 2], col = cols[classnames[[l]][j]], pch = pch[classnames[[l]][j]], type = "l") + colvect <- append(colvect, cols[classnames[[l]][j]]) + } + legend("topright", paste(gsub("(^.+)\\..*$", "\\1", basename(files[c(classnames[[k]], classnames[[l]])]))), col = colvect, lty = lty, pch = pch) + } # end length ==2 + + if (length(class) == 1) { + k <- 1 + ylim <- range(sapply(TIC, function(x) range(x[, 2]))) + plot(0, 0, type = "n", xlim = xlim / 60, ylim = ylim, main = paste("Total Ion Chromatograms \n", "TICs_", class[k], sep = ""), xlab = "Retention Time (min)", ylab = "TIC") + colvect <- NULL + for (j in 1:length(classnames[[k]])) { + tic <- TIC[[classnames[[k]][j]]] + # points(tic[,1]/60, tic[,2], col = cols[i], pch = pch[i], type="l") + points(tic[, 1] / 60, tic[, 2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type = "l") + colvect <- append(colvect, cols[classnames[[k]][j]]) + } + legend("topright", paste(gsub("(^.+)\\..*$", "\\1", basename(files[c(classnames[[k]])]))), col = colvect, lty = lty, pch = pch) + } # end length ==1 + dev.off() } -#Update by J.Saint-Vanne -##addition for quality control of peak picking -#metaMS EIC and pspectra plotting option -#version 20190520 -#only for Galaxy -plotUnknowns<-function(resGC, unkn="", DB=NULL, fileFrom=NULL){ - - ##Annotation table each value is a pcgrp associated to the unknown - ##NOTE pcgrp index are different between xcmsSet and resGC due to filtering steps in metaMS - ##R. Wehrens give me some clues on that and we found a correction - - #correction of annotation matrix due to pcgrp removal by quality check in runGCresult - #matrix of correspondance between an@pspectra and filtered pspectra from runGC - #Select only pspectra which correpond to them select in metaMS - # col1 = filtered spectra from runGC and col2 = an@spectra - allPCGRPs <-lapply(1:length(resGC$xset), - function(i) { - an <- resGC$xset[[i]] - huhn <- an@pspectra[which(sapply(an@pspectra, length) >= - metaSetting(resGC$settings, - "DBconstruction.minfeat"))] - matCORR<-cbind(1:length(huhn), match(huhn, an@pspectra)) - }) - #Build a new annotation list with sampnames and pseudospectra number from xset - helpannotation <- list() - for(j in 1:length(resGC$xset)){ - helpannotation[[j]] <- resGC$annotation[[j]][1:2] - pspvector <- vector() - for(i in 1: nrow(helpannotation[[j]])){ - #Find corresponding pspec - psplink <- allPCGRPs[[j]][match(helpannotation[[j]][i,1],allPCGRPs[[j]]),2] - pspvector <- c(pspvector,psplink) - #Change the annotation column into sampname column - if(helpannotation[[j]][i,2] < 0){ - #It's an unknown - new_name <- paste("Unknown",abs(as.integer(helpannotation[[j]][i,2]))) - helpannotation[[j]][i,2] <- new_name - }else{ - #It has been found in local database - for(k in 1:length(DB)){ - if(helpannotation[[j]][i,2] == k){ - helpannotation[[j]][i,2] <- DB[[k]]$Name - break - } - } +# Update by J.Saint-Vanne +## addition for quality control of peak picking +# metaMS EIC and pspectra plotting option +# version 20190520 +# only for Galaxy +plotUnknowns <- function(resGC, unkn = "", DB = NULL, fileFrom = NULL) { + ## Annotation table each value is a pcgrp associated to the unknown + ## NOTE pcgrp index are different between xcmsSet and resGC due to filtering steps in metaMS + ## R. Wehrens give me some clues on that and we found a correction + + # correction of annotation matrix due to pcgrp removal by quality check in runGCresult + # matrix of correspondance between an@pspectra and filtered pspectra from runGC + # Select only pspectra which correpond to them select in metaMS + # col1 = filtered spectra from runGC and col2 = an@spectra + allPCGRPs <- lapply( + 1:length(resGC$xset), + function(i) { + an <- resGC$xset[[i]] + huhn <- an@pspectra[which(sapply(an@pspectra, length) >= + metaSetting( + resGC$settings, + "DBconstruction.minfeat" + ))] + matCORR <- cbind(1:length(huhn), match(huhn, an@pspectra)) + } + ) + # Build a new annotation list with sampnames and pseudospectra number from xset + helpannotation <- list() + for (j in 1:length(resGC$xset)) { + helpannotation[[j]] <- resGC$annotation[[j]][1:2] + pspvector <- vector() + for (i in 1:nrow(helpannotation[[j]])) { + # Find corresponding pspec + psplink <- allPCGRPs[[j]][match(helpannotation[[j]][i, 1], allPCGRPs[[j]]), 2] + pspvector <- c(pspvector, psplink) + # Change the annotation column into sampname column + if (helpannotation[[j]][i, 2] < 0) { + # It's an unknown + new_name <- paste("Unknown", abs(as.integer(helpannotation[[j]][i, 2]))) + helpannotation[[j]][i, 2] <- new_name + } else { + # It has been found in local database + for (k in 1:length(DB)) { + if (helpannotation[[j]][i, 2] == k) { + helpannotation[[j]][i, 2] <- DB[[k]]$Name + break + } + } + } + } + helpannotation[[j]] <- cbind(helpannotation[[j]], pspvector) + names(helpannotation)[j] <- names(resGC$annotation[j]) + } + peaktable <- resGC$PeakTable + + par(mar = c(5, 4, 4, 2) + 0.1) + # For each unknown + for (l in 1:length(unkn)) { + # recordPlot + perpage <- 3 # if change change layout also! + dev.new(width = 21 / 2.54, height = 29.7 / 2.54, file = paste("Unknown_", unkn[l], ".pdf", sep = "")) # A4 pdf + # par(mfrow=c(perpage,2)) + layout(matrix(c(1, 1, 2, 3, 4, 4, 5, 6, 7, 7, 8, 9), 6, 2, byrow = TRUE), widths = rep(c(1, 1), perpage), heights = rep(c(1, 5), perpage)) + # layout.show(6) + oma.saved <- par("oma") + par(oma = rep.int(0, 4)) + par(oma = oma.saved) + o.par <- par(mar = rep.int(0, 4)) + on.exit(par(o.par)) + # For each sample + for (c in 1:length(resGC$xset)) { + # get sample name + sampname <- basename(resGC$xset[[c]]@xcmsSet@filepaths) + # remove .cdf, .mzXML filepattern + filepattern <- c( + "[Cc][Dd][Ff]", "[Nn][Cc]", "([Mm][Zz])?[Xx][Mm][Ll]", + "[Mm][Zz][Dd][Aa][Tt][Aa]", "[Mm][Zz][Mm][Ll]" + ) + filepattern <- paste(paste("\\.", filepattern, "$", sep = ""), collapse = "|") + sampname <- gsub(filepattern, "", sampname) + title1 <- paste(peaktable[unkn[l], 1], "from", sampname, sep = " ") + an <- resGC$xset[[c]] + if (fileFrom == "zipfile") { + an@xcmsSet@filepaths <- paste0("./", an@xcmsSet@phenoData[, "class"], "/", basename(an@xcmsSet@filepaths)) + } # else { + # print(an@xcmsSet@filepaths) + # an@xcmsSet@filepaths <- paste0("./",basename(an@xcmsSet@filepaths)) + # } + # Find the good annotation for this sample + for (a in 1:length(helpannotation)) { + if (gsub(filepattern, "", names(helpannotation)[a]) == paste0("./", sampname)) { + # Find the unkn or the matched std in this sample + findunkn <- FALSE + for (r in 1:nrow(helpannotation[[a]])) { + if (helpannotation[[a]][r, "annotation"] == peaktable[unkn[l], 1]) { + findunkn <- TRUE + pcgrp <- helpannotation[[a]][r, "pspvector"] + par(mar = c(0, 0, 0, 0) + 0.1) + # Write title + plot.new() + box() + text(0.5, 0.5, title1, cex = 2) + par(mar = c(3, 2.5, 3, 1.5) + 0.1) + # Window for EIC + plotEICs(an, pspec = pcgrp, maxlabel = 2) + # Window for pseudospectra + plotPsSpectrum(an, pspec = pcgrp, maxlabel = 2) } + } + if (!findunkn) { + par(mar = c(0, 0, 0, 0) + 0.1) + # Write title + plot.new() + box() + text(0.5, 0.5, title1, cex = 2) + # Window for EIC + plot.new() + box() + text(0.5, 0.5, "NOT FOUND", cex = 2) + # Window for pseudospectra + plot.new() + box() + text(0.5, 0.5, "NOT FOUND", cex = 2) + } + break } - helpannotation[[j]] <- cbind(helpannotation[[j]],pspvector) - names(helpannotation)[j] <- names(resGC$annotation[j]) + } } - peaktable <- resGC$PeakTable - - par (mar=c(5, 4, 4, 2) + 0.1) - #For each unknown - for (l in 1:length(unkn)){ - #recordPlot - perpage=3 #if change change layout also! - dev.new(width=21/2.54, height=29.7/2.54, file=paste("Unknown_",unkn[l],".pdf", sep="")) #A4 pdf - # par(mfrow=c(perpage,2)) - layout(matrix(c(1,1,2,3,4,4,5,6,7,7,8,9), 6, 2, byrow = TRUE), widths=rep(c(1,1),perpage), heights=rep(c(1,5),perpage)) - # layout.show(6) - oma.saved <- par("oma") - par(oma = rep.int(0, 4)) - par(oma = oma.saved) - o.par <- par(mar = rep.int(0, 4)) - on.exit(par(o.par)) - #For each sample - for (c in 1:length(resGC$xset)) { - #get sample name - sampname<-basename(resGC$xset[[c]]@xcmsSet@filepaths) - #remove .cdf, .mzXML filepattern - filepattern <- c("[Cc][Dd][Ff]", "[Nn][Cc]", "([Mm][Zz])?[Xx][Mm][Ll]", - "[Mm][Zz][Dd][Aa][Tt][Aa]", "[Mm][Zz][Mm][Ll]") - filepattern <- paste(paste("\\.", filepattern, "$", sep = ""), collapse = "|") - sampname<-gsub(filepattern, "",sampname) - title1<-paste(peaktable[unkn[l],1],"from",sampname, sep = " ") - an<-resGC$xset[[c]] - if(fileFrom == "zipfile") { - an@xcmsSet@filepaths <- paste0("./",an@xcmsSet@phenoData[,"class"],"/",basename(an@xcmsSet@filepaths)) - }#else { - #print(an@xcmsSet@filepaths) - #an@xcmsSet@filepaths <- paste0("./",basename(an@xcmsSet@filepaths)) - #} - #Find the good annotation for this sample - for(a in 1:length(helpannotation)){ - if(gsub(filepattern, "", names(helpannotation)[a]) == paste0("./",sampname)){ - #Find the unkn or the matched std in this sample - findunkn <- FALSE - for(r in 1:nrow(helpannotation[[a]])){ - if(helpannotation[[a]][r,"annotation"] == peaktable[unkn[l],1]){ - findunkn <- TRUE - pcgrp <- helpannotation[[a]][r,"pspvector"] - par (mar=c(0, 0, 0, 0) + 0.1) - #Write title - plot.new() - box() - text(0.5, 0.5, title1, cex=2) - par (mar=c(3, 2.5, 3, 1.5) + 0.1) - #Window for EIC - plotEICs(an, pspec=pcgrp, maxlabel=2) - #Window for pseudospectra - plotPsSpectrum(an, pspec=pcgrp, maxlabel=2) - } - } - if(!findunkn) { - par (mar=c(0, 0, 0, 0) + 0.1) - #Write title - plot.new() - box() - text(0.5, 0.5, title1, cex=2) - #Window for EIC - plot.new() - box() - text(0.5, 0.5, "NOT FOUND", cex=2) - #Window for pseudospectra - plot.new() - box() - text(0.5, 0.5, "NOT FOUND", cex=2) - } - break - } - } - } - graphics.off() - }#end for unkn[l] -}#end function + graphics.off() + } # end for unkn[l] +} # end function