Mercurial > repos > marie-tremblay-metatoul > nmr_annotation

comparison nmr_preprocessing/NmrPreprocessing_wrapper.R @ 2:7304ec2c9ab7 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author marie-tremblay-metatoul
date Mon, 30 Jul 2018 10:33:03 -0400
parents
children
comparison
equal deleted inserted replaced
1:b55559a2854f 2:7304ec2c9ab7
1 #!/usr/local/public/bin/Rscript --vanilla --slave --no-site-file
2
3 ## 170116_NmrPreprocessing.R
4 ## Manon Martin and Marie Tremblay-Franco
5
6 ##======================================================
7 ##======================================================
8 # Preamble
9 ##======================================================
10 ##======================================================
11
12 runExampleL <- FALSE
13
14
15 ##------------------------------
16 ## Options
17 ##------------------------------
18 strAsFacL <- options()$stringsAsFactors
19 options(stringsAsFactors = FALSE)
20
21 ##------------------------------
22 ## Libraries laoding
23 ##------------------------------
24 library(batch)
25 library(ptw)
26 library(Matrix)
27 library(ggplot2)
28 library(gridExtra)
29 library(reshape2)
30
31
32 # R script call
33 source_local <- function(fname)
34 {
35 argv <- commandArgs(trailingOnly = FALSE)
36 base_dir <- dirname(substring(argv[grep("--file=", argv)], 8))
37 source(paste(base_dir, fname, sep="/"))
38 }
39 #Import the different functions
40 source_local("NmrPreprocessing_script.R")
41 source_local("DrawFunctions.R")
42
43 ##------------------------------
44 ## Script
45 ##------------------------------
46 runExampleL <- FALSE
47
48
49 if(!runExampleL)
50 argLs <- parseCommandArgs(evaluate=FALSE)
51
52 sink(argLs$logOut)
53
54
55 ##------------------------------
56 ## Errors ?????????????????????
57 ##------------------------------
58
59
60 ##------------------------------
61 ## Constants
62 ##------------------------------
63 topEnvC <- environment()
64 flagC <- "\n"
65
66
67
68
69 # log file
70 # print(argLs[["logOut"]])
71
72 ## Starting
73 cat("\nStart of 'Preprocessing' Galaxy module call: ", as.character(Sys.time()), "\n", sep = "")
74
75
76 ##======================================================
77 ##======================================================
78 ## Parameters Loading
79 ##======================================================
80 ##======================================================
81
82 # graphical inputs
83 FirstOPCGraph <- argLs[["FirstOPCGraph"]]
84 SSGraph <- argLs[["SSGraph"]]
85 ApodGraph <- argLs[["ApodGraph"]]
86 FTGraph <- argLs[["FTGraph"]]
87 SRGraph <- argLs[["SRGraph"]]
88 ZeroOPCGraph <- argLs[["ZeroOPCGraph"]]
89 BCGraph <- argLs[["BCGraph"]]
90 FinalGraph <- argLs[["FinalGraph"]]
91
92
93 # 1rst order phase correction ------------------------
94 # Inputs
95 ## Data matrix
96 Fid_data0 <- read.table(argLs[["dataMatrixFid"]],header=TRUE, check.names=FALSE, sep='\t')
97 # Fid_data0 <- Fid_data0[,-1]
98 Fid_data0 <- as.matrix(Fid_data0)
99
100 ## Samplemetadata
101 samplemetadataFid <- read.table(argLs[["sampleMetadataFid"]],check.names=FALSE,header=TRUE,sep="\t")
102 samplemetadataFid <- as.matrix(samplemetadataFid)
103
104
105 # water and solvent(s) correction ------------------------
106 # Inputs
107 lambda <- argLs[["lambda"]]
108
109
110
111 # apodization -----------------------------------------
112 # Inputs
113 phase=0
114 rectRatio=1/2
115 gaussLB=1
116 expLB=1
117 apodization <- argLs[["apodizationMethod"]]
118
119 if (apodization=='exp'){
120 expLB <- argLs[["expLB"]]
121 } else if (apodization=='cos2'){
122 phase <- argLs[["phase"]]
123 } else if (apodization=='hanning'){
124 phase <- argLs[["phase"]]
125 } else if (apodization=='hamming'){
126 phase <- argLs[["phase"]]
127 } else if (apodization=='blockexp'){
128 rectRatio <- argLs[["rectRatio"]]
129 expLB <- argLs[["expLB"]]
130 } else if (apodization=='blockcos2'){
131 rectRatio <- argLs[["rectRatio"]]
132 } else if (apodization=='gauss'){
133 rectRatio <- argLs[["rectRatio"]]
134 gaussLB <- argLs[["gaussLB"]]
135 }
136
137
138 # Fourier transform ----------------------------------
139 # Inputs
140
141
142 # Zero Order Phase Correction -------------------------------
143 # Inputs
144
145 angle = NULL
146 excludeZOPC = NULL
147
148
149 zeroOrderPhaseMethod <- argLs[["zeroOrderPhaseMethod"]]
150
151 if (zeroOrderPhaseMethod=='manual'){
152 angle <- argLs[["angle"]]
153 }
154
155 excludeZoneZeroPhase <- argLs[["excludeZoneZeroPhase.choice"]]
156 if (excludeZoneZeroPhase == 'YES') {
157 excludeZoneZeroPhaseList <- list()
158 for(i in which(names(argLs)=="excludeZoneZeroPhase_left")) {
159 excludeZoneZeroPhaseLeft <- argLs[[i]]
160 excludeZoneZeroPhaseRight <- argLs[[i+1]]
161 excludeZoneZeroPhaseList <- c(excludeZoneZeroPhaseList,list(c(excludeZoneZeroPhaseLeft,excludeZoneZeroPhaseRight)))
162 }
163 excludeZOPC <- excludeZoneZeroPhaseList
164 }
165
166
167 # Internal referencering ----------------------------------
168 # Inputs
169 shiftTreshold = 2 # c
170 ppm = TRUE
171 shiftReferencingRangeList = NULL # fromto.RC
172 pctNearValue = 0.02 # pc
173 rowindex_graph = NULL
174 ppm_ref = 0 # ppm.ref
175
176 #
177 # shiftReferencing <- argLs[["shiftReferencing"]]
178 # print(shiftReferencing)
179 #
180 # if (shiftReferencing=="YES")
181 # {
182 #
183 # shiftReferencingMethod <- argLs[["shiftReferencingMethod"]]
184 #
185 # if (shiftReferencingMethod == "thres") {
186 # shiftTreshold <- argLs[["shiftTreshold"]]
187 # }
188
189 shiftReferencingRange <- argLs[["shiftReferencingRange"]]
190
191 if (shiftReferencingRange == "near0"){
192 pctNearValue <- argLs[["pctNearValue"]]
193 }
194
195 if (shiftReferencingRange == "window"){
196 shiftReferencingRangeList <- list()
197 for(i in which(names(argLs)=="shiftReferencingRangeLeft"))
198 {
199 shiftReferencingRangeLeft <- argLs[[i]]
200 shiftReferencingRangeRight <- argLs[[i+1]]
201 shiftReferencingRangeList <- c(shiftReferencingRangeList,list(c(shiftReferencingRangeLeft,shiftReferencingRangeRight)))
202 }
203 }
204
205 shiftHandling <- argLs[["shiftHandling"]]
206
207 ppmvalue <- argLs[["ppmvalue"]]
208
209
210
211 # }
212
213
214 # Baseline Correction -------------------------------
215 # Inputs
216 lambdaBc <- argLs[["lambdaBc"]]
217 pBc <- argLs[["pBc"]]
218 epsilon <- argLs[["epsilon"]]
219
220 excludeBC = NULL
221
222 excludeZoneBC <- argLs[["excludeZoneBC.choice"]]
223 if (excludeZoneBC == 'YES') {
224 excludeZoneBCList <- list()
225 for(i in which(names(argLs)=="excludeZoneBC_left")) {
226 excludeZoneBCLeft <- argLs[[i]]
227 excludeZoneBCRight <- argLs[[i+1]]
228 excludeZoneBCList <- c(excludeZoneBCList,list(c(excludeZoneBCLeft,excludeZoneBCRight)))
229 }
230 excludeBC <- excludeZoneBCList
231 }
232
233 # transformation of negative values -------------------------------
234 # Inputs
235 NegativetoZero <- argLs[["NegativetoZero"]]
236
237
238 # Outputs
239 nomGraphe <- argLs[["graphOut"]]
240 # dataMatrixOut <- argLs[["dataMatrixOut"]]
241 log <- argLs[["logOut"]]
242
243
244
245 ## Checking arguments
246 ##-------------------
247 error.stock <- "\n"
248
249 if(length(error.stock) > 1)
250 stop(error.stock)
251
252
253 ##======================================================
254 ##======================================================
255 ## Computation
256 ##======================================================
257 ##======================================================
258
259 pdf(nomGraphe, onefile = TRUE, width = 13, height = 13)
260
261 # FirstOrderPhaseCorrection ---------------------------------
262 Fid_data <- GroupDelayCorrection(Fid_data0, Fid_info = samplemetadataFid, group_delay = NULL)
263
264 if (FirstOPCGraph == "YES") {
265 title = "FIDs after Group Delay Correction"
266 DrawSignal(Fid_data, subtype = "stacked",
267 ReImModArg = c(TRUE, FALSE, FALSE, FALSE), vertical = T,
268 xlab = "Frequency", num.stacked = 4,
269 main = title, createWindow=FALSE)
270 }
271
272 # SolventSuppression ---------------------------------
273 Fid_data <- SolventSuppression(Fid_data, lambda.ss = lambda, ptw.ss = TRUE, plotSolvent = F, returnSolvent = F)
274
275 if (SSGraph == "YES") {
276 title = "FIDs after Solvent Suppression "
277 DrawSignal(Fid_data, subtype = "stacked",
278 ReImModArg = c(TRUE, FALSE, FALSE, FALSE), vertical = T,
279 xlab = "Frequency", num.stacked = 4,
280 main = title, createWindow=FALSE)
281 }
282
283
284 # Apodization ---------------------------------
285 Fid_data <- Apodization(Fid_data, Fid_info = samplemetadataFid, DT = NULL,
286 type.apod = apodization, phase = phase, rectRatio = rectRatio, gaussLB = gaussLB, expLB = expLB, plotWindow = F, returnFactor = F)
287
288 if (ApodGraph == "YES") {
289 title = "FIDs after Apodization"
290 DrawSignal(Fid_data, subtype = "stacked",
291 ReImModArg = c(TRUE, FALSE, FALSE, FALSE), vertical = T,
292 xlab = "Frequency", num.stacked = 4,
293 main = title, createWindow=FALSE)
294 }
295
296
297 # FourierTransform ---------------------------------
298 Spectrum_data <- FourierTransform(Fid_data, Fid_info = samplemetadataFid, reverse.axis = TRUE)
299
300
301 if (FTGraph == "YES") {
302 title = "Fourier transformed spectra"
303 DrawSignal(Spectrum_data, subtype = "stacked",
304 ReImModArg = c(TRUE, FALSE, FALSE, FALSE), vertical = T,
305 xlab = "Frequency", num.stacked = 4,
306 main = title, createWindow=FALSE)
307 }
308
309
310
311 # ZeroOrderPhaseCorrection ---------------------------------
312 Spectrum_data <- ZeroOrderPhaseCorrection(Spectrum_data, type.zopc = zeroOrderPhaseMethod,
313 plot_rms = NULL, returnAngle = FALSE,
314 createWindow = TRUE,angle = angle,
315 plot_spectra = FALSE,
316 ppm.zopc = TRUE, exclude.zopc = excludeZOPC)
317
318
319 # InternalReferencing ---------------------------------
320 # if (shiftReferencing=="YES") {
321 Spectrum_data <- InternalReferencing(Spectrum_data, samplemetadataFid, method = "max", range = shiftReferencingRange,
322 ppm.value = ppmvalue, shiftHandling = shiftHandling, ppm.ir = TRUE,
323 fromto.RC = shiftReferencingRangeList, pc = pctNearValue)
324
325 if (SRGraph == "YES") {
326 title = "Spectra after Shift Referencing"
327 DrawSignal(Spectrum_data, subtype = "stacked",
328 ReImModArg = c(TRUE, FALSE, FALSE, FALSE), vertical = T,
329 xlab = "Frequency", num.stacked = 4,
330 main = title, createWindow=FALSE)
331 }
332
333 # }
334
335 if (ZeroOPCGraph == "YES") {
336 title = "Spectra after Zero Order Phase Correction"
337 DrawSignal(Spectrum_data, subtype = "stacked",
338 ReImModArg = c(TRUE, FALSE, FALSE, FALSE), vertical = T,
339 xlab = "Frequency", num.stacked = 4,
340 main = title, createWindow=FALSE)
341 }
342
343
344 # BaselineCorrection ---------------------------------
345 Spectrum_data <- BaselineCorrection(Spectrum_data, ptw.bc = TRUE, lambda.bc = lambdaBc,
346 p.bc = pBc, eps = epsilon, ppm.bc = TRUE,
347 exclude.bc = excludeBC,
348 returnBaseline = F)
349
350
351
352 if (BCGraph == "YES") {
353 title = "Spectra after Baseline Correction"
354 DrawSignal(Spectrum_data, subtype = "stacked",
355 ReImModArg = c(TRUE, FALSE, FALSE, FALSE), vertical = T,
356 xlab = "Frequency", num.stacked = 4,
357 main = title, createWindow=FALSE)
358 }
359
360
361 # NegativeValuesZeroing ---------------------------------
362 if (NegativetoZero=="YES") {
363 Spectrum_data <- NegativeValuesZeroing(Spectrum_data)
364 }
365
366 if (FinalGraph == "YES") {
367 title = "Final preprocessed spectra"
368 DrawSignal(Spectrum_data, subtype = "stacked",
369 ReImModArg = c(TRUE, FALSE, FALSE, FALSE), vertical = T,
370 xlab = "Frequency", num.stacked = 4,
371 main = title, createWindow=FALSE)
372 }
373
374 invisible(dev.off())
375
376
377 data_variable <- matrix(NA, nrow = 1, ncol = dim(Spectrum_data)[2], dimnames = list("ID", NULL))
378 colnames(data_variable) <- colnames(Spectrum_data)
379 data_variable[1,] <- colnames(data_variable)
380
381
382 ##======================================================
383 ##======================================================
384 ## Saving
385 ##======================================================
386 ##======================================================
387
388 # Data Matrix
389 write.table(t(Re(Spectrum_data)),file=argLs$dataMatrix, quote=FALSE, row.names=TRUE, sep="\t", col.names=TRUE)
390
391 # Variable metadata
392 write.table(data_variable,file=argLs$variableMetadata, quote=FALSE, row.names=TRUE, sep="\t", col.names=TRUE)
393
394 # log file
395 # write.table(t(data.frame(argLs)), file = argLs$logOut, col.names = FALSE, quote=FALSE)
396
397 # input arguments
398 cat("\n INPUT and OUTPUT ARGUMENTS :\n")
399
400 argLs
401
402
403 ## Ending
404
405 cat("\nEnd of 'Preprocessing' Galaxy module call: ", as.character(Sys.time()), sep = "")
406
407 sink()
408
409 options(stringsAsFactors = strAsFacL)
410
411 rm(list = ls())