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
|
3
|
20 library(batch)
|
1
|
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
|
3
|
54
|
1
|
55 ## Parameters Loading
|
|
56 ##-------------------
|
|
57 # Inputs
|
|
58 ## Library of spectra to align
|
|
59 if (!is.null(argLs[["zipfile"]])){
|
2
|
60 fileType="zip"
|
1
|
61 zipfile= argLs[["zipfile"]]
|
|
62 directory=unzip(zipfile, list=F)
|
|
63 directory=paste(getwd(),strsplit(directory[1],"/")[[1]][2],sep="/")
|
2
|
64 } else if (!is.null(argLs[["tsvfile"]])){
|
|
65 fileType="tsv"
|
|
66 directory <- read.table(argLs[["tsvfile"]],check.names=FALSE,header=TRUE,sep="\t")
|
1
|
67 }
|
|
68
|
|
69
|
|
70 ## Spectral width
|
|
71 leftBorder <- argLs[["left_border"]]
|
|
72 rightBorder <- argLs[["right_border"]]
|
3
|
73
|
1
|
74 ##Exclusion zone(s)
|
|
75 exclusionZones <- argLs[["zone_exclusion_choices.choice"]]
|
|
76 exclusionZonesBorders <- NULL
|
|
77 if (!is.null(argLs$zone_exclusion_left))
|
|
78 {
|
|
79 for(i in which(names(argLs)=="zone_exclusion_left"))
|
|
80 {
|
|
81 exclusionZonesBorders <- c(exclusionZonesBorders,list(c(argLs[[i]],argLs[[i+1]])))
|
|
82 }
|
|
83 }
|
3
|
84
|
1
|
85 ## Reference spectrum
|
|
86 reference <- argLs[["reference"]]
|
|
87
|
|
88 ## Size of a small nDivRange
|
|
89 nDivRange <- argLs[["nDivRange"]]
|
3
|
90
|
1
|
91 ## Intensity threshold for peak removal
|
|
92 baselineThresh <- argLs[["baselineThresh"]]
|
|
93
|
|
94
|
|
95 # Outputs
|
|
96 logOut <- argLs[["logOut"]]
|
|
97 alignedSpectra <- argLs[["alignedSpectra"]]
|
|
98 graphOut <- argLs[["graphOut"]]
|
|
99
|
|
100
|
|
101 ## Checking arguments
|
|
102 ##-------------------
|
|
103 error.stock <- "\n"
|
|
104 if(length(error.stock) > 1)
|
|
105 stop(error.stock)
|
3
|
106
|
|
107
|
1
|
108 ## Computation
|
|
109 ##------------
|
2
|
110 directory.alignement <- nmr.alignment(fileType=fileType,directory=directory,leftBorder=leftBorder,rightBorder=rightBorder,exclusionZones=exclusionZones,
|
3
|
111 exclusionZonesBorders=exclusionZonesBorders, reference=reference, nDivRange=nDivRange,
|
1
|
112 baselineThresh=baselineThresh, maxshift=50, verbose=FALSE)
|
|
113 directory.raw <- directory.alignement[[1]]
|
|
114 directory.aligned <- directory.alignement[[2]]
|
|
115
|
|
116 ## Saving
|
|
117 ##-------
|
|
118 # Aligned spectra
|
|
119 t.directory.aligned <- t(directory.aligned)
|
|
120 rownames(t.directory.aligned) <- colnames(directory.aligned)
|
|
121 # colnames(t.directory.aligned) <- c("Bucket",colnames(t.directory.aligned))
|
3
|
122 write.table(t.directory.aligned,file=alignedSpectra,row.names=TRUE,quote=FALSE,sep="\t")
|
1
|
123
|
|
124
|
|
125 excludedZone <- NULL
|
|
126 for (c in 1:length(exclusionZonesBorders))
|
|
127 {
|
|
128 excludedZone <- c(excludedZone,exclusionZonesBorders[[c]])
|
|
129 excludedZone <- sort(excludedZone)
|
|
130 }
|
|
131
|
|
132 ## Graphical output: overlay of raw and estimated spectra
|
|
133 pdf(graphOut,onefile=TRUE)
|
|
134 par(mfrow=c(2,1))
|
|
135
|
|
136 raw.spectra <- data.frame(directory.raw)
|
|
137 colnames(raw.spectra) <- substr(colnames(raw.spectra),2,7)
|
|
138
|
|
139 aligned.spectra <- data.frame(directory.aligned)
|
|
140 colnames(aligned.spectra) <- substr(colnames(aligned.spectra),2,7)
|
|
141
|
|
142 drawSpec(raw.spectra,xlab="", ylab="Raw spectra", main="")
|
|
143 drawSpec(aligned.spectra,xlab="", ylab="Aligned spectra", main="")
|
|
144
|
|
145 nbZones <- length(excludedZone)/2
|
|
146 if (nbZones != 0)
|
|
147 {
|
|
148 n <- length(excludedZone)
|
|
149 drawSpec(raw.spectra[,1:which(round(as.numeric(colnames(raw.spectra)),2) == excludedZone[n])[1]],xlab="", ylab="Raw spectra", main="")
|
|
150 drawSpec(aligned.spectra[,1:which(round(as.numeric(colnames(aligned.spectra)),2) == excludedZone[n])[1]],xlab="", ylab="Aligned spectra", main="")
|
|
151
|
|
152 n <- n - 1
|
|
153 while (n >= nbZones & nbZones > 1)
|
|
154 {
|
|
155 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="")
|
|
156 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="")
|
|
157 n <- n - 2
|
|
158 }
|
3
|
159
|
1
|
160 drawSpec(raw.spectra[,(which(round(as.numeric(colnames(raw.spectra)),2) == excludedZone[1])[1]):ncol(raw.spectra)],xlab="", ylab="Raw spectra", main="")
|
|
161 drawSpec(aligned.spectra[,(which(round(as.numeric(colnames(aligned.spectra)),2) == excludedZone[1])[1]):ncol(aligned.spectra)],xlab="", ylab="Aligned spectra", main="")
|
|
162 }
|
|
163 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="")
|
|
164 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="")
|
|
165
|
|
166 dev.off()
|
|
167
|
|
168
|
|
169 ## Ending
|
|
170 ##---------------------
|
|
171 cat("\nEnd of 'NMR alignment' Galaxy module call: ", as.character(Sys.time()), sep = "")
|
|
172 options(stringsAsFactors = strAsFacL)
|
|
173 rm(list = ls())
|