Mercurial > repos > pieterlukasse > prims_metabolomics

annotate xcms_get_alignment_eic.r @ 58:a79599cd2a9c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
fix for xcms support in msclust
author pieter.lukasse@wur.nl
date Fri, 12 Dec 2014 12:07:30 +0100
parents f772a5caa86a
children 35f506f30ae4
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 # shows all alignment results in a rt region
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
2 ## read args:
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
3 args <- commandArgs(TRUE)
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
4 # xset data:
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
5 args.xsetData <- args[1]
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
6
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
7 args.rtStart <- strtoi(args[2])
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
8 args.rtEnd <- strtoi(args[3])
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
9
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
10 # limit max diff to 600 and minNrSamples to at least 2 :
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
11 if (args.rtEnd - args.rtStart > 600)
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
12 stop("The RT window should be <= 600")
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.minNrSamples <- strtoi(args[4]) #min nr samples
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
15 if (args.minNrSamples < 2)
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
16 stop("Set 'Minimum number of samples' to 2 or higher")
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
17
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
18
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
19 args.sampleNames <- strsplit(args[5], ",")[[1]]
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
20 # trim leading and trailing spaces:
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
21 args.sampleNames <- gsub("^\\s+|\\s+$", "", args.sampleNames)
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
22
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
23 ## report files
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
24 args.htmlReportFile <- args[6]
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
25 args.htmlReportFile.files_path <- args[7]
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
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
28 if (length(args) == 8)
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
29 {
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
30 args.outLogFile <- args[8]
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
31 # suppress messages:
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
32 # 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
33 msg <- file(args.outLogFile, open="wt")
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
34 sink(msg, type="message")
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
35 sink(msg, type="output")
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
36 }
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
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
39
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
40 tryCatch(
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
41 {
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
42 library(metaMS)
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
43
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
44 # load the xset data :
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
45 xsetData <- readRDS(args.xsetData)
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
46 # if here to support both scenarios:
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
47 if ("xcmsSet" %in% slotNames(xsetData) )
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
48 {
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
49 xsetData <- xsetData@xcmsSet
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
50 }
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
51
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
52 # report
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
53 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
54 message(paste("\nGenerating report.........in ", args.htmlReportFile.files_path))
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 write(paste("<html><body><h1>Extracted Ion Chromatograms (EIC) of alignments with peaks in ",args.minNrSamples, " or more samples</h1>"),
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
57 file=args.htmlReportFile)
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 gt <- groups(xsetData)
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
60 message("\nParsed groups... ")
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
61 groupidx1 <- which(gt[,"rtmed"] > args.rtStart & gt[,"rtmed"] < args.rtEnd & gt[,"npeaks"] >= args.minNrSamples) # this should be only on samples selected
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 if (length(groupidx1) > 0)
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
64 {
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
65 message("\nGetting EIC... ")
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
66 eiccor <- getEIC(xsetData, groupidx = c(groupidx1), sampleidx=args.sampleNames)
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
67 #eicraw <- getEIC(xsetData, groupidx = c(groupidx1), sampleidx=args.sampleNames, rt = "raw")
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
68
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
69 #sampleNamesIdx <- which(sampnames(LC$xset@xcmsSet) %in% args.sampleNames, arr.ind = TRUE)
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
70 #or (from bioconductor code for getEIC): NB: this is assuming sample indexes used in data are ordered after the order of sample names!!
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
71 sampNames <- sampnames(xsetData)
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
72 sampleNamesIdx <- match( args.sampleNames, sampNames)
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
73 message(paste("Samples: ", sampleNamesIdx))
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 #TODO html <- paste(html, "<table><tbody>")
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
76 message(paste("\nPlotting figures... "))
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
77
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
78 #get the mz list (interestingly, this [,"mz"] is relatively slow):
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
79 peakMzList <- xsetData@peaks[,"mz"]
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
80 peakSampleList <- xsetData@peaks[,"sample"]
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
81 #signal to noise list:
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
82 peakSnList <- xsetData@peaks[,"sn"]
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
83
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
84 message(paste("Total nr of peaks: ", length(peakMzList)))
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
85
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
86 for (i in 1:length(groupidx1))
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
87 {
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
88 groupMembers <- xsetData@groupidx[[groupidx1[i]]]
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
89
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
90 groupMzList <- ""
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
91 groupSampleList <- ""
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
92 finalGroupSize <- 0
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 for (j in 1:length(groupMembers))
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 #get only the peaks from the selected samples:
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
97 memberSample <- peakSampleList[groupMembers[j]]
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
98 memberSn <- peakSnList[groupMembers[j]]
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
99 if (!is.na(memberSn) && memberSample %in% sampleNamesIdx)
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
100 {
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
101 message(paste("Group: ", groupidx1[i], " / Member sample: ", memberSample))
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
102 memberMz <- peakMzList[groupMembers[j]]
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
103
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 if (finalGroupSize == 0)
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
106 {
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
107 groupMzList <- memberMz
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
108 groupSampleList <- sampNames[memberSample]
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
109 } else {
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
110 groupMzList <- paste(groupMzList,",",memberMz, sep="")
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
111 groupSampleList <- paste(groupSampleList,",",sampNames[memberSample], sep="")
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
112 }
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
113
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
114 finalGroupSize <- finalGroupSize +1
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
115 }
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 # here we do the real check on group size and only the groups that have enough
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
118 # members with signal to noise > 0 will be plotted here:
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
119 if (finalGroupSize >= args.minNrSamples)
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
120 {
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
121 message(paste("Plotting figure ",i, " of max ", length(groupidx1)," figures... "))
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 figureName <- paste(args.htmlReportFile.files_path, "/figure", i,".png", sep="")
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
124 write(paste("<img src='", "figure", i,".png' />", sep="") ,
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
125 file=args.htmlReportFile, append=TRUE)
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
126
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
127 png( figureName, type="cairo" )
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
128 plot(eiccor, xsetData, groupidx = i)
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
129 devname = dev.off()
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 write(paste("<p>Alignment id: [", groupidx1[i], "]. m/z list of peaks in alignment with signal/noise <> NA:<br/>", groupMzList,"</p>", sep="") ,
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
132 file=args.htmlReportFile, append=TRUE)
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
133 write(paste("<p>m/z values found in the following samples respectively: <br/>", groupSampleList,"</p>", sep="") ,
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
134 file=args.htmlReportFile, append=TRUE)
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
135 }
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
136 }
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
137
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
138 }
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
139
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
140 write("</body><html>",
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
141 file=args.htmlReportFile, append=TRUE)
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
142 message("finished generating report")
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
143 # unlink(args.htmlReportFile)
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
144 cat("\nWarnings================:\n")
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
145 str( warnings() )
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 error=function(cond) {
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
148 sink(NULL, type="message") # default setting
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
149 sink(stderr(), type="output")
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
150 message("\nERROR: ===========\n")
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
151 print(cond)
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
152 }
f772a5caa86a Added more options and better documentation.
pieter.lukasse@wur.nl
parents:
diff changeset
153 )