Mercurial > repos > marie-tremblay-metatoul > normalization

comparison nmr_bucketing/NmrBucketing_wrapper.R @ 2:3cd762aac7a4 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author marie-tremblay-metatoul
date Thu, 20 Apr 2017 08:53:46 -0400
parents
children
comparison
equal deleted inserted replaced
1:f9554d5bb47e 2:3cd762aac7a4
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())