1
|
1 #!/usr/local/public/bin/Rscript --vanilla --slave --no-site-file
|
|
2
|
|
3 ## 06102016_NmrAlignment_wrapper.R
|
|
4 ## Marie Tremblay-Franco
|
|
5 ## marie.tremblay-franco@toulouse.inra.fr
|
|
6
|
|
7 runExampleL <- FALSE
|
|
8
|
|
9 ##------------------------------
|
|
10 ## Options
|
|
11 ##------------------------------
|
|
12 strAsFacL <- options()$stringsAsFactors
|
|
13 options(stringsAsFactors = FALSE)
|
|
14
|
|
15
|
|
16 ##------------------------------
|
|
17 ## Libraries loading
|
|
18 ##------------------------------
|
|
19 # ParseCommandArgs function
|
|
20 library(batch)
|
|
21 # Alignment
|
|
22 library(speaq)
|
|
23
|
|
24
|
|
25 # R script call
|
|
26 source_local <- function(fname)
|
|
27 {
|
|
28 argv <- commandArgs(trailingOnly = FALSE)
|
|
29 base_dir <- dirname(substring(argv[grep("--file=", argv)], 8))
|
|
30 source(paste(base_dir, fname, sep="/"))
|
|
31 }
|
|
32 # Function import
|
|
33 source_local("NmrAlignment_script.R")
|
|
34
|
|
35
|
|
36 ##------------------------------
|
|
37 ## Errors ?????????????????????
|
|
38 ##------------------------------
|
|
39
|
|
40
|
|
41 ##------------------------------
|
|
42 ## Constants
|
|
43 ##------------------------------
|
|
44 topEnvC <- environment()
|
|
45 flagC <- "\n"
|
|
46
|
|
47
|
|
48 ##------------------------------
|
|
49 ## Script
|
|
50 ##------------------------------
|
|
51 if(!runExampleL)
|
|
52 argLs <- parseCommandArgs(evaluate=FALSE)
|
|
53
|
|
54
|
|
55 ## Parameters Loading
|
|
56 ##-------------------
|
|
57 # Inputs
|
|
58 ## Library of spectra to align
|
|
59 if (!is.null(argLs[["zipfile"]])){
|
|
60 zipfile= argLs[["zipfile"]]
|
|
61 directory=unzip(zipfile, list=F)
|
|
62 directory=paste(getwd(),strsplit(directory[1],"/")[[1]][2],sep="/")
|
|
63 } else if (!is.null(argLs[["library"]])){
|
|
64 directory=argLs[["library"]]
|
|
65 if(!file.exists(directory)){
|
|
66 error_message=paste("Cannot access the directory :",directory,".Please verify if the directory exists or not.")
|
|
67 print(error_message)
|
|
68 stop(error_message)
|
|
69 }
|
|
70 }
|
|
71
|
|
72
|
|
73 ## Spectral width
|
|
74 leftBorder <- argLs[["left_border"]]
|
|
75 rightBorder <- argLs[["right_border"]]
|
|
76
|
|
77 ##Exclusion zone(s)
|
|
78 exclusionZones <- argLs[["zone_exclusion_choices.choice"]]
|
|
79 exclusionZonesBorders <- NULL
|
|
80 if (!is.null(argLs$zone_exclusion_left))
|
|
81 {
|
|
82 for(i in which(names(argLs)=="zone_exclusion_left"))
|
|
83 {
|
|
84 exclusionZonesBorders <- c(exclusionZonesBorders,list(c(argLs[[i]],argLs[[i+1]])))
|
|
85 }
|
|
86 }
|
|
87
|
|
88 ## Reference spectrum
|
|
89 reference <- argLs[["reference"]]
|
|
90
|
|
91 ## Size of a small nDivRange
|
|
92 nDivRange <- argLs[["nDivRange"]]
|
|
93
|
|
94 ## Intensity threshold for peak removal
|
|
95 baselineThresh <- argLs[["baselineThresh"]]
|
|
96
|
|
97
|
|
98 # Outputs
|
|
99 logOut <- argLs[["logOut"]]
|
|
100 alignedSpectra <- argLs[["alignedSpectra"]]
|
|
101 graphOut <- argLs[["graphOut"]]
|
|
102
|
|
103
|
|
104 ## Checking arguments
|
|
105 ##-------------------
|
|
106 error.stock <- "\n"
|
|
107 if(length(error.stock) > 1)
|
|
108 stop(error.stock)
|
|
109
|
|
110
|
|
111 ## Computation
|
|
112 ##------------
|
|
113 directory.alignement <- nmr.alignment(directory=directory,leftBorder=leftBorder,rightBorder=rightBorder,exclusionZones=exclusionZones,
|
|
114 exclusionZonesBorders=exclusionZonesBorders, reference=reference, nDivRange=nDivRange,
|
|
115 baselineThresh=baselineThresh, maxshift=50, verbose=FALSE)
|
|
116 directory.raw <- directory.alignement[[1]]
|
|
117 directory.aligned <- directory.alignement[[2]]
|
|
118
|
|
119 ## Saving
|
|
120 ##-------
|
|
121 # Aligned spectra
|
|
122 t.directory.aligned <- t(directory.aligned)
|
|
123 rownames(t.directory.aligned) <- colnames(directory.aligned)
|
|
124 # colnames(t.directory.aligned) <- c("Bucket",colnames(t.directory.aligned))
|
|
125 write.table(t.directory.aligned,file=alignedSpectra,row.names=TRUE,quote=FALSE,sep="\t")
|
|
126
|
|
127
|
|
128 excludedZone <- NULL
|
|
129 for (c in 1:length(exclusionZonesBorders))
|
|
130 {
|
|
131 excludedZone <- c(excludedZone,exclusionZonesBorders[[c]])
|
|
132 excludedZone <- sort(excludedZone)
|
|
133 }
|
|
134
|
|
135 ## Graphical output: overlay of raw and estimated spectra
|
|
136 pdf(graphOut,onefile=TRUE)
|
|
137 par(mfrow=c(2,1))
|
|
138
|
|
139 raw.spectra <- data.frame(directory.raw)
|
|
140 colnames(raw.spectra) <- substr(colnames(raw.spectra),2,7)
|
|
141
|
|
142 aligned.spectra <- data.frame(directory.aligned)
|
|
143 colnames(aligned.spectra) <- substr(colnames(aligned.spectra),2,7)
|
|
144
|
|
145 drawSpec(raw.spectra,xlab="", ylab="Raw spectra", main="")
|
|
146 drawSpec(aligned.spectra,xlab="", ylab="Aligned spectra", main="")
|
|
147
|
|
148 nbZones <- length(excludedZone)/2
|
|
149 if (nbZones != 0)
|
|
150 {
|
|
151 n <- length(excludedZone)
|
|
152 drawSpec(raw.spectra[,1:which(round(as.numeric(colnames(raw.spectra)),2) == excludedZone[n])[1]],xlab="", ylab="Raw spectra", main="")
|
|
153 drawSpec(aligned.spectra[,1:which(round(as.numeric(colnames(aligned.spectra)),2) == excludedZone[n])[1]],xlab="", ylab="Aligned spectra", main="")
|
|
154
|
|
155 n <- n - 1
|
|
156 while (n >= nbZones & nbZones > 1)
|
|
157 {
|
|
158 drawSpec(raw.spectra[,(which(round(as.numeric(colnames(raw.spectra)),2) == excludedZone[n])[1]):(which(round(as.numeric(colnames(raw.spectra)),2) == excludedZone[n-1])[1])],xlab="", ylab="Raw spectra", main="")
|
|
159 drawSpec(aligned.spectra[,(which(round(as.numeric(colnames(aligned.spectra)),2) == excludedZone[n])[1]):(which(round(as.numeric(colnames(aligned.spectra)),2) == excludedZone[n-1])[1])],xlab="", ylab="Aligned spectra", main="")
|
|
160 n <- n - 2
|
|
161 }
|
|
162
|
|
163 drawSpec(raw.spectra[,(which(round(as.numeric(colnames(raw.spectra)),2) == excludedZone[1])[1]):ncol(raw.spectra)],xlab="", ylab="Raw spectra", main="")
|
|
164 drawSpec(aligned.spectra[,(which(round(as.numeric(colnames(aligned.spectra)),2) == excludedZone[1])[1]):ncol(aligned.spectra)],xlab="", ylab="Aligned spectra", main="")
|
|
165 }
|
|
166 drawSpec(raw.spectra[,(which(round(as.numeric(colnames(raw.spectra)),2) == 2.4)[1]):(which(round(as.numeric(colnames(raw.spectra)),2) == 2.8)[1])],xlab="", ylab="Raw spectra", main="")
|
|
167 drawSpec(aligned.spectra[,(which(round(as.numeric(colnames(aligned.spectra)),2) == 2.4)[1]):(which(round(as.numeric(colnames(aligned.spectra)),2) == 2.8)[1])],xlab="", ylab="Aligned spectra", main="")
|
|
168
|
|
169 dev.off()
|
|
170
|
|
171
|
|
172 ## Ending
|
|
173 ##---------------------
|
|
174 cat("\nEnd of 'NMR alignment' Galaxy module call: ", as.character(Sys.time()), sep = "")
|
|
175 options(stringsAsFactors = strAsFacL)
|
|
176 rm(list = ls())
|