Mercurial > repos > yguitton > metams_rungc
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 |