Mercurial > repos > marie-tremblay-metatoul > 2dnmrannotation

comparison annotationRmn2DWrapper.R @ 3:546c7ccd2ed4 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload for repository https://github.com/workflow4metabolomics/tools-metabolomics commit 911f4beba3dcb25c1033e8239426f8f763683523"
author workflow4metabolomics
date Fri, 04 Feb 2022 09:01:11 +0000
parents dff7bde22102
children
comparison
equal deleted inserted replaced
2:dff7bde22102 3:546c7ccd2ed4
1 #!/usr/local/public/bin/Rscript --vanilla --slave --no-site-file 1 #!/usr/local/public/bin/Rscript --vanilla --slave --no-site-file
2 2
3 ## 201919016 2DNmrAnnotation_1.0.0.R 3 ## 201919016 2DNmrAnnotation_1.0.0.R
4 ## Marie Tremblay-Franco 4 ## Marie Tremblay-Franco
5 ## MetaboHUB: The French Infrastructure for Metabolomics and Fluxomics 5 ## MetaboHUB: The French Infrastructure for Metabolomics and Fluxomics
6 ## www.metabohub.fr/en 6 ## marie.tremblay-franco@inrae.fr
7 ## marie.tremblay-franco@toulouse.inra.fr
8 7
9 runExampleL <- FALSE 8 runExampleL <- FALSE
10 9
11 if(runExampleL) { 10 if (runExampleL) {
12 ##------------------------------ 11 ##------------------------------
13 ## Example of arguments 12 ## Example of arguments
14 ##------------------------------ 13 ##------------------------------
15 } 14 }
16 15
18 ##------------------------------ 17 ##------------------------------
19 ## Options 18 ## Options
20 ##------------------------------ 19 ##------------------------------
21 strAsFacL <- options()$stringsAsFactors 20 strAsFacL <- options()$stringsAsFactors
22 options(stringsAsFactors = FALSE) 21 options(stringsAsFactors = FALSE)
22 options(digits = 8, scipen = 3)
23 23
24 ##------------------------------ 24 ##------------------------------
25 ## Constants 25 ## Constants
26 ##------------------------------ 26 ##------------------------------
27 topEnvC <- environment() 27 topEnvC <- environment()
39 library(dplyr) 39 library(dplyr)
40 library(ggplot2) 40 library(ggplot2)
41 library(openxlsx) 41 library(openxlsx)
42 library(stringr) 42 library(stringr)
43 library(tidyr) 43 library(tidyr)
44 44 library(curl)
45 if(!runExampleL) 45 library(jsonlite)
46 argLs <- parseCommandArgs(evaluate=FALSE) 46 library(stringi)
47
48 if (!runExampleL)
49 argLs <- parseCommandArgs(evaluate = FALSE)
47 logFile <- argLs[["logOut"]] 50 logFile <- argLs[["logOut"]]
48 sink(logFile) 51 sink(logFile)
49 52
50 cat("\tPACKAGE INFO\n") 53 cat("\tPACKAGE INFO\n")
51 sessionInfo() 54 sessionInfo()
52 55
53 ##------------------------------ 56 ##------------------------------
54 ## Functions 57 ## Functions
55 ##------------------------------ 58 ##------------------------------
56 source_local <- function(fname) 59 source_local <- function(fname) {
57 { 60 argv <- commandArgs(trailingOnly = FALSE)
58 argv <- commandArgs(trailingOnly = FALSE) 61 base_dir <- dirname(substring(argv[grep("--file=", argv)], 8))
59 base_dir <- dirname(substring(argv[grep("--file=", argv)], 8)) 62 source(paste(base_dir, fname, sep = "/"))
60 source(paste(base_dir, fname, sep="/"))
61 } 63 }
62 #Import the different functions 64 #Import the different functions
63 source_local("annotationRmn2D.R") 65 source_local("annotationRmn2D.R")
64 source_local("annotationRmn2DGlobale.R") 66 source_local("annotationRmn2DGlobale.R")
65 source_local("viridis.R") 67 source_local("viridis.R")
66 68
67 ## Input parameter values 69 ## Input parameter values
68 fileToAnnotate <- argLs[[1]] 70 fileToAnnotate <- argLs[[1]]
71 # Constraints values
72 ph <- argLs$pH
73 field <- argLs$magneticField
74
69 # Chosen sequence(s) 75 # Chosen sequence(s)
70 cosy <- 0 76 cosy <- 0
71 hmbc <- 0 77 hmbc <- 0
72 hsqc <- 0 78 hsqc <- 0
73 jres <- 0 79 jres <- 0
74 tocsy <- 0 80 tocsy <- 0
75 ## sequences <- str_split(argLs[[2]], ",")[[1]] 81
76 ## for (s in 1:length(sequences)) 82 if (argLs$cosy_2dsequences == "yes") {
77 ## {
78 ## argv <- commandArgs(trailingOnly = FALSE)
79 ## currentDir <- dirname(substring(argv[grep("--file=", argv)], 8))
80 ## if (sequences[s]=="cosy"){
81 ## cosy <- 1
82 ## load(paste(currentDir, "BdDReference_COSY.RData", sep="/"))
83 ## }else if(sequences[s]=="hmbc"){
84 ## hmbc <- 1
85 ## load(paste(currentDir, "BdDReference_HMBC.RData", sep="/"))
86 ## }else if(sequences[s]=="hsqc"){
87 ## hsqc <- 1
88 ## load(paste(currentDir, "BdDReference_HSQC.RData", sep="/"))
89 ## }else if(sequences[s]=="jres"){
90 ## jres <- 1
91 ## load(paste(currentDir, "BdDReference_JRES.RData", sep="/"))
92 ## }else if(sequences[s]=="tocsy"){
93 ## tocsy <- 1
94 ## load(paste(currentDir, "BdDReference_TOCSY.RData", sep="/"))
95 ## }else
96 ## stop("No chosen sequence", call.=FALSE)
97 ## }
98
99 if (argLs[[2]]=='yes')
100 {
101 argv <- commandArgs(trailingOnly = FALSE)
102 currentDir <- dirname(substring(argv[grep("--file=", argv)], 8))
103 cosy <- 1 83 cosy <- 1
104 load(paste(currentDir, "BdDReference_COSY.RData", sep="/")) 84 peakforestSpectra <- readLines(curl("https://metabohub.peakforest.org/rest/v1/spectra/nmr2d/search?query=cosy&token=9131jq9l8gsjn1j14t351h716u&max=500"))
105 } 85 peakforestSpectra <- fromJSON(peakforestSpectra, simplifyDataFrame = TRUE)
106 86 if (ph != 0)
107 if (argLs[[3]]=='yes') 87 peakforestSpectra <- peakforestSpectra[peakforestSpectra$sampleNMRTubeConditionsMetadata$potentiaHydrogenii == ph, ]
108 { 88 if (field != 0)
109 argv <- commandArgs(trailingOnly = FALSE) 89 peakforestSpectra <- peakforestSpectra[peakforestSpectra$analyzerNMRSpectrometerDevice$magneticFieldStrenght == field, ]
110 currentDir <- dirname(substring(argv[grep("--file=", argv)], 8)) 90
91 if (nrow(peakforestSpectra) != 0) {
92 BdDReference_COSY <- peakforestSpectra$peaks
93 names(BdDReference_COSY) <- str_split(peakforestSpectra[, 2], simplify = TRUE, pattern = ";")[, 1]
94 names(BdDReference_COSY) <- enc2utf8(names(BdDReference_COSY))
95 names(BdDReference_COSY) <- str_replace_all(names(BdDReference_COSY), "\u00e9", "e")
96
97 for (k in seq_len(length(BdDReference_COSY))) {
98 peakforestSpectra_df <- data.frame(ppm.dim1 = BdDReference_COSY[[k]][, 2], ppm.dim2 = BdDReference_COSY[[k]][, 1],
99 BdDReference_COSY[[k]][, 3:ncol(BdDReference_COSY[[k]])])
100 BdDReference_COSY[[k]] <- peakforestSpectra_df
101 }
102 } else {
103 stop("No COSY spectra correspond to requested pH and/or magnetic field", call. = FALSE)
104 }
105 rm(peakforestSpectra)
106 rm(peakforestSpectra_df)
107 }
108
109 if (argLs$hmbc_2dsequences == "yes") {
110 hmbc <- 1
111 peakforestSpectra <- readLines(curl("https://metabohub.peakforest.org/rest/v1/spectra/nmr2d/search?query=hmbc&token=9131jq9l8gsjn1j14t351h716u&max=500"))
112 peakforestSpectra <- fromJSON(peakforestSpectra, simplifyDataFrame = TRUE)
113 if (ph != 0)
114 peakforestSpectra <- peakforestSpectra[peakforestSpectra$sampleNMRTubeConditionsMetadata$potentiaHydrogenii == ph, ]
115 if (field != 0)
116 peakforestSpectra <- peakforestSpectra[peakforestSpectra$analyzerNMRSpectrometerDevice$magneticFieldStrenght == field, ]
117
118 if (nrow(peakforestSpectra) != 0) {
119
120 BdDReference_HMBC <- peakforestSpectra$peaks
121 names(BdDReference_HMBC) <- str_split(peakforestSpectra[, 2], simplify = TRUE, pattern = ";")[, 1]
122 names(BdDReference_HMBC) <- enc2utf8(names(BdDReference_HMBC))
123 names(BdDReference_HMBC) <- str_replace_all(names(BdDReference_HMBC), "\u00e9", "e")
124
125 peakforestSpectra_df <- data.frame()
126 for (k in seq_len(length(BdDReference_HMBC))) {
127 peakforestSpectra_df <- data.frame(ppm.dim1 = BdDReference_HMBC[[k]][, 2], ppm.dim2 = BdDReference_HMBC[[k]][, 1],
128 BdDReference_HMBC[[k]][, 3:ncol(BdDReference_HMBC[[k]])])
129 BdDReference_HMBC[[k]] <- peakforestSpectra_df
130 }
131 } else {
132 stop("No HMBC spectra correspond to requested pH and/or magnetic field", call. = FALSE)
133 }
134 rm(peakforestSpectra)
135 rm(peakforestSpectra_df)
136 }
137
138 if (argLs$hsqc_2dsequences == "yes") {
139 hsqc <- 1
140 peakforestSpectra <- readLines(curl("https://metabohub.peakforest.org/rest/v1/spectra/nmr2d/search?query=hsqc&token=9131jq9l8gsjn1j14t351h716u&max=500"))
141 peakforestSpectra <- fromJSON(peakforestSpectra, simplifyDataFrame = TRUE)
142
143 if (ph != 0)
144 peakforestSpectra <- peakforestSpectra[peakforestSpectra$sampleNMRTubeConditionsMetadata$potentiaHydrogenii == ph, ]
145 if (field != 0)
146 peakforestSpectra <- peakforestSpectra[peakforestSpectra$analyzerNMRSpectrometerDevice$magneticFieldStrenght == field, ]
147
148 if (nrow(peakforestSpectra) != 0) {
149 BdDReference_HSQC <- peakforestSpectra$peaks
150 names(BdDReference_HSQC) <- str_split(peakforestSpectra[, 2], simplify = TRUE, pattern = ";")[, 1]
151 names(BdDReference_HSQC) <- enc2utf8(names(BdDReference_HSQC))
152 names(BdDReference_HSQC) <- str_replace_all(names(BdDReference_HSQC), "\u00e9", "e")
153
154 for (k in seq_len(length(BdDReference_HSQC))) {
155 peakforestSpectra_df <- data.frame(ppm.dim1 = BdDReference_HSQC[[k]][, 2], ppm.dim2 = BdDReference_HSQC[[k]][, 1],
156 BdDReference_HSQC[[k]][, 3:ncol(BdDReference_HSQC[[k]])])
157 BdDReference_HSQC[[k]] <- peakforestSpectra_df
158 }
159 } else {
160 stop("No HSQC spectra correspond to requested pH and/or magnetic field", call. = FALSE)
161 }
162 rm(peakforestSpectra)
163 rm(peakforestSpectra_df)
164 }
165
166 if (argLs$jres_2dsequences == "yes") {
111 jres <- 1 167 jres <- 1
112 load(paste(currentDir, "BdDReference_JRES.RData", sep="/")) 168 peakforestSpectra <- readLines(curl("https://metabohub.peakforest.org/rest/v1/spectra/nmr2d/search?query=jres&token=9131jq9l8gsjn1j14t351h716u&max=500"))
113 } 169 peakforestSpectra <- fromJSON(peakforestSpectra, simplifyDataFrame = TRUE)
114 170
115 if (argLs[[4]]=='yes') 171 if (ph != 0)
116 { 172 peakforestSpectra <- peakforestSpectra[peakforestSpectra$sampleNMRTubeConditionsMetadata$potentiaHydrogenii == ph, ]
117 argv <- commandArgs(trailingOnly = FALSE) 173 if (field != 0)
118 currentDir <- dirname(substring(argv[grep("--file=", argv)], 8)) 174 peakforestSpectra <- peakforestSpectra[peakforestSpectra$analyzerNMRSpectrometerDevice$magneticFieldStrenght == field, ]
119 hmbc <- 1 175
120 load(paste(currentDir, "BdDReference_HMBC.RData", sep="/")) 176 if (nrow(peakforestSpectra) != 0) {
121 } 177 BdDReference_JRES <- peakforestSpectra$peaks
122 178 names(BdDReference_JRES) <- str_split(peakforestSpectra[, 2], simplify = TRUE, pattern = ";")[, 1]
123 if (argLs[[5]]=='yes') 179 names(BdDReference_JRES) <- enc2utf8(names(BdDReference_JRES))
124 { 180 names(BdDReference_JRES) <- str_replace_all(names(BdDReference_JRES), "\u00e9", "e")
125 argv <- commandArgs(trailingOnly = FALSE) 181
126 currentDir <- dirname(substring(argv[grep("--file=", argv)], 8)) 182 for (k in seq_len(length(BdDReference_JRES))) {
127 hsqc <- 1 183 peakforestSpectra_df <- data.frame(ppm.dim1 = BdDReference_JRES[[k]][, 2], ppm.dim2 = BdDReference_JRES[[k]][, 1],
128 load(paste(currentDir, "BdDReference_HSQC.RData", sep="/")) 184 BdDReference_JRES[[k]][, 3:ncol(BdDReference_JRES[[k]])])
129 } 185 BdDReference_JRES[[k]] <- peakforestSpectra_df
130 186 }
131 if (argLs[[6]]=='yes') 187 } else {
132 { 188 stop("No JRES spectra correspond to requested pH and/or magnetic field", call. = FALSE)
133 argv <- commandArgs(trailingOnly = FALSE) 189 }
134 currentDir <- dirname(substring(argv[grep("--file=", argv)], 8)) 190 rm(peakforestSpectra)
191 rm(peakforestSpectra_df)
192 }
193
194 if (argLs$tocsy_2dsequences == "yes") {
135 tocsy <- 1 195 tocsy <- 1
136 load(paste(currentDir, "BdDReference_TOCSY.RData", sep="/")) 196 peakforestSpectra <- readLines(curl("https://metabohub.peakforest.org/rest/v1/spectra/nmr2d/search?query=tocsy&token=9131jq9l8gsjn1j14t351h716u&max=500"))
137 } 197 peakforestSpectra <- fromJSON(peakforestSpectra, simplifyDataFrame = TRUE)
138 198
139 if (argLs[[2]]=='no' & argLs[[3]]=='no' & argLs[[4]]=='no' & argLs[[5]]=='no' & argLs[[6]]=='no') 199 if (ph != 0)
140 stop("No chosen sequence", call.=FALSE) 200 peakforestSpectra <- peakforestSpectra[peakforestSpectra$sampleNMRTubeConditionsMetadata$potentiaHydrogenii == ph, ]
201 if (field != 0)
202 peakforestSpectra <- peakforestSpectra[peakforestSpectra$analyzerNMRSpectrometerDevice$magneticFieldStrenght == field, ]
203
204 if (nrow(peakforestSpectra) != 0) {
205 BdDReference_TOCSY <- peakforestSpectra$peaks
206 names(BdDReference_TOCSY) <- str_split(peakforestSpectra[, 2], simplify = TRUE, pattern = ";")[, 1]
207 names(BdDReference_TOCSY) <- enc2utf8(names(BdDReference_TOCSY))
208 names(BdDReference_TOCSY) <- str_replace_all(names(BdDReference_TOCSY), "\u00e9", "e")
209
210 for (k in seq_len(length(BdDReference_TOCSY))) {
211 peakforestSpectra_df <- data.frame(ppm.dim1 = BdDReference_TOCSY[[k]][, 2], ppm.dim2 = BdDReference_TOCSY[[k]][, 1],
212 BdDReference_TOCSY[[k]][, 3:ncol(BdDReference_TOCSY[[k]])])
213 BdDReference_TOCSY[[k]] <- peakforestSpectra_df
214 }
215 } else {
216 stop("No TOCSY spectra correspond to requested pH and/or magnetic field", call. = FALSE)
217 }
218 rm(peakforestSpectra)
219 rm(peakforestSpectra_df)
220 }
221
222 if (argLs$cosy_2dsequences == "no" & argLs$hmbc_2dsequences == "no" & argLs$hsqc_2dsequences == "no" & argLs$jres_2dsequences == "no" & argLs$tocsy_2dsequences == "no")
223 stop("No chosen sequence. You have to choose at least 1 sequence", call. = FALSE)
141 224
142 225
143 # User database 226 # User database
144 227
145 228
156 AnnotationGraph <- argLs[["AnnotationGraph"]] 239 AnnotationGraph <- argLs[["AnnotationGraph"]]
157 240
158 print(argLs) 241 print(argLs)
159 242
160 ## ANNOTATION 243 ## ANNOTATION
161 st0=Sys.time() 244 st0 <- Sys.time()
162 pdf(AnnotationGraph,onefile=TRUE) 245 pdf(AnnotationGraph, onefile = TRUE)
163 annotationMelange <- annotationRmn2DGlobale(fileToAnnotate, tolPpm1=tolPpm1, tolPpm2HJRes=tolPpm2HJRes, 246 annotationMelange <- annotationRmn2DGlobale(fileToAnnotate, tolPpm1 = tolPpm1, tolPpm2HJRes = tolPpm2HJRes,
164 tolPpm2C=tolPpm2C, cosy=cosy, hmbc=hmbc, hsqc=hsqc, 247 tolPpm2C = tolPpm2C, cosy = cosy, hmbc = hmbc, hsqc = hsqc,
165 jres=jres, tocsy=tocsy, seuil=seuil, unicite=unicite) 248 jres = jres, tocsy = tocsy, seuil = seuil, unicite = unicite)
166 dev.off() 249 dev.off()
167 250
168 if (cosy==1) 251 if (cosy == 1) {
169 { 252 write.table(annotationMelange$COSY$liste_resultat, file = argLs[["annotationCOSY"]], quote = FALSE,
170 write.table(annotationMelange$COSY$liste_resultat, file=argLs[["annotationCOSY"]], quote=FALSE, 253 row.names = FALSE, sep = "\t")
171 row.names=FALSE,sep="\t")
172 if (nrow(annotationMelange$COSY$listing_ppm_commun) != 0) 254 if (nrow(annotationMelange$COSY$listing_ppm_commun) != 0)
173 write.table(annotationMelange$COSY$listing_ppm_commun, file=argLs[["ppmCommunCOSY"]], quote=FALSE, 255 write.table(annotationMelange$COSY$listing_ppm_commun, file = argLs[["ppmCommunCOSY"]], quote = FALSE,
174 row.names=FALSE,sep="\t") 256 row.names = FALSE, sep = "\t")
175 } 257 }
176 258
177 if (hmbc==1) 259 if (hmbc == 1) {
178 { 260 write.table(annotationMelange$HMBC$liste_resultat, file = argLs[["annotationHMBC"]], quote = FALSE,
179 write.table(annotationMelange$HMBC$liste_resultat, file=argLs[["annotationHMBC"]], quote=FALSE, 261 row.names = FALSE, sep = "\t")
180 row.names=FALSE,sep="\t")
181 if (nrow(annotationMelange$HMBC$listing_ppm_commun) != 0) 262 if (nrow(annotationMelange$HMBC$listing_ppm_commun) != 0)
182 write.table(annotationMelange$HMBC$listing_ppm_commun, file=argLs[["ppmCommunHMBC"]], quote=FALSE, 263 write.table(annotationMelange$HMBC$listing_ppm_commun, file = argLs[["ppmCommunHMBC"]], quote = FALSE,
183 row.names=FALSE,sep="\t") 264 row.names = FALSE, sep = "\t")
184 } 265 }
185 266
186 if (hsqc==1) 267 if (hsqc == 1) {
187 { 268 write.table(annotationMelange$HSQC$liste_resultat, file = argLs[["annotationHSQC"]], quote = FALSE,
188 write.table(annotationMelange$HSQC$liste_resultat, file=argLs[["annotationHSQC"]], quote=FALSE, 269 row.names = FALSE, sep = "\t")
189 row.names=FALSE,sep="\t")
190 if (nrow(annotationMelange$HSQC$listing_ppm_commun) != 0) 270 if (nrow(annotationMelange$HSQC$listing_ppm_commun) != 0)
191 write.table(annotationMelange$HSQC$listing_ppm_commun, file=argLs[["ppmCommunHSQC"]], quote=FALSE, 271 write.table(annotationMelange$HSQC$listing_ppm_commun, file = argLs[["ppmCommunHSQC"]], quote = FALSE,
192 row.names=FALSE,sep="\t") 272 row.names = FALSE, sep = "\t")
193 } 273 }
194 274
195 if (jres==1) 275 if (jres == 1) {
196 { 276 write.table(annotationMelange$JRES$liste_resultat, file = argLs[["annotationJRES"]], quote = FALSE,
197 write.table(annotationMelange$JRES$liste_resultat, file=argLs[["annotationJRES"]], quote=FALSE, 277 row.names = FALSE, sep = "\t")
198 row.names=FALSE,sep="\t")
199 if (nrow(annotationMelange$JRES$listing_ppm_commun) != 0) 278 if (nrow(annotationMelange$JRES$listing_ppm_commun) != 0)
200 write.table(annotationMelange$JRES$listing_ppm_commun, file=argLs[["ppmCommunJRES"]], quote=FALSE, 279 write.table(annotationMelange$JRES$listing_ppm_commun, file = argLs[["ppmCommunJRES"]], quote = FALSE,
201 row.names=FALSE,sep="\t") 280 row.names = FALSE, sep = "\t")
202 } 281 }
203 282
204 if (tocsy==1) 283 if (tocsy == 1) {
205 { 284 write.table(annotationMelange$TOCSY$liste_resultat, file = argLs[["annotationTOCSY"]], quote = FALSE,
206 write.table(annotationMelange$TOCSY$liste_resultat, file=argLs[["annotationTOCSY"]], quote=FALSE, 285 row.names = FALSE, sep = "\t")
207 row.names=FALSE,sep="\t")
208 if (nrow(annotationMelange$TOCSY$listing_ppm_commun) != 0) 286 if (nrow(annotationMelange$TOCSY$listing_ppm_commun) != 0)
209 write.table(annotationMelange$TOCSY$listing_ppm_commun, file=argLs[["ppmCommunTOCSY"]], quote=FALSE, 287 write.table(annotationMelange$TOCSY$listing_ppm_commun, file = argLs[["ppmCommunTOCSY"]], quote = FALSE,
210 row.names=FALSE,sep="\t") 288 row.names = FALSE, sep = "\t")
211 } 289 }
212 290
213 ## Combinaison de sequences 291 ## Combinaison de sequences
214 if (cosy + jres + hmbc + hsqc + tocsy > 1) 292 if (cosy + jres + hmbc + hsqc + tocsy > 1) {
215 { 293 write.table(annotationMelange$combination, file = argLs[["annotationCombination"]], quote = FALSE,
216 write.table(annotationMelange$combination, file=argLs[["annotationCombination"]], quote=FALSE, 294 row.names = FALSE, sep = "\t")
217 row.names=FALSE,sep="\t") 295 }
218 } 296 st1 <- Sys.time()
219 st1=Sys.time() 297 print(st1 - st0)
220 print(st1-st0)
221 298
222 ## Ending 299 ## Ending
223 ##-------- 300 ##--------
224 cat("\nEnd of '2D NMR annotation' Galaxy module call: ", as.character(Sys.time()), sep = "") 301 cat("\nEnd of '2D NMR annotation' Galaxy module call: ", as.character(Sys.time()), sep = "")
225 sink() 302 sink()