|
6
|
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", width=800 )
|
|
|
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 )
|
