comparison lib_metams.r @ 4:c10824185547 draft

planemo upload for repository https://github.com/workflow4metabolomics/metaMS commit 174a1713024f246c1485cbd75218577e89353522
author workflow4metabolomics
date Wed, 03 Jul 2019 05:14:32 -0400
parents c75532b75ba1
children b8d4129dd2a6
comparison
equal deleted inserted replaced
3:c75532b75ba1 4:c10824185547
1 # lib_metams.r version 2.0.0 1 # lib_metams.r version 2.1.1
2 # R function for metaMS runGC under W4M 2 # R function for metaMS runGC under W4M
3 # author Yann GUITTON CNRS IRISA/LINA Idealg project 2014-2015 3 # author Yann GUITTON CNRS IRISA/LINA Idealg project 2014-2015
4 # author Yann GUITTON Oniris Laberca 2015-2017 4 # author Yann GUITTON Oniris Laberca 2015-2017
5 5
6 6
7 #@author G. Le Corguille
8 # This function will
9 # - load the packages
10 # - display the sessionInfo
11 loadAndDisplayPackages <- function(pkgs) {
12 for(pkg in pkgs) suppressPackageStartupMessages( stopifnot( library(pkg, quietly=TRUE, logical.return=TRUE, character.only=TRUE)))
13
14 sessioninfo = sessionInfo()
15 cat(sessioninfo$R.version$version.string,"\n")
16 cat("Main packages:\n")
17 for (pkg in names(sessioninfo$otherPkgs)) { cat(paste(pkg,packageVersion(pkg)),"\t") }; cat("\n")
18 cat("Other loaded packages:\n")
19 for (pkg in names(sessioninfo$loadedOnly)) { cat(paste(pkg,packageVersion(pkg)),"\t") }; cat("\n")
20 }
21
22 #This function list the compatible files within the directory as xcms did
23 #@author Gildas Le Corguille lecorguille@sb-roscoff.fr ABiMS TEAM
24 getMSFiles <- function (directory) {
25 filepattern <- c("[Cc][Dd][Ff]", "[Nn][Cc]", "([Mm][Zz])?[Xx][Mm][Ll]","[Mm][Zz][Dd][Aa][Tt][Aa]", "[Mm][Zz][Mm][Ll]")
26 filepattern <- paste(paste("\\.", filepattern, "$", sep=""),collapse="|")
27 info <- file.info(directory)
28 listed <- list.files(directory[info$isdir], pattern=filepattern,recursive=TRUE, full.names=TRUE)
29 files <- c(directory[!info$isdir], listed)
30 exists <- file.exists(files)
31 files <- files[exists]
32 return(files)
33 }
34
35 # This function retrieve a xset like object
36 #@author Gildas Le Corguille lecorguille@sb-roscoff.fr
37 getxcmsSetObject <- function(xobject) {
38 # XCMS 1.x
39 if (class(xobject) == "xcmsSet")
40 return (xobject)
41 # XCMS 3.x
42 if (class(xobject) == "XCMSnExp") {
43 # Get the legacy xcmsSet object
44 suppressWarnings(xset <- as(xobject, 'xcmsSet'))
45 if (!is.null(xset@phenoData$sample_group))
46 sampclass(xset) <- xset@phenoData$sample_group
47 else
48 sampclass(xset) <- "."
49 return (xset)
50 }
51 }
52
53 #J.Saint-Vanne
54 #Function to correct the file names which can be written like "..alg8.mzData" and we just want "alg8"
55 getCorrectFileName <- function(peaktable,sampleMetadata){
56
57 #Try to start for the first sample, avoid description of line with colnamesdontwant
58 i <- 1
59 while(!(sampleMetadata[1,1] %in% strsplit(colnames(peaktable)[i],"\\.")[[1]])) {
60 if(i < length(peaktable)) {
61 i <- i + 1
62 } else {
63 break
64 }
65 }
66 #i now correspond to the first column with a sample
67 for(j in 1:(nrow(sampleMetadata))) {
68 col <- j + i-1 #minus 1 cause i is the good column to start and j start at 1
69 if(col <= length(colnames(peaktable))) {
70 newname <- gsub("(^.*)(\\..*$)","\\1",colnames(peaktable)[col])
71 if(newname != sampleMetadata[j,1]){
72 #Correction for 2 points starting the name (I don't know why they are here...)
73 if(".." == gsub("(^\\.+)(.*)","\\1",newname)){
74 newname <- sub("(^\\.+)(.*)","\\2",newname)
75 }
76 }
77 colnames(peaktable)[col] <- newname
78 }
79 }
80 return(peaktable)
81 }
82
83
84 #@author Gildas Le Corguille lecorguille@sb-roscoff.fr
85 # This function get the raw file path from the arguments
86 getRawfilePathFromArguments <- function(singlefile, zipfile, listArguments) {
87 if (!is.null(listArguments[["zipfile"]])) zipfile = listArguments[["zipfile"]]
88 if (!is.null(listArguments[["zipfilePositive"]])) zipfile = listArguments[["zipfilePositive"]]
89 if (!is.null(listArguments[["zipfileNegative"]])) zipfile = listArguments[["zipfileNegative"]]
90
91 if (!is.null(listArguments[["singlefile_galaxyPath"]])) {
92 singlefile_galaxyPaths = listArguments[["singlefile_galaxyPath"]];
93 singlefile_sampleNames = listArguments[["singlefile_sampleName"]]
94 }
95 if (!is.null(listArguments[["singlefile_galaxyPathPositive"]])) {
96 singlefile_galaxyPaths = listArguments[["singlefile_galaxyPathPositive"]];
97 singlefile_sampleNames = listArguments[["singlefile_sampleNamePositive"]]
98 }
99 if (!is.null(listArguments[["singlefile_galaxyPathNegative"]])) {
100 singlefile_galaxyPaths = listArguments[["singlefile_galaxyPathNegative"]];
101 singlefile_sampleNames = listArguments[["singlefile_sampleNameNegative"]]
102 }
103 if (exists("singlefile_galaxyPaths")){
104 singlefile_galaxyPaths = unlist(strsplit(singlefile_galaxyPaths,","))
105 singlefile_sampleNames = unlist(strsplit(singlefile_sampleNames,","))
106
107 singlefile=NULL
108 for (singlefile_galaxyPath_i in seq(1:length(singlefile_galaxyPaths))) {
109 singlefile_galaxyPath=singlefile_galaxyPaths[singlefile_galaxyPath_i]
110 singlefile_sampleName=singlefile_sampleNames[singlefile_galaxyPath_i]
111 singlefile[[singlefile_sampleName]] = singlefile_galaxyPath
112 }
113 }
114 for (argument in c("zipfile","zipfilePositive","zipfileNegative","singlefile_galaxyPath","singlefile_sampleName","singlefile_galaxyPathPositive","singlefile_sampleNamePositive","singlefile_galaxyPathNegative","singlefile_sampleNameNegative")) {
115 listArguments[[argument]]=NULL
116 }
117 return(list(zipfile=zipfile, singlefile=singlefile, listArguments=listArguments))
118 }
119
120
121 #@author Gildas Le Corguille lecorguille@sb-roscoff.fr
122 # This function retrieve the raw file in the working directory
123 # - if zipfile: unzip the file with its directory tree
124 # - if singlefiles: set symlink with the good filename
125 retrieveRawfileInTheWorkingDirectory <- function(singlefile, zipfile) {
126 if(!is.null(singlefile) && (length("singlefile")>0)) {
127 for (singlefile_sampleName in names(singlefile)) {
128 singlefile_galaxyPath = singlefile[[singlefile_sampleName]]
129 if(!file.exists(singlefile_galaxyPath)){
130 error_message=paste("Cannot access the sample:",singlefile_sampleName,"located:",singlefile_galaxyPath,". Please, contact your administrator ... if you have one!")
131 print(error_message); stop(error_message)
132 }
133 file.symlink(singlefile_galaxyPath,singlefile_sampleName)
134 }
135 directory = "."
136 }
137 if(!is.null(zipfile) && (zipfile!="")) {
138 if(!file.exists(zipfile)){
139 error_message=paste("Cannot access the Zip file:",zipfile,". Please, contact your administrator ... if you have one!")
140 print(error_message)
141 stop(error_message)
142 }
143
144 #list all file in the zip file
145 #zip_files=unzip(zipfile,list=T)[,"Name"]
146
147 #unzip
148 suppressWarnings(unzip(zipfile, unzip="unzip"))
149
150 #get the directory name
151 filesInZip=unzip(zipfile, list=T);
152 directories=unique(unlist(lapply(strsplit(filesInZip$Name,"/"), function(x) x[1])));
153 directories=directories[!(directories %in% c("__MACOSX")) & file.info(directories)$isdir]
154 directory = "."
155 if (length(directories) == 1) directory = directories
156
157 cat("files_root_directory\t",directory,"\n")
158 }
159 return (directory)
160 }
161
7 ##ADDITIONS FROM Y. Guitton 162 ##ADDITIONS FROM Y. Guitton
8 getBPC <- function(file,rtcor=NULL, ...) { 163 getBPC <- function(file,rtcor=NULL, ...) {
9 object <- xcmsRaw(file) 164 object <- xcmsRaw(file)
10 sel <- profRange(object, ...) 165 sel <- profRange(object, ...)
11 cbind(if (is.null(rtcor)) object@scantime[sel$scanidx] else rtcor ,xcms:::colMax(object@env$profile[sel$massidx,sel$scanidx,drop=FALSE])) 166 cbind(if (is.null(rtcor)) object@scantime[sel$scanidx] else rtcor ,xcms:::colMax(object@env$profile[sel$massidx,sel$scanidx,drop=FALSE]))
12 167 }
13 } 168
14 169 getBPC2s <- function (files, xset = NULL, pdfname="BPCs.pdf", rt = c("raw","corrected"), scanrange=NULL) {
15 getBPC2s <- function (files, pdfname="BPCs.pdf", rt = c("raw","corrected"), scanrange=NULL) {
16 require(xcms) 170 require(xcms)
17 171
18
19
20 #create sampleMetadata, get sampleMetadata and class 172 #create sampleMetadata, get sampleMetadata and class
21 sampleMetadata<-xcms:::phenoDataFromPaths(files) 173 if(!is.null(xset)) {
22 class<-class<-as.vector(levels(sampleMetadata[,"class"])) #create phenoData like table 174 #When files come from XCMS3 directly before metaMS
175 sampleMetadata <- xset@phenoData
176 } else {
177 #When files come from a zip directory with raw files before metaMS
178 sampleMetadata<-xcms:::phenoDataFromPaths(files)
179 }
180 class<-unique(sampleMetadata[,"class"]) #create phenoData like table
23 classnames<-vector("list",length(class)) 181 classnames<-vector("list",length(class))
24 for (i in 1:length(class)){ 182 for (i in 1:length(class)){
25 classnames[[i]]<-which( sampleMetadata[,1]==class[i]) 183 classnames[[i]]<-which( sampleMetadata[,"class"]==class[i])
26 } 184 }
27 185
28 N <- dim(sampleMetadata)[1] 186 N <- dim(sampleMetadata)[1]
29
30
31
32 TIC <- vector("list",N) 187 TIC <- vector("list",N)
33 188
34
35 for (j in 1:N) { 189 for (j in 1:N) {
36
37 cat(files[j],"\n")
38 TIC[[j]] <- getBPC(files[j]) 190 TIC[[j]] <- getBPC(files[j])
39 #good for raw 191 #good for raw
40 # seems strange for corrected 192 # seems strange for corrected
41 #errors if scanrange used in xcmsSetgeneration 193 #errors if scanrange used in xcmsSetgeneration
42 if (!is.null(xcmsSet) && rt == "corrected") 194 if (!is.null(xcmsSet) && rt == "corrected"){
43 rtcor <- xcmsSet@rt$corrected[[j]] else 195 rtcor <- xcmsSet@rt$corrected[[j]]
44 rtcor <- NULL 196 }else{
45 TIC[[j]] <- getBPC(files[j],rtcor=rtcor) 197 rtcor <- NULL
198 }
199 TIC[[j]] <- getBPC(files[j],rtcor=rtcor)
46 } 200 }
47 201
48 pdf(pdfname,w=16,h=10) 202 pdf(pdfname,w=16,h=10)
49 cols <- rainbow(N) 203 cols <- rainbow(N)
50 lty = 1:N 204 lty = 1:N
52 #search for max x and max y in BPCs 206 #search for max x and max y in BPCs
53 xlim = range(sapply(TIC, function(x) range(x[,1]))) 207 xlim = range(sapply(TIC, function(x) range(x[,1])))
54 ylim = range(sapply(TIC, function(x) range(x[,2]))) 208 ylim = range(sapply(TIC, function(x) range(x[,2])))
55 ylim = c(-ylim[2], ylim[2]) 209 ylim = c(-ylim[2], ylim[2])
56 210
57
58 ##plot start
59
60 if (length(class)>2){
61 for (k in 1:(length(class)-1)){
62 for (l in (k+1):length(class)){
63 print(paste(class[k],"vs",class[l],sep=" "))
64 plot(0, 0, type="n", xlim = xlim/60, ylim = ylim, main = paste("Base Peak Chromatograms \n","BPCs_",class[k]," vs ",class[l], sep=""), xlab = "Retention Time (min)", ylab = "BPC")
65 colvect<-NULL
66 for (j in 1:length(classnames[[k]])) {
67
68 tic <- TIC[[classnames[[k]][j]]]
69 # points(tic[,1]/60, tic[,2], col = cols[i], pch = pch[i], type="l")
70 points(tic[,1]/60, tic[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l")
71 colvect<-append(colvect,cols[classnames[[k]][j]])
72 }
73 for (j in 1:length(classnames[[l]])) {
74 # i=class2names[j]
75 tic <- TIC[[classnames[[l]][j]]]
76 points(tic[,1]/60, -tic[,2], col = cols[classnames[[l]][j]], pch = pch[classnames[[l]][j]], type="l")
77 colvect<-append(colvect,cols[classnames[[l]][j]])
78 }
79 legend("topright",paste(basename(files[c(classnames[[k]],classnames[[l]])])), col = colvect, lty = lty, pch = pch)
80 }
81 }
82 }#end if length >2
83 if (length(class)==2){
84 k=1
85 l=2
86 colvect<-NULL
87 plot(0, 0, type="n", xlim = xlim/60, ylim = ylim, main = paste("Base Peak Chromatograms \n","BPCs_",class[k],"vs",class[l], sep=""), xlab = "Retention Time (min)", ylab = "BPC")
88
89 for (j in 1:length(classnames[[k]])) {
90
91 tic <- TIC[[classnames[[k]][j]]]
92 # points(tic[,1]/60, tic[,2], col = cols[i], pch = pch[i], type="l")
93 points(tic[,1]/60, tic[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l")
94 colvect<-append(colvect,cols[classnames[[k]][j]])
95 }
96 for (j in 1:length(classnames[[l]])) {
97 # i=class2names[j]
98 tic <- TIC[[classnames[[l]][j]]]
99 points(tic[,1]/60, -tic[,2], col = cols[classnames[[l]][j]], pch = pch[classnames[[l]][j]], type="l")
100 colvect<-append(colvect,cols[classnames[[l]][j]])
101 }
102 legend("topright",paste(basename(files[c(classnames[[k]],classnames[[l]])])), col = colvect, lty = lty, pch = pch)
103
104 }#end length ==2
105 if (length(class)==1){
106 k=1
107 ylim = range(sapply(TIC, function(x) range(x[,2])))
108 colvect<-NULL
109 plot(0, 0, type="n", xlim = xlim/60, ylim = ylim, main = paste("Base Peak Chromatograms \n","BPCs_",class[k], sep=""), xlab = "Retention Time (min)", ylab = "BPC")
110
111 for (j in 1:length(classnames[[k]])) {
112
113 tic <- TIC[[classnames[[k]][j]]]
114 # points(tic[,1]/60, tic[,2], col = cols[i], pch = pch[i], type="l")
115 points(tic[,1]/60, tic[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l")
116 colvect<-append(colvect,cols[classnames[[k]][j]])
117 }
118
119 legend("topright",paste(basename(files[c(classnames[[k]])])), col = colvect, lty = lty, pch = pch)
120
121 }#end length ==1
122 dev.off()
123
124 # invisible(TIC)
125 }
126
127 getTIC <- function(file,rtcor=NULL) {
128 object <- xcmsRaw(file)
129 cbind(if (is.null(rtcor)) object@scantime else rtcor, rawEIC(object,mzrange=range(object@env$mz))$intensity)
130 }
131
132 ##
133 ## overlay TIC from all files in current folder or from xcmsSet, create pdf
134 ##
135 getTIC2s <- function(files, pdfname="TICs.pdf", rt=c("raw","corrected")) {
136
137 #create sampleMetadata, get sampleMetadata and class
138 sampleMetadata<-xcms:::phenoDataFromPaths(files)
139 class<-class<-as.vector(levels(sampleMetadata[,"class"])) #create phenoData like table
140 classnames<-vector("list",length(class))
141 for (i in 1:length(class)){
142 classnames[[i]]<-which( sampleMetadata[,1]==class[i])
143 }
144
145 N <- dim(sampleMetadata)[1]
146 TIC <- vector("list",N)
147
148 for (i in 1:N) {
149 cat(files[i],"\n")
150 if (!is.null(xcmsSet) && rt == "corrected")
151 rtcor <- xcmsSet@rt$corrected[[i]]
152 else
153 rtcor <- NULL
154 TIC[[i]] <- getTIC(files[i],rtcor=rtcor)
155 }
156
157 pdf(pdfname,w=16,h=10)
158 cols <- rainbow(N)
159 lty = 1:N
160 pch = 1:N
161 #search for max x and max y in TICs
162 xlim = range(sapply(TIC, function(x) range(x[,1])))
163 ylim = range(sapply(TIC, function(x) range(x[,2])))
164 ylim = c(-ylim[2], ylim[2])
165
166
167 ##plot start 211 ##plot start
168 if (length(class)>2){ 212 if (length(class)>2){
169 for (k in 1:(length(class)-1)){ 213 for (k in 1:(length(class)-1)){
170 for (l in (k+1):length(class)){ 214 for (l in (k+1):length(class)){
171 print(paste(class[k],"vs",class[l],sep=" ")) 215 cat(paste(class[k],"vs",class[l],sep=" ","\n"))
172 plot(0, 0, type="n", xlim = xlim/60, ylim = ylim, main = paste("Total Ion Chromatograms \n","TICs_",class[k]," vs ",class[l], sep=""), xlab = "Retention Time (min)", ylab = "TIC") 216 plot(0, 0, type="n", xlim = xlim/60, ylim = ylim, main = paste("Base Peak Chromatograms \n","BPCs_",class[k]," vs ",class[l], sep=""), xlab = "Retention Time (min)", ylab = "BPC")
173 colvect<-NULL 217 colvect<-NULL
174 for (j in 1:length(classnames[[k]])) { 218 for (j in 1:length(classnames[[k]])) {
175
176 tic <- TIC[[classnames[[k]][j]]] 219 tic <- TIC[[classnames[[k]][j]]]
177 # points(tic[,1]/60, tic[,2], col = cols[i], pch = pch[i], type="l") 220 # points(tic[,1]/60, tic[,2], col = cols[i], pch = pch[i], type="l")
178 points(tic[,1]/60, tic[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l") 221 points(tic[,1]/60, tic[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l")
179 colvect<-append(colvect,cols[classnames[[k]][j]]) 222 colvect<-append(colvect,cols[classnames[[k]][j]])
180 } 223 }
182 # i=class2names[j] 225 # i=class2names[j]
183 tic <- TIC[[classnames[[l]][j]]] 226 tic <- TIC[[classnames[[l]][j]]]
184 points(tic[,1]/60, -tic[,2], col = cols[classnames[[l]][j]], pch = pch[classnames[[l]][j]], type="l") 227 points(tic[,1]/60, -tic[,2], col = cols[classnames[[l]][j]], pch = pch[classnames[[l]][j]], type="l")
185 colvect<-append(colvect,cols[classnames[[l]][j]]) 228 colvect<-append(colvect,cols[classnames[[l]][j]])
186 } 229 }
187 legend("topright",paste(basename(files[c(classnames[[k]],classnames[[l]])])), col = colvect, lty = lty, pch = pch) 230 legend("topright",paste(gsub("(^.+)\\..*$","\\1",basename(files[c(classnames[[k]],classnames[[l]])]))), col = colvect, lty = lty, pch = pch)
188 } 231 }
189 } 232 }
190 }#end if length >2 233 }#end if length >2
234
191 if (length(class)==2){ 235 if (length(class)==2){
192
193
194 k=1 236 k=1
195 l=2 237 l=2
196
197 plot(0, 0, type="n", xlim = xlim/60, ylim = ylim, main = paste("Total Ion Chromatograms \n","TICs_",class[k],"vs",class[l], sep=""), xlab = "Retention Time (min)", ylab = "TIC")
198 colvect<-NULL 238 colvect<-NULL
239 plot(0, 0, type="n", xlim = xlim/60, ylim = ylim, main = paste("Base Peak Chromatograms \n","BPCs_",class[k],"vs",class[l], sep=""), xlab = "Retention Time (min)", ylab = "BPC")
240
199 for (j in 1:length(classnames[[k]])) { 241 for (j in 1:length(classnames[[k]])) {
200
201 tic <- TIC[[classnames[[k]][j]]] 242 tic <- TIC[[classnames[[k]][j]]]
202 # points(tic[,1]/60, tic[,2], col = cols[i], pch = pch[i], type="l") 243 # points(tic[,1]/60, tic[,2], col = cols[i], pch = pch[i], type="l")
203 points(tic[,1]/60, tic[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l") 244 points(tic[,1]/60, tic[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l")
204 colvect<-append(colvect,cols[classnames[[k]][j]]) 245 colvect<-append(colvect,cols[classnames[[k]][j]])
205 } 246 }
207 # i=class2names[j] 248 # i=class2names[j]
208 tic <- TIC[[classnames[[l]][j]]] 249 tic <- TIC[[classnames[[l]][j]]]
209 points(tic[,1]/60, -tic[,2], col = cols[classnames[[l]][j]], pch = pch[classnames[[l]][j]], type="l") 250 points(tic[,1]/60, -tic[,2], col = cols[classnames[[l]][j]], pch = pch[classnames[[l]][j]], type="l")
210 colvect<-append(colvect,cols[classnames[[l]][j]]) 251 colvect<-append(colvect,cols[classnames[[l]][j]])
211 } 252 }
212 legend("topright",paste(basename(files[c(classnames[[k]],classnames[[l]])])), col = colvect, lty = lty, pch = pch) 253 legend("topright",paste(gsub("(^.+)\\..*$","\\1",basename(files[c(classnames[[k]],classnames[[l]])]))), col = colvect, lty = lty, pch = pch)
213
214 }#end length ==2 254 }#end length ==2
255
215 if (length(class)==1){ 256 if (length(class)==1){
216 k=1 257 k=1
217 ylim = range(sapply(TIC, function(x) range(x[,2]))) 258 ylim = range(sapply(TIC, function(x) range(x[,2])))
218
219 plot(0, 0, type="n", xlim = xlim/60, ylim = ylim, main = paste("Total Ion Chromatograms \n","TICs_",class[k], sep=""), xlab = "Retention Time (min)", ylab = "TIC")
220 colvect<-NULL 259 colvect<-NULL
260 plot(0, 0, type="n", xlim = xlim/60, ylim = ylim, main = paste("Base Peak Chromatograms \n","BPCs_",class[k], sep=""), xlab = "Retention Time (min)", ylab = "BPC")
261
221 for (j in 1:length(classnames[[k]])) { 262 for (j in 1:length(classnames[[k]])) {
222
223 tic <- TIC[[classnames[[k]][j]]] 263 tic <- TIC[[classnames[[k]][j]]]
224 # points(tic[,1]/60, tic[,2], col = cols[i], pch = pch[i], type="l") 264 # points(tic[,1]/60, tic[,2], col = cols[i], pch = pch[i], type="l")
225 points(tic[,1]/60, tic[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l") 265 points(tic[,1]/60, tic[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l")
226 colvect<-append(colvect,cols[classnames[[k]][j]]) 266 colvect<-append(colvect,cols[classnames[[k]][j]])
227 } 267 }
228 268 legend("topright",paste(gsub("(^.+)\\..*$","\\1",basename(files[c(classnames[[k]])]))), col = colvect, lty = lty, pch = pch)
229 legend("topright",paste(basename(files[c(classnames[[k]])])), col = colvect, lty = lty, pch = pch)
230
231 }#end length ==1 269 }#end length ==1
232 dev.off() 270 dev.off()
233 271 }
234 # invisible(TIC) 272
235 } 273 getTIC <- function(file,rtcor=NULL) {
236 274 object <- xcmsRaw(file)
237 275 cbind(if (is.null(rtcor)) object@scantime else rtcor, rawEIC(object,mzrange=range(object@env$mz))$intensity)
276 }
277
278 ## overlay TIC from all files in current folder or from xcmsSet, create pdf
279 getTIC2s <- function(files, xset=NULL, pdfname="TICs.pdf", rt=c("raw","corrected")) {
280 require(xcms)
281
282 #create sampleMetadata, get sampleMetadata and class
283 if(!is.null(xset)){
284 #When files come from XCMS3 before metaMS treatment
285 sampleMetadata<-xset@phenoData
286 } else {
287 #When files come from a zip directory with raw files before metaMS
288 sampleMetadata<-xcms:::phenoDataFromPaths(files)
289 }
290 class<-as.vector(levels(sampleMetadata[,"class"])) #create phenoData like table
291 classnames<-vector("list",length(class))
292 for (i in 1:length(class)){
293 classnames[[i]]<-which( sampleMetadata[,"class"]==class[i])
294 }
295
296 N <- dim(sampleMetadata)[1]
297 TIC <- vector("list",N)
298
299 for (i in 1:N) {
300 cat(files[i],"\n")
301 if (!is.null(xcmsSet) && rt == "corrected")
302 rtcor <- xcmsSet@rt$corrected[[i]]
303 else
304 rtcor <- NULL
305 TIC[[i]] <- getTIC(files[i],rtcor=rtcor)
306 }
307
308 pdf(pdfname,w=16,h=10)
309 cols <- rainbow(N)
310 lty = 1:N
311 pch = 1:N
312 #search for max x and max y in TICs
313 xlim = range(sapply(TIC, function(x) range(x[,1])))
314 ylim = range(sapply(TIC, function(x) range(x[,2])))
315 ylim = c(-ylim[2], ylim[2])
316
317 ##plot start
318 if (length(class)>2){
319 for (k in 1:(length(class)-1)){
320 for (l in (k+1):length(class)){
321 print(paste(class[k],"vs",class[l],sep=" "))
322 plot(0, 0, type="n", xlim = xlim/60, ylim = ylim, main = paste("Total Ion Chromatograms \n","TICs_",class[k]," vs ",class[l], sep=""), xlab = "Retention Time (min)", ylab = "TIC")
323 colvect<-NULL
324 for (j in 1:length(classnames[[k]])) {
325 tic <- TIC[[classnames[[k]][j]]]
326 # points(tic[,1]/60, tic[,2], col = cols[i], pch = pch[i], type="l")
327 points(tic[,1]/60, tic[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l")
328 colvect<-append(colvect,cols[classnames[[k]][j]])
329 }
330 for (j in 1:length(classnames[[l]])) {
331 # i=class2names[j]
332 tic <- TIC[[classnames[[l]][j]]]
333 points(tic[,1]/60, -tic[,2], col = cols[classnames[[l]][j]], pch = pch[classnames[[l]][j]], type="l")
334 colvect<-append(colvect,cols[classnames[[l]][j]])
335 }
336 legend("topright",paste(gsub("(^.+)\\..*$","\\1",basename(files[c(classnames[[k]],classnames[[l]])]))), col = colvect, lty = lty, pch = pch)
337 }
338 }
339 }#end if length >2
340
341 if (length(class)==2){
342 k=1
343 l=2
344 plot(0, 0, type="n", xlim = xlim/60, ylim = ylim, main = paste("Total Ion Chromatograms \n","TICs_",class[k],"vs",class[l], sep=""), xlab = "Retention Time (min)", ylab = "TIC")
345 colvect<-NULL
346 for (j in 1:length(classnames[[k]])) {
347 tic <- TIC[[classnames[[k]][j]]]
348 # points(tic[,1]/60, tic[,2], col = cols[i], pch = pch[i], type="l")
349 points(tic[,1]/60, tic[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l")
350 colvect<-append(colvect,cols[classnames[[k]][j]])
351 }
352 for (j in 1:length(classnames[[l]])) {
353 # i=class2names[j]
354 tic <- TIC[[classnames[[l]][j]]]
355 points(tic[,1]/60, -tic[,2], col = cols[classnames[[l]][j]], pch = pch[classnames[[l]][j]], type="l")
356 colvect<-append(colvect,cols[classnames[[l]][j]])
357 }
358 legend("topright",paste(gsub("(^.+)\\..*$","\\1",basename(files[c(classnames[[k]],classnames[[l]])]))), col = colvect, lty = lty, pch = pch)
359 }#end length ==2
360
361 if (length(class)==1){
362 k=1
363 ylim = range(sapply(TIC, function(x) range(x[,2])))
364 plot(0, 0, type="n", xlim = xlim/60, ylim = ylim, main = paste("Total Ion Chromatograms \n","TICs_",class[k], sep=""), xlab = "Retention Time (min)", ylab = "TIC")
365 colvect<-NULL
366 for (j in 1:length(classnames[[k]])) {
367 tic <- TIC[[classnames[[k]][j]]]
368 # points(tic[,1]/60, tic[,2], col = cols[i], pch = pch[i], type="l")
369 points(tic[,1]/60, tic[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l")
370 colvect<-append(colvect,cols[classnames[[k]][j]])
371 }
372 legend("topright",paste(gsub("(^.+)\\..*$","\\1",basename(files[c(classnames[[k]])]))), col = colvect, lty = lty, pch = pch)
373 }#end length ==1
374 dev.off()
375 }
376
377
378 #Update by J.Saint-Vanne
238 ##addition for quality control of peak picking 379 ##addition for quality control of peak picking
239 #metaMS EIC and pspectra plotting option 380 #metaMS EIC and pspectra plotting option
240 #version 20150512 381 #version 20190520
241 #only for Galaxy 382 #only for Galaxy
242 383 plotUnknowns<-function(resGC, unkn="", DB=NULL, fileFrom=NULL){
243 plotUnknowns<-function(resGC, unkn=""){
244
245 384
246 ##Annotation table each value is a pcgrp associated to the unknown 385 ##Annotation table each value is a pcgrp associated to the unknown
247 ##NOTE pcgrp index are different between xcmsSet and resGC due to filtering steps in metaMS 386 ##NOTE pcgrp index are different between xcmsSet and resGC due to filtering steps in metaMS
248 ##R. Wehrens give me some clues on that and we found a correction 387 ##R. Wehrens give me some clues on that and we found a correction
249 #if unkn="none" 388
250 389 #correction of annotation matrix due to pcgrp removal by quality check in runGCresult
251 if(unkn=="none") { 390 #matrix of correspondance between an@pspectra and filtered pspectra from runGC
252 pdf("Unknown_Empty.pdf") 391 #Select only pspectra which correpond to them select in metaMS
253 plot.new() 392 # col1 = filtered spectra from runGC and col2 = an@spectra
254 text(x=0.5,y=1,pos=1, labels="No EIC ploting required") 393 allPCGRPs <-lapply(1:length(resGC$xset),
255 dev.off()
256 }else {
257
258 mat<-matrix(ncol=length(resGC$xset), nrow=dim(resGC$PeakTable)[1])
259
260 for (j in 1: length(resGC$xset)){
261 test<-resGC$annotation[[j]]
262 print(paste("j=",j))
263 for (i in 1:dim(test)[1]){
264 if (as.numeric(row.names(test)[i])>dim(mat)[1]){
265 next
266 } else {
267 mat[as.numeric(row.names(test)[i]),j]<-test[i,1]
268 }
269 }
270 }
271 colnames(mat)<-colnames(resGC$PeakTable[,c((which(colnames(resGC$PeakTable)=="rt"|colnames(resGC$PeakTable)=="RI")[length(which(colnames(resGC$PeakTable)=="rt"|colnames(resGC$PeakTable)=="RI"))]+1):dim(resGC$PeakTable)[2])])
272
273 #debug
274
275 # print(dim(mat))
276 # print(mat[1:3,])
277 # write.table(mat, file="myannotationtable.tsv", sep="\t", row.names=FALSE)
278 #correction of annotation matrix due to pcgrp removal by quality check in runGCresult
279 #matrix of correspondance between an@pspectra and filtered pspectra from runGC
280
281 allPCGRPs <-
282 lapply(1:length(resGC$xset),
283 function(i) { 394 function(i) {
284 an <- resGC$xset[[i]] 395 an <- resGC$xset[[i]]
285 huhn <- an@pspectra[which(sapply(an@pspectra, length) >= 396 huhn <- an@pspectra[which(sapply(an@pspectra, length) >=
286 metaSetting(resGC$settings, 397 metaSetting(resGC$settings,
287 "DBconstruction.minfeat"))] 398 "DBconstruction.minfeat"))]
288 matCORR<-cbind(1:length(huhn), match(huhn, an@pspectra)) 399 matCORR<-cbind(1:length(huhn), match(huhn, an@pspectra))
289 }) 400 })
290 401 #Build a new annotation list with sampnames and pseudospectra number from xset
291 if (unkn[1]==""){ 402 helpannotation <- list()
292 #plot EIC and spectra for all unknown for comparative purpose 403 for(j in 1:length(resGC$xset)){
293 404 helpannotation[[j]] <- resGC$annotation[[j]][1:2]
294 405 pspvector <- vector()
295 par (mar=c(5, 4, 4, 2) + 0.1) 406 for(i in 1: nrow(helpannotation[[j]])){
296 for (l in 1:dim(resGC$PeakTable)[1]){ #l=2 407 #Find corresponding pspec
297 #recordPlot 408 psplink <- allPCGRPs[[j]][match(helpannotation[[j]][i,1],allPCGRPs[[j]]),2]
298 perpage=3 #if change change layout also! 409 pspvector <- c(pspvector,psplink)
299 num.plots <- ceiling(dim(mat)[2]/perpage) #three pcgroup per page 410 #Change the annotation column into sampname column
300 my.plots <- vector(num.plots, mode='list') 411 if(helpannotation[[j]][i,2] < 0){
301 dev.new(width=21/2.54, height=29.7/2.54, file=paste("Unknown_",l,".pdf", sep="")) #A4 pdf 412 #It's an unknown
302 # par(mfrow=c(perpage,2)) 413 new_name <- paste("Unknown",abs(as.integer(helpannotation[[j]][i,2])))
303 layout(matrix(c(1,1,2,3,4,4,5,6,7,7,8,9), 6, 2, byrow = TRUE), widths=rep(c(1,1),perpage), heights=rep(c(1,5),perpage)) 414 helpannotation[[j]][i,2] <- new_name
304 # layout.show(6) 415 }else{
305 oma.saved <- par("oma") 416 #It has been found in local database
306 par(oma = rep.int(0, 4)) 417 for(k in 1:length(DB)){
307 par(oma = oma.saved) 418 if(helpannotation[[j]][i,2] == k){
308 o.par <- par(mar = rep.int(0, 4)) 419 helpannotation[[j]][i,2] <- DB[[k]]$Name
309 on.exit(par(o.par)) 420 break
310 stop=0 #initialize 421 }
311 for (i in 1:num.plots) { 422 }
312 start=stop+1 423 }
313 stop=start+perpage-1 # 424 }
314 for (c in start:stop){ 425 helpannotation[[j]] <- cbind(helpannotation[[j]],pspvector)
315 if (c <=dim(mat)[2]){ 426 names(helpannotation)[j] <- names(resGC$annotation[j])
316 427 }
317 #get sample name 428 peaktable <- resGC$PeakTable
318 sampname<-basename(resGC$xset[[c]]@xcmsSet@filepaths) 429
319 430 par (mar=c(5, 4, 4, 2) + 0.1)
320 #remove .cdf, .mzXML filepattern 431 #For each unknown
321 filepattern <- c("[Cc][Dd][Ff]", "[Nn][Cc]", "([Mm][Zz])?[Xx][Mm][Ll]", 432 for (l in 1:length(unkn)){
433 #recordPlot
434 perpage=3 #if change change layout also!
435 dev.new(width=21/2.54, height=29.7/2.54, file=paste("Unknown_",unkn[l],".pdf", sep="")) #A4 pdf
436 # par(mfrow=c(perpage,2))
437 layout(matrix(c(1,1,2,3,4,4,5,6,7,7,8,9), 6, 2, byrow = TRUE), widths=rep(c(1,1),perpage), heights=rep(c(1,5),perpage))
438 # layout.show(6)
439 oma.saved <- par("oma")
440 par(oma = rep.int(0, 4))
441 par(oma = oma.saved)
442 o.par <- par(mar = rep.int(0, 4))
443 on.exit(par(o.par))
444 #For each sample
445 for (c in 1:length(resGC$xset)) {
446 #get sample name
447 sampname<-basename(resGC$xset[[c]]@xcmsSet@filepaths)
448 #remove .cdf, .mzXML filepattern
449 filepattern <- c("[Cc][Dd][Ff]", "[Nn][Cc]", "([Mm][Zz])?[Xx][Mm][Ll]",
322 "[Mm][Zz][Dd][Aa][Tt][Aa]", "[Mm][Zz][Mm][Ll]") 450 "[Mm][Zz][Dd][Aa][Tt][Aa]", "[Mm][Zz][Mm][Ll]")
323 filepattern <- paste(paste("\\.", filepattern, "$", sep = ""), 451 filepattern <- paste(paste("\\.", filepattern, "$", sep = ""), collapse = "|")
324 collapse = "|") 452 sampname<-gsub(filepattern, "",sampname)
325 sampname<-gsub(filepattern, "",sampname) 453 title1<-paste(peaktable[unkn[l],1],"from",sampname, sep = " ")
326 454 an<-resGC$xset[[c]]
327 title1<-paste("unknown", l,"from",sampname, sep=" ") 455 if(fileFrom == "zipfile") {
328 an<-resGC$xset[[c]] 456 an@xcmsSet@filepaths <- paste0("./",an@xcmsSet@phenoData[,"class"],"/",basename(an@xcmsSet@filepaths))
329 457 }#else {
458 #print(an@xcmsSet@filepaths)
459 #an@xcmsSet@filepaths <- paste0("./",basename(an@xcmsSet@filepaths))
460 #}
461 #Find the good annotation for this sample
462 for(a in 1:length(helpannotation)){
463 if(gsub(filepattern, "", names(helpannotation)[a]) == paste0("./",sampname)){
464 #Find the unkn or the matched std in this sample
465 findunkn <- FALSE
466 for(r in 1:nrow(helpannotation[[a]])){
467 if(helpannotation[[a]][r,"annotation"] == peaktable[unkn[l],1]){
468 findunkn <- TRUE
469 pcgrp <- helpannotation[[a]][r,"pspvector"]
330 par (mar=c(0, 0, 0, 0) + 0.1) 470 par (mar=c(0, 0, 0, 0) + 0.1)
471 #Write title
331 plot.new() 472 plot.new()
332 box() 473 box()
333 text(0.5, 0.5, title1, cex=2) 474 text(0.5, 0.5, title1, cex=2)
334 if (!is.na(mat[l,c])){ 475 par (mar=c(3, 2.5, 3, 1.5) + 0.1)
335 pcgrp=allPCGRPs[[c]][which(allPCGRPs[[c]][,1]==mat[l,c]),2] 476 #Window for EIC
336 if (pcgrp!=mat[l,c]) print ("pcgrp changed") 477 plotEICs(an, pspec=pcgrp, maxlabel=2)
337 par (mar=c(3, 2.5, 3, 1.5) + 0.1) 478 #Window for pseudospectra
338 plotEICs(an, pspec=pcgrp, maxlabel=2) 479 plotPsSpectrum(an, pspec=pcgrp, maxlabel=2)
339 plotPsSpectrum(an, pspec=pcgrp, maxlabel=2)
340 } else {
341 plot.new()
342 box()
343 text(0.5, 0.5, "NOT FOUND", cex=2)
344 plot.new()
345 box()
346 text(0.5, 0.5, "NOT FOUND", cex=2)
347 }
348 }
349 }
350 # my.plots[[i]] <- recordPlot()
351 }
352 graphics.off()
353
354 # pdf(file=paste("Unknown_",l,".pdf", sep=""), onefile=TRUE)
355 # for (my.plot in my.plots) {
356 # replayPlot(my.plot)
357 # }
358 # my.plots
359 # graphics.off()
360
361 }#end for l
362 }#end if unkn=""
363 else{
364 par (mar=c(5, 4, 4, 2) + 0.1)
365 l=unkn
366 if (length(l)==1){
367 #recordPlot
368 perpage=3 #if change change layout also!
369 num.plots <- ceiling(dim(mat)[2]/perpage) #three pcgroup per page
370 my.plots <- vector(num.plots, mode='list')
371
372 dev.new(width=21/2.54, height=29.7/2.54, file=paste("Unknown_",l,".pdf", sep="")) #A4 pdf
373 # par(mfrow=c(perpage,2))
374 layout(matrix(c(1,1,2,3,4,4,5,6,7,7,8,9), 6, 2, byrow = TRUE), widths=rep(c(1,1),perpage), heights=rep(c(1,5),perpage))
375 # layout.show(6)
376 oma.saved <- par("oma")
377 par(oma = rep.int(0, 4))
378 par(oma = oma.saved)
379 o.par <- par(mar = rep.int(0, 4))
380 on.exit(par(o.par))
381 stop=0 #initialize
382 for (i in 1:num.plots) {
383 start=stop+1
384 stop=start+perpage-1 #
385 for (c in start:stop){
386 if (c <=dim(mat)[2]){
387
388 #get sample name
389 sampname<-basename(resGC$xset[[c]]@xcmsSet@filepaths)
390
391 #remove .cdf, .mzXML filepattern
392 filepattern <- c("[Cc][Dd][Ff]", "[Nn][Cc]", "([Mm][Zz])?[Xx][Mm][Ll]", "[Mm][Zz][Dd][Aa][Tt][Aa]", "[Mm][Zz][Mm][Ll]")
393 filepattern <- paste(paste("\\.", filepattern, "$", sep = ""), collapse = "|")
394 sampname<-gsub(filepattern, "",sampname)
395
396 title1<-paste("unknown", l,"from",sampname, sep=" ")
397 an<-resGC$xset[[c]]
398
399 par (mar=c(0, 0, 0, 0) + 0.1)
400 plot.new()
401 box()
402 text(0.5, 0.5, title1, cex=2)
403 if (!is.na(mat[l,c])){
404 pcgrp=allPCGRPs[[c]][which(allPCGRPs[[c]][,1]==mat[l,c]),2]
405 if (pcgrp!=mat[l,c]) print ("pcgrp changed")
406 par (mar=c(3, 2.5, 3, 1.5) + 0.1)
407 plotEICs(an, pspec=pcgrp, maxlabel=2)
408 plotPsSpectrum(an, pspec=pcgrp, maxlabel=2)
409 } else {
410 plot.new()
411 box()
412 text(0.5, 0.5, "NOT FOUND", cex=2)
413 plot.new()
414 box()
415 text(0.5, 0.5, "NOT FOUND", cex=2)
416 }
417 } 480 }
418 } 481 }
419 # my.plots[[i]] <- recordPlot() 482 if(!findunkn) {
483 par (mar=c(0, 0, 0, 0) + 0.1)
484 #Write title
485 plot.new()
486 box()
487 text(0.5, 0.5, title1, cex=2)
488 #Window for EIC
489 plot.new()
490 box()
491 text(0.5, 0.5, "NOT FOUND", cex=2)
492 #Window for pseudospectra
493 plot.new()
494 box()
495 text(0.5, 0.5, "NOT FOUND", cex=2)
496 }
497 break
420 } 498 }
421 graphics.off() 499 }
422
423 # pdf(file=paste("Unknown_",l,".pdf", sep=""), onefile=TRUE)
424 # for (my.plot in my.plots) {
425 # replayPlot(my.plot)
426 # }
427 # my.plots
428 # graphics.off()
429 } else {
430 par (mar=c(5, 4, 4, 2) + 0.1)
431 for (l in 1:length(unkn)){ #l=2
432 #recordPlot
433 perpage=3 #if change change layout also!
434 num.plots <- ceiling(dim(mat)[2]/perpage) #three pcgroup per page
435 my.plots <- vector(num.plots, mode='list')
436 dev.new(width=21/2.54, height=29.7/2.54, file=paste("Unknown_",unkn[l],".pdf", sep="")) #A4 pdf
437 # par(mfrow=c(perpage,2))
438 layout(matrix(c(1,1,2,3,4,4,5,6,7,7,8,9), 6, 2, byrow = TRUE), widths=rep(c(1,1),perpage), heights=rep(c(1,5),perpage))
439 # layout.show(6)
440 oma.saved <- par("oma")
441 par(oma = rep.int(0, 4))
442 par(oma = oma.saved)
443 o.par <- par(mar = rep.int(0, 4))
444 on.exit(par(o.par))
445 stop=0 #initialize
446 for (i in 1:num.plots) {
447 start=stop+1
448 stop=start+perpage-1 #
449 for (c in start:stop){
450 if (c <=dim(mat)[2]){
451
452 #get sample name
453 sampname<-basename(resGC$xset[[c]]@xcmsSet@filepaths)
454
455 #remove .cdf, .mzXML filepattern
456 filepattern <- c("[Cc][Dd][Ff]", "[Nn][Cc]", "([Mm][Zz])?[Xx][Mm][Ll]", "[Mm][Zz][Dd][Aa][Tt][Aa]", "[Mm][Zz][Mm][Ll]")
457 filepattern <- paste(paste("\\.", filepattern, "$", sep = ""), collapse = "|")
458 sampname<-gsub(filepattern, "",sampname)
459
460 title1<-paste("unknown",unkn[l],"from",sampname, sep=" ")
461 an<-resGC$xset[[c]]
462
463 par (mar=c(0, 0, 0, 0) + 0.1)
464 plot.new()
465 box()
466 text(0.5, 0.5, title1, cex=2)
467 if (!is.na(mat[unkn[l],c])){
468 pcgrp=allPCGRPs[[c]][which(allPCGRPs[[c]][,1]==mat[unkn[l],c]),2]
469 if (pcgrp!=mat[unkn[l],c]) print ("pcgrp changed")
470 par (mar=c(3, 2.5, 3, 1.5) + 0.1)
471 plotEICs(an, pspec=pcgrp, maxlabel=2)
472 plotPsSpectrum(an, pspec=pcgrp, maxlabel=2)
473 } else {
474 plot.new()
475 box()
476 text(0.5, 0.5, "NOT FOUND", cex=2)
477 plot.new()
478 box()
479 text(0.5, 0.5, "NOT FOUND", cex=2)
480 }
481 }
482 }
483 # my.plots[[i]] <- recordPlot()
484 }
485 graphics.off()
486
487 # pdf(file=paste("Unknown_",unkn[l],".pdf", sep=""), onefile=TRUE)
488 # for (my.plot in my.plots) {
489 # replayPlot(my.plot)
490 # }
491 # my.plots
492 # graphics.off()
493
494 }#end for unkn[l]
495
496 }
497
498 } 500 }
499 } 501 graphics.off()
500 } #end function 502 }#end for unkn[l]
501 503 }#end function
502 # This function get the raw file path from the arguments
503 getRawfilePathFromArguments <- function(singlefile, zipfile, listArguments) {
504 if (!is.null(listArguments[["zipfile"]])) zipfile = listArguments[["zipfile"]]
505 if (!is.null(listArguments[["zipfilePositive"]])) zipfile = listArguments[["zipfilePositive"]]
506 if (!is.null(listArguments[["zipfileNegative"]])) zipfile = listArguments[["zipfileNegative"]]
507
508 if (!is.null(listArguments[["singlefile_galaxyPath"]])) {
509 singlefile_galaxyPaths = listArguments[["singlefile_galaxyPath"]];
510 singlefile_sampleNames = listArguments[["singlefile_sampleName"]]
511 }
512 if (!is.null(listArguments[["singlefile_galaxyPathPositive"]])) {
513 singlefile_galaxyPaths = listArguments[["singlefile_galaxyPathPositive"]];
514 singlefile_sampleNames = listArguments[["singlefile_sampleNamePositive"]]
515 }
516 if (!is.null(listArguments[["singlefile_galaxyPathNegative"]])) {
517 singlefile_galaxyPaths = listArguments[["singlefile_galaxyPathNegative"]];
518 singlefile_sampleNames = listArguments[["singlefile_sampleNameNegative"]]
519 }
520 if (exists("singlefile_galaxyPaths")){
521 singlefile_galaxyPaths = unlist(strsplit(singlefile_galaxyPaths,","))
522 singlefile_sampleNames = unlist(strsplit(singlefile_sampleNames,","))
523
524 singlefile=NULL
525 for (singlefile_galaxyPath_i in seq(1:length(singlefile_galaxyPaths))) {
526 singlefile_galaxyPath=singlefile_galaxyPaths[singlefile_galaxyPath_i]
527 singlefile_sampleName=singlefile_sampleNames[singlefile_galaxyPath_i]
528 singlefile[[singlefile_sampleName]] = singlefile_galaxyPath
529 }
530 }
531 for (argument in c("zipfile","zipfilePositive","zipfileNegative","singlefile_galaxyPath","singlefile_sampleName","singlefile_galaxyPathPositive","singlefile_sampleNamePositive","singlefile_galaxyPathNegative","singlefile_sampleNameNegative")) {
532 listArguments[[argument]]=NULL
533 }
534 return(list(zipfile=zipfile, singlefile=singlefile, listArguments=listArguments))
535 }
536
537
538 # This function retrieve the raw file in the working directory
539 # - if zipfile: unzip the file with its directory tree
540 # - if singlefiles: set symlink with the good filename
541 retrieveRawfileInTheWorkingDirectory <- function(singlefile, zipfile) {
542 if(!is.null(singlefile) && (length("singlefile")>0)) {
543 for (singlefile_sampleName in names(singlefile)) {
544 singlefile_galaxyPath = singlefile[[singlefile_sampleName]]
545 if(!file.exists(singlefile_galaxyPath)){
546 error_message=paste("Cannot access the sample:",singlefile_sampleName,"located:",singlefile_galaxyPath,". Please, contact your administrator ... if you have one!")
547 print(error_message); stop(error_message)
548 }
549
550 file.symlink(singlefile_galaxyPath,singlefile_sampleName)
551 }
552 directory = "."
553
554 }
555 if(!is.null(zipfile) && (zipfile!="")) {
556 if(!file.exists(zipfile)){
557 error_message=paste("Cannot access the Zip file:",zipfile,". Please, contact your administrator ... if you have one!")
558 print(error_message)
559 stop(error_message)
560 }
561
562 #list all file in the zip file
563 #zip_files=unzip(zipfile,list=T)[,"Name"]
564
565 #unzip
566 suppressWarnings(unzip(zipfile, unzip="unzip"))
567
568 #get the directory name
569 filesInZip=unzip(zipfile, list=T);
570 directories=unique(unlist(lapply(strsplit(filesInZip$Name,"/"), function(x) x[1])));
571 directories=directories[!(directories %in% c("__MACOSX")) & file.info(directories)$isdir]
572 directory = "."
573 if (length(directories) == 1) directory = directories
574
575 cat("files_root_directory\t",directory,"\n")
576
577 }
578 return (directory)
579 }
580