|
2
|
1 #!/usr/local/public/bin/Rscript --vanilla --slave --no-site-file
|
|
|
2
|
|
|
3 ## 070115_NmrBucketing2galaxy_v1.R
|
|
|
4 ## Marie Tremblay-Franco
|
|
|
5 ## MetaboHUB: The French Infrastructure for Metabolomics and Fluxomics
|
|
|
6 ## www.metabohub.fr/en
|
|
|
7 ## marie.tremblay-franco@toulouse.inra.fr
|
|
|
8
|
|
|
9 runExampleL <- FALSE
|
|
|
10
|
|
|
11 if(runExampleL) {
|
|
|
12 ##------------------------------
|
|
|
13 ## Example of arguments
|
|
|
14 ##------------------------------
|
|
|
15 argLs <- list(StudyDir = "Tlse_BPASourisCerveau",
|
|
|
16 upper = "10.0",
|
|
|
17 lower = "0.50",
|
|
|
18 bucket.width = "0.01",
|
|
|
19 exclusion = "TRUE",
|
|
|
20 exclusion.zone = list(c(6.5,4.5)),
|
|
|
21 graph="Overlay")
|
|
|
22
|
|
|
23 argLs <- c(argLs,
|
|
|
24 list(dataMatrixOut = paste(directory,"_NmrBucketing_dataMatrix.tsv",sep=""),
|
|
|
25 sampleMetadataOut = paste(directory,"_NmrBucketing_sampleMetadata.tsv",sep=""),
|
|
|
26 variableMetadataOut = paste(directory,"_NmrBucketing_variableMetadata.tsv",sep=""),
|
|
|
27 graphOut = paste(directory,"_NmrBucketing_graph.pdf",sep=""),
|
|
|
28 logOut = paste(directory,"_NmrBucketing_log.txt",sep="")))
|
|
|
29 }
|
|
|
30
|
|
|
31 ##------------------------------
|
|
|
32 ## Options
|
|
|
33 ##------------------------------
|
|
|
34 strAsFacL <- options()$stringsAsFactors
|
|
|
35 options(stringsAsFactors = FALSE)
|
|
|
36
|
|
|
37
|
|
|
38 ##------------------------------
|
|
|
39 ## Libraries laoding
|
|
|
40 ##------------------------------
|
|
|
41 # For parseCommandArgs function
|
|
|
42 library(batch)
|
|
|
43 # For cumtrapz function
|
|
|
44 library(pracma)
|
|
|
45
|
|
|
46 # R script call
|
|
|
47 source_local <- function(fname)
|
|
|
48 {
|
|
|
49 argv <- commandArgs(trailingOnly = FALSE)
|
|
|
50 base_dir <- dirname(substring(argv[grep("--file=", argv)], 8))
|
|
|
51 source(paste(base_dir, fname, sep="/"))
|
|
|
52 }
|
|
|
53 #Import the different functions
|
|
|
54 source_local("NmrBucketing_script.R")
|
|
|
55 source_local("DrawSpec.R")
|
|
|
56
|
|
|
57 ##------------------------------
|
|
|
58 ## Errors ?????????????????????
|
|
|
59 ##------------------------------
|
|
|
60
|
|
|
61
|
|
|
62 ##------------------------------
|
|
|
63 ## Constants
|
|
|
64 ##------------------------------
|
|
|
65 topEnvC <- environment()
|
|
|
66 flagC <- "\n"
|
|
|
67
|
|
|
68
|
|
|
69 ##------------------------------
|
|
|
70 ## Script
|
|
|
71 ##------------------------------
|
|
|
72 if(!runExampleL)
|
|
|
73 argLs <- parseCommandArgs(evaluate=FALSE)
|
|
|
74
|
|
|
75
|
|
|
76 ## Parameters Loading
|
|
|
77 ##-------------------
|
|
|
78 # Inputs
|
|
|
79 if (!is.null(argLs[["zipfile"]])){
|
|
|
80 fileType="zip"
|
|
|
81 zipfile= argLs[["zipfile"]]
|
|
|
82 directory=unzip(zipfile, list=F)
|
|
|
83 directory=paste(getwd(),strsplit(directory[1],"/")[[1]][2],sep="/")
|
|
|
84 } else if (!is.null(argLs[["tsvfile"]])){
|
|
|
85 fileType="tsv"
|
|
|
86 directory <- read.table(argLs[["tsvfile"]],check.names=FALSE,header=TRUE,sep="\t")
|
|
|
87 }
|
|
|
88
|
|
|
89 leftBorder <- argLs[["left_border"]]
|
|
|
90 rightBorder <- argLs[["right_border"]]
|
|
|
91 bucketSize <- argLs[["bucket_width"]]
|
|
|
92 exclusionZones <- argLs[["zone_exclusion_choices.choice"]]
|
|
|
93
|
|
|
94 exclusionZonesBorders <- NULL
|
|
|
95 if (!is.null(argLs$zone_exclusion_left))
|
|
|
96 {
|
|
|
97 for(i in which(names(argLs)=="zone_exclusion_left"))
|
|
|
98 {
|
|
|
99 exclusionZonesBorders <- c(exclusionZonesBorders,list(c(argLs[[i]],argLs[[i+1]])))
|
|
|
100 }
|
|
|
101 }
|
|
|
102
|
|
|
103 graphique <- argLs[["graphType"]]
|
|
|
104
|
|
|
105 # Outputs
|
|
|
106 nomGraphe <- argLs[["graphOut"]]
|
|
|
107 dataMatrixOut <- argLs[["dataMatrixOut"]]
|
|
|
108 logFile <- argLs[["logOut"]]
|
|
|
109 if (fileType=="zip")
|
|
|
110 {
|
|
|
111 sampleMetadataOut <- argLs[["sampleOut"]]
|
|
|
112 variableMetadataOut <- argLs[["variableOut"]]
|
|
|
113 }
|
|
|
114
|
|
|
115 ## Checking arguments
|
|
|
116 ##-------------------
|
|
|
117 error.stock <- "\n"
|
|
|
118
|
|
|
119 if(length(error.stock) > 1)
|
|
|
120 stop(error.stock)
|
|
|
121
|
|
|
122
|
|
|
123 ## Computation
|
|
|
124 ##------------
|
|
|
125 outputs <- NmrBucketing(fileType=fileType, fileName=directory, leftBorder=leftBorder, rightBorder=rightBorder, bucketSize=bucketSize,
|
|
|
126 exclusionZones=exclusionZones, exclusionZonesBorders=exclusionZonesBorders, graph=graphique, nomFichier=nomGraphe,
|
|
|
127 savLog.txtC=logFile)
|
|
|
128 data_bucket <- outputs[[1]]
|
|
|
129 data_sample <- outputs[[2]]
|
|
|
130 data_variable <- outputs[[3]]
|
|
|
131 ppm <- outputs[[4]]
|
|
|
132 ppm <- round(ppm,2)
|
|
|
133
|
|
|
134 ## Graphical outputs
|
|
|
135 ##------------------
|
|
|
136 if (graphique != "None")
|
|
|
137 {
|
|
|
138 excludedZone <- NULL
|
|
|
139 for (c in 1:length(exclusionZonesBorders))
|
|
|
140 {
|
|
|
141 excludedZone <- c(excludedZone,exclusionZonesBorders[[c]])
|
|
|
142 excludedZone <- sort(excludedZone)
|
|
|
143 }
|
|
|
144 nbZones <- length(excludedZone)/2
|
|
|
145 n <- length(excludedZone)
|
|
|
146
|
|
|
147 # Graphic Device opening
|
|
|
148 pdf(nomGraphe,onefile=TRUE)
|
|
|
149
|
|
|
150 if (graphique == "Overlay")
|
|
|
151 {
|
|
|
152 # Global spectral window
|
|
|
153 spectra <- data.frame(t(data_bucket))
|
|
|
154 drawSpec(spectra,xlab="", ylab="Intensity", main="")
|
|
|
155
|
|
|
156 ## Zoomed spectral window depending on exclusion zone(s)
|
|
|
157 if (nbZones != 0)
|
|
|
158 {
|
|
|
159 BInf <- excludedZone[n]
|
|
|
160 if (round(BInf,1) == BInf)
|
|
|
161 {
|
|
|
162 BInf <- BInf+0.01
|
|
|
163 }
|
|
|
164 spectra <- data.frame(t(data_bucket[1:(which(ppm == BInf)[[1]]),]))
|
|
|
165 drawSpec(spectra,xlab="", ylab="Intensity", main="")
|
|
|
166 n <- n - 1
|
|
|
167
|
|
|
168 while (n >= nbZones & nbZones > 1)
|
|
|
169 {
|
|
|
170 BInf <- excludedZone[n-1]
|
|
|
171 if (round(BInf,1) > BInf)
|
|
|
172 {
|
|
|
173 BInf <- BInf+0.01
|
|
|
174 }
|
|
|
175 spectra <- data.frame(t(data_bucket[(which(ppm == excludedZone[n])[[1]]):(which(ppm == BInf)[[1]]),]))
|
|
|
176 drawSpec(spectra,xlab="", ylab="Intensity", main="")
|
|
|
177 n <- n - 2
|
|
|
178 }
|
|
|
179
|
|
|
180 BInf <- excludedZone[1]
|
|
|
181 if (round(BInf,1) <= BInf)
|
|
|
182 {
|
|
|
183 BInf <- BInf+0.01
|
|
|
184 }
|
|
|
185 spectra <- data.frame(t(data_bucket[(which(ppm == BInf)[[1]]):nrow(data_bucket),]))
|
|
|
186 drawSpec(spectra,xlab="", ylab="Intensity", main="")
|
|
|
187 }
|
|
|
188 }
|
|
|
189 else
|
|
|
190 {
|
|
|
191 for (i in 1:ncol(data_bucket))
|
|
|
192 {
|
|
|
193 par(mfrow=c((nbZones+2),1))
|
|
|
194 n <- length(excludedZone)
|
|
|
195 spectra <- t(data_bucket[,i])
|
|
|
196 names(spectra) <- rownames(data_bucket)
|
|
|
197 plot(1:length(spectra), spectra, type='l', xlab="", ylab="Intensity", main=colnames(data_bucket)[i], xaxt = "n")
|
|
|
198 xPos <- 1
|
|
|
199 nAxisPos <- 4
|
|
|
200 startP <- length(nAxisPos)
|
|
|
201 endP <- nrow(data_bucket)
|
|
|
202 GraphRange <- c(startP:endP)
|
|
|
203 tempVal = trunc(length(GraphRange)/nAxisPos)
|
|
|
204 xPos = c(0:nAxisPos) * tempVal
|
|
|
205 axis(1, at = xPos, labels = rownames(data_bucket)[xPos + startP])
|
|
|
206
|
|
|
207 ## Zoomed spectral window depending on exclusion zone(s)
|
|
|
208 if (nbZones != 0)
|
|
|
209 {
|
|
|
210 BInf <- excludedZone[n]
|
|
|
211 if (round(BInf,1) == BInf)
|
|
|
212 {
|
|
|
213 BInf <- BInf+0.01
|
|
|
214 }
|
|
|
215 spectra <- t(data_bucket[1:(which(ppm == BInf)[[1]]),i])
|
|
|
216 names(spectra) <- rownames(data_bucket)[1:(which(ppm == BInf)[[1]])]
|
|
|
217 plot(1:length(spectra), spectra, type='l',xlab="", ylab="Intensity", main="", xaxt = "n")
|
|
|
218 xPos <- 1
|
|
|
219 nAxisPos <- 4
|
|
|
220 startP <- length(nAxisPos)
|
|
|
221 endP <- length(spectra)
|
|
|
222 GraphRange <- c(startP:endP)
|
|
|
223 tempVal = trunc(length(GraphRange)/nAxisPos)
|
|
|
224 xPos = c(0:nAxisPos) * tempVal
|
|
|
225 axis(1, at = xPos, labels = rownames(data_bucket)[xPos + startP])
|
|
|
226 n <- n - 1
|
|
|
227
|
|
|
228 while (n >= nbZones & nbZones > 1)
|
|
|
229 {
|
|
|
230 BInf <- excludedZone[n-1]
|
|
|
231 if (round(BInf,1) > BInf)
|
|
|
232 {
|
|
|
233 BInf <- BInf+0.01
|
|
|
234 }
|
|
|
235 spectra <- t(data_bucket[(which(ppm == excludedZone[n])[[1]]):(which(ppm == BInf)[[1]]),i])
|
|
|
236 names(spectra) <- rownames(data_bucket)[(which(ppm == excludedZone[n])[[1]]):(which(ppm == BInf)[[1]])]
|
|
|
237 plot(1:length(spectra), spectra, type='l',xlab="", ylab="Intensity", main="", xaxt = "n")
|
|
|
238 xPos <- 1
|
|
|
239 nAxisPos <- 4
|
|
|
240 startP <- length(nAxisPos)
|
|
|
241 endP <- length(spectra)
|
|
|
242 GraphRange <- c(startP:endP)
|
|
|
243 tempVal = trunc(length(GraphRange)/nAxisPos)
|
|
|
244 xPos = c(0:nAxisPos) * tempVal
|
|
|
245 axis(1, at = xPos, labels = rownames(data_bucket)[xPos + startP])
|
|
|
246 n <- n - 2
|
|
|
247 }
|
|
|
248
|
|
|
249 BInf <- excludedZone[1]
|
|
|
250 if (round(BInf,1) <= BInf)
|
|
|
251 {
|
|
|
252 BInf <- BInf+0.01
|
|
|
253 }
|
|
|
254 spectra <- t(data_bucket[(which(ppm == BInf)[[1]]):nrow(data_bucket),i])
|
|
|
255 names(spectra) <- rownames(data_bucket)[(which(ppm == BInf)[[1]]):nrow(data_bucket)]
|
|
|
256 plot(1:length(spectra), spectra, type='l',xlab="", ylab="Intensity", main="", xaxt = "n")
|
|
|
257 xPos <- 1
|
|
|
258 nAxisPos <- 4
|
|
|
259 startP <- length(nAxisPos)
|
|
|
260 endP <- length(spectra)
|
|
|
261 GraphRange <- c(startP:endP)
|
|
|
262 tempVal = trunc(length(GraphRange)/nAxisPos)
|
|
|
263 xPos = c(0:nAxisPos) * tempVal
|
|
|
264 axis(1, at = xPos, labels = rownames(data_bucket)[xPos + startP])
|
|
|
265 }
|
|
|
266 }
|
|
|
267 }
|
|
|
268 dev.off()
|
|
|
269 }
|
|
|
270 ## Saving
|
|
|
271 ##-------
|
|
|
272 # Data
|
|
|
273 data_bucket <- cbind(rownames(data_bucket),data_bucket)
|
|
|
274 colnames(data_bucket) <- c("Bucket",colnames(data_bucket)[-1])
|
|
|
275 write.table(data_bucket,file=argLs$dataMatrixOut,quote=FALSE,row.names=FALSE,sep="\t")
|
|
|
276 # Sample
|
|
|
277 data_sample <- cbind(rownames(data_sample),data_sample)
|
|
|
278 colnames(data_sample) <- c("Sample",colnames(data_sample)[-1])
|
|
|
279 write.table(data_sample,file=argLs$sampleOut,quote=FALSE,row.names=FALSE,sep="\t")
|
|
|
280 # Variable
|
|
|
281 data_variable <- cbind(rownames(data_variable),data_variable)
|
|
|
282 colnames(data_variable) <- c("Bucket",colnames(data_variable)[-1])
|
|
|
283 write.table(data_variable,file=argLs$variableOut,quote=FALSE,row.names=FALSE,sep="\t")
|
|
|
284
|
|
|
285
|
|
|
286 ## Ending
|
|
|
287 ##---------------------
|
|
|
288
|
|
|
289 cat("\nEnd of 'NMR bucketing' Galaxy module call: ", as.character(Sys.time()), sep = "")
|
|
|
290
|
|
|
291 ## sink(NULL)
|
|
|
292
|
|
|
293 options(stringsAsFactors = strAsFacL)
|
|
|
294
|
|
|
295 rm(list = ls())
|
