Mercurial > repos > pieterlukasse > prims_metabolomics

annotate xcms_get_mass_eic.r @ 50:93102202ab79

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