comparison asics_wrapper.R @ 4:0ff2d9211ebe draft default tip

planemo upload for repository https://github.com/workflow4metabolomics/nmr_annotation commit 3791815505685d0038e392a702860843fe1d443e
author marie-tremblay-metatoul
date Fri, 21 Sep 2018 07:42:53 -0400
parents b55559a2854f
children
comparison
equal deleted inserted replaced
3:9b45c0a33573 4:0ff2d9211ebe
1 #!/usr/local/public/bin/Rscript --vanilla --slave --no-site-file
2
3 ## 29122017_asics_wrapper.R
4 ## Remi Servien, Patrick Tardivel, Marie Tremblay-Franco and Gaelle Lefort
5 ## marie.tremblay-franco@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 library(ASICS)
22
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 #Import the different functions
33 source_local("DrawSpec.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 # Standards loading
55 load(argLs[["standards"]])
56
57 ## Parameters Loading
58 ##-------------------
59 # Inputs
60 ## Spectrum to annotate
61 zipfile= argLs[["zipfile"]]
62 directory=unzip(zipfile, list=F)
63 directory=paste(getwd(),strsplit(directory[1],"/")[[1]][2],sep="/")
64
65
66 ##Exclusion zone(s)
67 exclusionZones <- argLs[["zone_exclusion_choices.choice"]]
68 exclusionZonesBorders <- NULL
69 if (!is.null(argLs$zone_exclusion_left))
70 {
71 for(i in which(names(argLs)=="zone_exclusion_left"))
72 {
73 # exclusionZonesBorders <- c(exclusionZonesBorders,list(c(argLs[[i]],argLs[[i+1]])))
74 exclusionZonesBorders <- c(exclusionZonesBorders,argLs[[i]],argLs[[i+1]])
75 }
76 exclusionZonesBorders <- matrix(exclusionZonesBorders, byrow=T, ncol=2)
77 }
78 if (is.null(argLs$zone_exclusion_left))
79 {
80 exclusionZonesBorders <- matrix(c(0,0), byrow=T, ncol=2)
81 }
82 ## Maximal allowed shift
83 shift <- argLs[["shift"]]
84
85 ## Graphical zone(s)
86 graphicalZones <- argLs[["zone_graphical_choices.choice"]]
87 graphicalZonesBorders <- NULL
88 if (!is.null(argLs$zone_exclusion_left))
89 {
90 for(i in which(names(argLs)=="zone_graphical_left"))
91 {
92 graphicalZonesBorders <- c(graphicalZonesBorders,list(c(argLs[[i]],argLs[[i+1]])))
93 }
94 }
95
96 # Outputs
97 logOut <- argLs[["logOut"]]
98 proportionEstimation <- argLs[["proportionEstimation"]]
99 graphOut <- argLs[["graphOut"]]
100
101 sink(logOut)
102 cat("\tPACKAGE INFO\n")
103 # pkgs=c("batch", "ASICS")
104 pkgs=c("batch", "ASICS")
105 for(pkg in pkgs) {
106 suppressPackageStartupMessages( stopifnot( library(pkg, quietly=TRUE, logical.return=TRUE, character.only=TRUE)))
107 cat(pkg,"\t",as.character(packageVersion(pkg)),"\n",sep="")
108 }
109 cat("\n")
110
111
112 ## Checking arguments
113 ##-------------------
114 error.stock <- "\n"
115 if(length(error.stock) > 1)
116 stop(error.stock)
117
118
119 ## Computation
120 ##------------
121 # annotation.Asics <- ASICS(directory, exclusion.areas=matrix(exclusionZonesBorders, byrow=T, ncol=2),
122 # max.shift=shift, which.spectra="last", library.metabolites=NULL,
123 # threshold.noise=0.02, seed=1234, nb.iter.signif=400)
124 annotation.Asics <- ASICS(directory, exclusion.areas=exclusionZonesBorders,
125 max.shift=shift, which.spectra="last", library.metabolites=NULL,
126 threshold.noise=0.02, seed=1234, nb.iter.signif=400)
127
128
129 ## Saving
130 ##-------
131 # Identified metabolites
132 metabolites.estimation <- present_metabolites(annotation.Asics)
133 colnames(metabolites.estimation) <- c("Metabolite",colnames(metabolites.estimation)[-1])
134 write.table(metabolites.estimation,file=argLs$proportionEstimation,row.names=FALSE,quote=FALSE,sep="\t")
135
136
137 ## Graphical display
138 ##------------------
139 # Raw and annotated spectra comparison
140 pdf(graphOut,onefile=TRUE)
141
142 ## Graphical output: overlay of raw and estimated spectra
143 ppm.metabolites.estimation <- data.frame(round(ppm_grid(annotation.Asics),3),
144 original_mixture(annotation.Asics))
145 colnames(ppm.metabolites.estimation) <- c("PPM", "EstimatedProportion")
146 ppm.metabolites.estimation <- ppm.metabolites.estimation[order(ppm.metabolites.estimation[,1],decreasing=T), ]
147
148 mix <- data.frame(t(ppm.metabolites.estimation[,2]))
149 colnames(mix) <- ppm.metabolites.estimation[,1]
150 ppm <- ppm.metabolites.estimation[,1]
151
152 estimatedMix <- data.frame(round(ppm_grid(annotation.Asics),3), reconstituted_mixture(annotation.Asics))
153 colnames(estimatedMix) <- c("PPM","EstimatedProportion")
154 estimatedMix <- estimatedMix[order(estimatedMix[,1],decreasing=T), ]
155 estimatedMix <- estimatedMix[,2]
156
157 ## Whole spectra
158 GraphRange <- 1:ncol(mix)
159 tempVal <- trunc(length(GraphRange)/10)
160 xPos <- c(10:0) * tempVal
161 plot(1:ncol(mix), mix, type='l', xlab="", main="", xaxt="n", ylab="")
162 axis(1, at=xPos, labels=colnames(mix)[xPos + 1])
163 lines(estimatedMix, col="red")
164 legend("topleft",legend=c("Real Mixture","Estimated Composition"),lty=c(1,1),col=c("black","red"))
165
166 ## Zoomed spectral window depending on user-selected zone(s)
167 graphical.zone.length <- length(graphicalZonesBorders)
168 if (graphical.zone.length != 0)
169
170 # par(mfrow=c(2,1))
171 for (g in 1:graphical.zone.length)
172 {
173 print(g)
174 plot(1:length((which(round(as.numeric(colnames(mix)),2) == graphicalZonesBorders[[g]][1])[1]):(which(round(as.numeric(colnames(mix)),2) == max(graphicalZonesBorders[[g]][2],0.5))[1])),
175 mix[(which(round(as.numeric(colnames(mix)),2) == graphicalZonesBorders[[g]][1])[1]):(which(round(as.numeric(colnames(mix)),2) == max(graphicalZonesBorders[[g]][2],0.5))[1])], type='l', xlab="", ylab="Intensity", main="", xaxt="n")
176 lines(estimatedMix[(which(round(as.numeric(colnames(mix)),2) == graphicalZonesBorders[[g]][1])[1]):(which(round(as.numeric(colnames(mix)),2) == max(graphicalZonesBorders[[g]][2],0.5))[1])],col="red")
177
178 xPos <- 1
179 nAxisPos <- 4
180 startP <- length(nAxisPos)
181 endP <- length((which(round(as.numeric(colnames(mix)),2) == graphicalZonesBorders[[g]][1])[1]):(which(round(as.numeric(colnames(mix)),2) == max(graphicalZonesBorders[[g]][2],0.5))[1]))
182 GraphRange <- c(startP:endP)
183 tempVal <- trunc(length(GraphRange)/nAxisPos)
184 xPos <- c(0:nAxisPos) * tempVal
185 noms <- ppm.metabolites.estimation[xPos + which(ppm == round(graphicalZonesBorders[[g]][1],1))[1],1]
186 axis(1, at=xPos, labels=noms)
187 }
188
189 invisible(dev.off())
190
191
192 ## Ending
193 ##---------------------
194 cat("\nEnd of 'NMR annotation' Galaxy module call: ", as.character(Sys.time()), sep="")
195 options(stringsAsFactors=strAsFacL)
196 rm(list=ls())
197 sink()
198