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
|
4
|
74 ## Exclusion zone(s)
|
1
|
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
|
4
|
94 ## Graphical zone(s)
|
|
95 graphicalZones <- argLs[["zone_graphical_choices.choice"]]
|
|
96 graphicalZonesBorders <- NULL
|
|
97 if (!is.null(argLs$zone_exclusion_left))
|
|
98 {
|
|
99 for(i in which(names(argLs)=="zone_graphical_left"))
|
|
100 {
|
|
101 graphicalZonesBorders <- c(graphicalZonesBorders,list(c(argLs[[i]],argLs[[i+1]])))
|
|
102 }
|
|
103 }
|
|
104
|
1
|
105
|
|
106 # Outputs
|
|
107 logOut <- argLs[["logOut"]]
|
|
108 alignedSpectra <- argLs[["alignedSpectra"]]
|
|
109 graphOut <- argLs[["graphOut"]]
|
|
110
|
|
111
|
4
|
112 ## Checking R packages
|
|
113 ##--------------------
|
|
114 sink(logOut)
|
|
115 cat("\tPACKAGE INFO\n")
|
|
116 #pkgs=c("xcms","batch")
|
|
117 pkgs=c("batch","speaq")
|
|
118 for(pkg in pkgs) {
|
|
119 suppressPackageStartupMessages( stopifnot( library(pkg, quietly=TRUE, logical.return=TRUE, character.only=TRUE)))
|
|
120 cat(pkg,"\t",as.character(packageVersion(pkg)),"\n",sep="")
|
|
121 }
|
|
122
|
1
|
123 ## Checking arguments
|
|
124 ##-------------------
|
|
125 error.stock <- "\n"
|
|
126 if(length(error.stock) > 1)
|
|
127 stop(error.stock)
|
3
|
128
|
|
129
|
1
|
130 ## Computation
|
|
131 ##------------
|
2
|
132 directory.alignement <- nmr.alignment(fileType=fileType,directory=directory,leftBorder=leftBorder,rightBorder=rightBorder,exclusionZones=exclusionZones,
|
3
|
133 exclusionZonesBorders=exclusionZonesBorders, reference=reference, nDivRange=nDivRange,
|
1
|
134 baselineThresh=baselineThresh, maxshift=50, verbose=FALSE)
|
|
135 directory.raw <- directory.alignement[[1]]
|
|
136 directory.aligned <- directory.alignement[[2]]
|
|
137
|
|
138 ## Saving
|
|
139 ##-------
|
|
140 # Aligned spectra
|
|
141 t.directory.aligned <- t(directory.aligned)
|
|
142 rownames(t.directory.aligned) <- colnames(directory.aligned)
|
|
143 # colnames(t.directory.aligned) <- c("Bucket",colnames(t.directory.aligned))
|
3
|
144 write.table(t.directory.aligned,file=alignedSpectra,row.names=TRUE,quote=FALSE,sep="\t")
|
1
|
145
|
|
146
|
|
147 ## Graphical output: overlay of raw and estimated spectra
|
|
148 pdf(graphOut,onefile=TRUE)
|
4
|
149
|
|
150 graphical.zone.length <- length(graphicalZonesBorders)
|
1
|
151 par(mfrow=c(2,1))
|
|
152
|
|
153 raw.spectra <- data.frame(directory.raw)
|
|
154 colnames(raw.spectra) <- substr(colnames(raw.spectra),2,7)
|
|
155
|
|
156 aligned.spectra <- data.frame(directory.aligned)
|
|
157 colnames(aligned.spectra) <- substr(colnames(aligned.spectra),2,7)
|
|
158
|
|
159 drawSpec(raw.spectra,xlab="", ylab="Raw spectra", main="")
|
|
160 drawSpec(aligned.spectra,xlab="", ylab="Aligned spectra", main="")
|
|
161
|
4
|
162 if (graphical.zone.length != 0)
|
1
|
163 {
|
4
|
164 # par(mfrow=c(2,1))
|
|
165 for (g in 1:graphical.zone.length)
|
1
|
166 {
|
4
|
167 print(paste(g, graphicalZonesBorders[[g]][1], graphicalZonesBorders[[g]][2]))
|
|
168 drawSpec(raw.spectra[,(which(round(as.numeric(colnames(raw.spectra)),2) == graphicalZonesBorders[[g]][1])[1]):(which(round(as.numeric(colnames(raw.spectra)),2) == graphicalZonesBorders[[g]][2])[1])],xlab="", main="")
|
|
169 drawSpec(aligned.spectra[,(which(round(as.numeric(colnames(aligned.spectra)),2) == graphicalZonesBorders[[g]][1])[1]):(which(round(as.numeric(colnames(aligned.spectra)),2) == graphicalZonesBorders[[g]][2])[1])],xlab="", main="")
|
1
|
170 }
|
|
171 }
|
|
172 dev.off()
|
|
173
|
|
174
|
|
175 ## Ending
|
|
176 ##---------------------
|
|
177 cat("\nEnd of 'NMR alignment' Galaxy module call: ", as.character(Sys.time()), sep = "")
|
4
|
178 sink()
|
1
|
179 options(stringsAsFactors = strAsFacL)
|
|
180 rm(list = ls())
|