Mercurial > repos > pieterlukasse > prims_metabolomics

comparison xcms_get_mass_eic.r @ 49:f772a5caa86a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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
comparison
equal deleted inserted replaced
48:26b93438f30e 49:f772a5caa86a
1 ## read args:
2 args <- commandArgs(TRUE)
3 # xset data:
4 args.xsetData <- args[1]
5
6 args.rtStart <- strtoi(args[2])
7 args.rtEnd <- strtoi(args[3])
8
9 args.mzStart <- as.double(args[4])
10 args.mzEnd <- as.double(args[5])
11 # there are 2 options: specify a mz range or a mz list:
12 if (args.mzStart < 0)
13 {
14 args.mzList <- as.double(strsplit(args[6], ",")[[1]])
15 cat(typeof(as.double(strsplit(args[6], ",")[[1]])))
16 args.mzTolPpm <- as.double(args[7])
17 # calculate mzends based on ppm tol:
18 mzListEnd <- c()
19 mzListStart <- c()
20 for (i in 1:length(args.mzList))
21 {
22 mzEnd <- args.mzList[i] + args.mzList[i]*args.mzTolPpm/1000000.0
23 mzStart <- args.mzList[i] - args.mzList[i]*args.mzTolPpm/1000000.0
24 mzListEnd <- c(mzListEnd, mzEnd)
25 mzListStart <- c(mzListStart, mzStart)
26 }
27 str(mzListStart)
28 str(mzListEnd)
29 } else {
30 mzListEnd <- c(args.mzEnd)
31 mzListStart <- c(args.mzStart)
32 }
33
34 args.sampleNames <- strsplit(args[8], ",")[[1]]
35 # trim leading and trailing spaces:
36 args.sampleNames <- gsub("^\\s+|\\s+$", "", args.sampleNames)
37
38 args.combineSamples <- args[9]
39 args.rtPlotMode <- args[10]
40
41 ## report files
42 args.htmlReportFile <- args[11]
43 args.htmlReportFile.files_path <- args[12]
44
45
46 if (length(args) == 13)
47 {
48 args.outLogFile <- args[13]
49 # suppress messages:
50 # Send all STDERR to STDOUT using sink() see http://mazamascience.com/WorkingWithData/?p=888
51 msg <- file(args.outLogFile, open="wt")
52 sink(msg, type="message")
53 sink(msg, type="output")
54 }
55
56 # TODO - add option to do masses in same plot (if given in same line oid) or in separate plots
57 # TODO2 - let it run in parallel
58
59 tryCatch(
60 {
61 library(metaMS)
62
63 # load the xset data :
64 xsetData <- readRDS(args.xsetData)
65 # if here to support both scenarios:
66 if ("xcmsSet" %in% slotNames(xsetData) )
67 {
68 xsetData <- xsetData@xcmsSet
69 }
70
71 # report
72 dir.create(file.path(args.htmlReportFile.files_path), showWarnings = FALSE, recursive = TRUE)
73 message(paste("\nGenerating report.........in ", args.htmlReportFile.files_path))
74
75 html <- "<html><body><h1>Extracted Ion Chromatograms (EIC) matching criteria</h1>"
76
77 if (args.combineSamples == "No")
78 {
79 if (length(args.sampleNames) > 1 && length(mzListStart) > 1 && length(args.sampleNames) != length(mzListStart))
80 stop(paste("The number of sample names should match the number of m/z values in the list. Found ", length(mzListStart),
81 " masses while ", length(args.sampleNames), " sample names were given."))
82
83 iterSize <- length(args.sampleNames)
84 # 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
85 fixSampleIdx <- 1
86 fixMzListIdx <- 1
87 if (length(args.sampleNames) == 1)
88 {
89 fixSampleIdx <- 0
90 iterSize <- length(mzListStart)
91 }
92 if (length(mzListStart) == 1)
93 {
94 fixMzListIdx <- 0
95 }
96 lineColors <- rainbow(iterSize)
97 for (i in 0:(iterSize-1))
98 {
99 message("\nGetting EIC... ")
100 eiccor <- getEIC(xsetData,
101 mzrange=matrix(c(mzListStart[i*fixMzListIdx+1],mzListEnd[i*fixMzListIdx+1]),nrow=1,ncol=2,byrow=TRUE),
102 rtrange=matrix(c(args.rtStart,args.rtEnd),nrow=1,ncol=2,byrow=TRUE),
103 sampleidx=c(args.sampleNames[i*fixSampleIdx+1]), rt=args.rtPlotMode)
104
105 message("\nPlotting figures... ")
106 figureName <- paste(args.htmlReportFile.files_path, "/figure", i,".png", sep="")
107 html <- paste(html,"<img src='", "figure", i,".png' /><br/>", sep="")
108 png( figureName, type="cairo", width=1100,height=250 )
109 #plot(eiccor, col=lineColors[i+1])
110 # black is better in this case:
111 plot(eiccor)
112 legend('topright', # places a legend at the appropriate place
113 legend=c(args.sampleNames[i*fixSampleIdx+1]), # puts text in the legend
114 lty=c(1,1), # gives the legend appropriate symbols (lines)
115 lwd=c(2.5,2.5))
116
117 devname = dev.off()
118 }
119
120 } else {
121 for (i in 1:length(mzListStart))
122 {
123 message("\nGetting EIC... ")
124 eiccor <- getEIC(xsetData,
125 mzrange=matrix(c(mzListStart[i],mzListEnd[i]),nrow=1,ncol=2,byrow=TRUE),
126 rtrange=matrix(c(args.rtStart,args.rtEnd),nrow=1,ncol=2,byrow=TRUE),
127 sampleidx=args.sampleNames, rt = args.rtPlotMode)
128
129 #set size, set option (plot per sample, plot per mass)
130
131 message("\nPlotting figures... ")
132 figureName <- paste(args.htmlReportFile.files_path, "/figure", i,".png", sep="")
133 html <- paste(html,"<img src='", "figure", i,".png' />", sep="")
134 png( figureName, type="cairo", width=1100,height=450 )
135 lineColors <- rainbow(length(args.sampleNames))
136 plot(eiccor, col=lineColors)
137 legend('topright', # places a legend at the appropriate place
138 legend=args.sampleNames, # puts text in the legend
139 lty=c(1,1), # gives the legend appropriate symbols (lines)
140 lwd=c(2.5,2.5),
141 col=lineColors)
142 devname = dev.off()
143 }
144 }
145 if (args.rtPlotMode == "corrected")
146 {
147 html <- paste(html,"<p>*rt values are corrected ones</p></body><html>")
148 }
149 html <- paste(html,"</body><html>")
150 message("finished generating report")
151 write(html,file=args.htmlReportFile)
152 # unlink(args.htmlReportFile)
153 cat("\nWarnings================:\n")
154 str( warnings() )
155 },
156 error=function(cond) {
157 sink(NULL, type="message") # default setting
158 sink(stderr(), type="output")
159 message("\nERROR: ===========\n")
160 print(cond)
161 }
162 )