|
2
|
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())
|
