Mercurial > repos > pieterlukasse > prims_metabolomics

annotate xcms_get_alignment_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 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 )