Mercurial > repos > pieterlukasse > prims_metabolomics
diff xcms_get_mass_eic.r @ 49:f772a5caa86a
Added more options and better documentation.
Added MsClust support for parsing XCMS alignment results.
Improved output reports for XCMS wrappers.
New tools.
| author | pieter.lukasse@wur.nl |
|---|---|
| date | Wed, 10 Dec 2014 22:03:27 +0100 |
| parents | |
| children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/xcms_get_mass_eic.r Wed Dec 10 22:03:27 2014 +0100 @@ -0,0 +1,162 @@ +## read args: +args <- commandArgs(TRUE) +# xset data: +args.xsetData <- args[1] + +args.rtStart <- strtoi(args[2]) +args.rtEnd <- strtoi(args[3]) + +args.mzStart <- as.double(args[4]) +args.mzEnd <- as.double(args[5]) +# there are 2 options: specify a mz range or a mz list: +if (args.mzStart < 0) +{ + args.mzList <- as.double(strsplit(args[6], ",")[[1]]) + cat(typeof(as.double(strsplit(args[6], ",")[[1]]))) + args.mzTolPpm <- as.double(args[7]) + # calculate mzends based on ppm tol: + mzListEnd <- c() + mzListStart <- c() + for (i in 1:length(args.mzList)) + { + mzEnd <- args.mzList[i] + args.mzList[i]*args.mzTolPpm/1000000.0 + mzStart <- args.mzList[i] - args.mzList[i]*args.mzTolPpm/1000000.0 + mzListEnd <- c(mzListEnd, mzEnd) + mzListStart <- c(mzListStart, mzStart) + } + str(mzListStart) + str(mzListEnd) +} else { + mzListEnd <- c(args.mzEnd) + mzListStart <- c(args.mzStart) +} + +args.sampleNames <- strsplit(args[8], ",")[[1]] +# trim leading and trailing spaces: +args.sampleNames <- gsub("^\\s+|\\s+$", "", args.sampleNames) + +args.combineSamples <- args[9] +args.rtPlotMode <- args[10] + +## report files +args.htmlReportFile <- args[11] +args.htmlReportFile.files_path <- args[12] + + +if (length(args) == 13) +{ + args.outLogFile <- args[13] + # suppress messages: + # Send all STDERR to STDOUT using sink() see http://mazamascience.com/WorkingWithData/?p=888 + msg <- file(args.outLogFile, open="wt") + sink(msg, type="message") + sink(msg, type="output") +} + +# TODO - add option to do masses in same plot (if given in same line oid) or in separate plots +# TODO2 - let it run in parallel + +tryCatch( + { + library(metaMS) + + # load the xset data : + xsetData <- readRDS(args.xsetData) + # if here to support both scenarios: + if ("xcmsSet" %in% slotNames(xsetData) ) + { + xsetData <- xsetData@xcmsSet + } + + # report + dir.create(file.path(args.htmlReportFile.files_path), showWarnings = FALSE, recursive = TRUE) + message(paste("\nGenerating report.........in ", args.htmlReportFile.files_path)) + + html <- "<html><body><h1>Extracted Ion Chromatograms (EIC) matching criteria</h1>" + + if (args.combineSamples == "No") + { + if (length(args.sampleNames) > 1 && length(mzListStart) > 1 && length(args.sampleNames) != length(mzListStart)) + stop(paste("The number of sample names should match the number of m/z values in the list. Found ", length(mzListStart), + " masses while ", length(args.sampleNames), " sample names were given.")) + + iterSize <- length(args.sampleNames) + # these can be set to 1 or 0 just as a trick to iterate OR not over the items. If the respective list is of length 1, only the first item should be used + fixSampleIdx <- 1 + fixMzListIdx <- 1 + if (length(args.sampleNames) == 1) + { + fixSampleIdx <- 0 + iterSize <- length(mzListStart) + } + if (length(mzListStart) == 1) + { + fixMzListIdx <- 0 + } + lineColors <- rainbow(iterSize) + for (i in 0:(iterSize-1)) + { + message("\nGetting EIC... ") + eiccor <- getEIC(xsetData, + mzrange=matrix(c(mzListStart[i*fixMzListIdx+1],mzListEnd[i*fixMzListIdx+1]),nrow=1,ncol=2,byrow=TRUE), + rtrange=matrix(c(args.rtStart,args.rtEnd),nrow=1,ncol=2,byrow=TRUE), + sampleidx=c(args.sampleNames[i*fixSampleIdx+1]), rt=args.rtPlotMode) + + message("\nPlotting figures... ") + figureName <- paste(args.htmlReportFile.files_path, "/figure", i,".png", sep="") + html <- paste(html,"<img src='", "figure", i,".png' /><br/>", sep="") + png( figureName, type="cairo", width=1100,height=250 ) + #plot(eiccor, col=lineColors[i+1]) + # black is better in this case: + plot(eiccor) + legend('topright', # places a legend at the appropriate place + legend=c(args.sampleNames[i*fixSampleIdx+1]), # puts text in the legend + lty=c(1,1), # gives the legend appropriate symbols (lines) + lwd=c(2.5,2.5)) + + devname = dev.off() + } + + } else { + for (i in 1:length(mzListStart)) + { + message("\nGetting EIC... ") + eiccor <- getEIC(xsetData, + mzrange=matrix(c(mzListStart[i],mzListEnd[i]),nrow=1,ncol=2,byrow=TRUE), + rtrange=matrix(c(args.rtStart,args.rtEnd),nrow=1,ncol=2,byrow=TRUE), + sampleidx=args.sampleNames, rt = args.rtPlotMode) + + #set size, set option (plot per sample, plot per mass) + + message("\nPlotting figures... ") + figureName <- paste(args.htmlReportFile.files_path, "/figure", i,".png", sep="") + html <- paste(html,"<img src='", "figure", i,".png' />", sep="") + png( figureName, type="cairo", width=1100,height=450 ) + lineColors <- rainbow(length(args.sampleNames)) + plot(eiccor, col=lineColors) + legend('topright', # places a legend at the appropriate place + legend=args.sampleNames, # puts text in the legend + lty=c(1,1), # gives the legend appropriate symbols (lines) + lwd=c(2.5,2.5), + col=lineColors) + devname = dev.off() + } + } + if (args.rtPlotMode == "corrected") + { + html <- paste(html,"<p>*rt values are corrected ones</p></body><html>") + } + html <- paste(html,"</body><html>") + message("finished generating report") + write(html,file=args.htmlReportFile) + # unlink(args.htmlReportFile) + cat("\nWarnings================:\n") + str( warnings() ) + }, + error=function(cond) { + sink(NULL, type="message") # default setting + sink(stderr(), type="output") + message("\nERROR: ===========\n") + print(cond) + } + )