Mercurial > repos > pieterlukasse > prims_metabolomics
view 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 source
## 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) } )