Mercurial > repos > pieterlukasse > prims_metabolomics2

comparison xcms_get_mass_eic.r @ 0:dffc38727496

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
initial commit
author pieter.lukasse@wur.nl
date Sat, 07 Feb 2015 22:02:00 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:dffc38727496
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 )