Mercurial > repos > pieterlukasse > prims_metabolomics

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