Mercurial > repos > marie-tremblay-metatoul > nmr_alignment
comparison nmr_alignement/NmrAlignment_wrapper.R @ 1:58eecef626da draft
Uploaded
author | marie-tremblay-metatoul |
---|---|
date | Thu, 02 Mar 2017 10:23:37 -0500 |
parents | |
children | 908e1345d7ca |
comparison
equal
deleted
inserted
replaced
0:d690c5ad932f | 1:58eecef626da |
---|---|
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()) |