Mercurial > repos > workflow4metabolomics > metams_plot

comparison lib_metams.r @ 2:96debae917e4 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/workflow4metabolomics/metaMS commit c7a518686137f6d62b7415715152e8d5a9953ed7
author workflow4metabolomics
date Fri, 06 Sep 2019 06:12:08 -0400
parents 708ab9928a70
children
comparison
equal deleted inserted replaced
1:708ab9928a70 2:96debae917e4
160 } 160 }
161 161
162 ##ADDITIONS FROM Y. Guitton 162 ##ADDITIONS FROM Y. Guitton
163 getBPC <- function(file,rtcor=NULL, ...) { 163 getBPC <- function(file,rtcor=NULL, ...) {
164 object <- xcmsRaw(file) 164 object <- xcmsRaw(file)
165 sel <- profRange(object, ...) 165 sel <- profRange(object, ...)
166 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]))
167 } 167 }
168 168
169 getBPC2s <- function (files, xset = NULL, pdfname="BPCs.pdf", rt = c("raw","corrected"), scanrange=NULL) { 169 getBPC2s <- function (files, xset = NULL, pdfname="BPCs.pdf", rt = c("raw","corrected"), scanrange=NULL) {
170 require(xcms) 170 require(xcms)
171
172 #Verification for cdf files
173 stop=FALSE
174 for(i in 1:length(files)){
175 extension <- unlist(strsplit(basename(files[i]),"\\."))[length(unlist(strsplit(basename(files[i]),"\\.")))]
176 if(extension == "CDF" || extension == "cdf"){
177 stop = TRUE
178 break
179 }
180 }
181 if(stop){
182 error_message <- "You have a CDF file and there is an issue to resolve on them for chromatograms.... !"
183 print(error_message)
184 stop(error_message)
185 }
186 171
187 #create sampleMetadata, get sampleMetadata and class 172 #create sampleMetadata, get sampleMetadata and class
188 if(!is.null(xset)) { 173 if(!is.null(xset)) {
189 #When files come from XCMS3 directly before metaMS 174 #When files come from XCMS3 directly before metaMS
190 sampleMetadata <- xset@phenoData 175 sampleMetadata <- xset@phenoData
197 for (i in 1:length(class)){ 182 for (i in 1:length(class)){
198 classnames[[i]]<-which( sampleMetadata[,"class"]==class[i]) 183 classnames[[i]]<-which( sampleMetadata[,"class"]==class[i])
199 } 184 }
200 185
201 N <- dim(sampleMetadata)[1] 186 N <- dim(sampleMetadata)[1]
202 TIC <- vector("list",N) 187 BPC <- vector("list",N)
203 188
204 for (j in 1:N) { 189 for (j in 1:N) {
205 TIC[[j]] <- getBPC(files[j]) 190 BPC[[j]] <- getBPC(files[j])
206 #good for raw 191 #good for raw
207 # seems strange for corrected 192 # seems strange for corrected
208 #errors if scanrange used in xcmsSetgeneration 193 #errors if scanrange used in xcmsSetgeneration
209 if (!is.null(xcmsSet) && rt == "corrected"){ 194 if (!is.null(xcmsSet) && rt == "corrected"){
210 rtcor <- xcmsSet@rt$corrected[[j]] 195 rtcor <- xcmsSet@rt$corrected[[j]]
211 }else{ 196 }else{
212 rtcor <- NULL 197 rtcor <- NULL
213 } 198 }
214 TIC[[j]] <- getBPC(files[j],rtcor=rtcor) 199 BPC[[j]] <- getBPC(files[j],rtcor=rtcor)
215 } 200 }
216 201
217 pdf(pdfname,w=16,h=10) 202 pdf(pdfname,w=16,h=10)
218 cols <- rainbow(N) 203 cols <- rainbow(N)
219 lty = 1:N 204 lty = 1:N
220 pch = 1:N 205 pch = 1:N
221 #search for max x and max y in BPCs 206 #search for max x and max y in BPCs
222 xlim = range(sapply(TIC, function(x) range(x[,1]))) 207
223 ylim = range(sapply(TIC, function(x) range(x[,2]))) 208 xlim = range(sapply(BPC, function(x) range(x[,1])))
209 ylim = range(sapply(BPC, function(x) range(x[,2])))
210
224 ylim = c(-ylim[2], ylim[2]) 211 ylim = c(-ylim[2], ylim[2])
225 212
226 ##plot start 213 ##plot start
227 if (length(class)>2){ 214 if (length(class)>2){
228 for (k in 1:(length(class)-1)){ 215 for (k in 1:(length(class)-1)){
229 for (l in (k+1):length(class)){ 216 for (l in (k+1):length(class)){
230 cat(paste(class[k],"vs",class[l],sep=" ","\n")) 217 cat(paste(class[k],"vs",class[l],sep=" ","\n"))
231 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") 218 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")
232 colvect<-NULL 219 colvect<-NULL
233 for (j in 1:length(classnames[[k]])) { 220 for (j in 1:length(classnames[[k]])) {
234 tic <- TIC[[classnames[[k]][j]]] 221 bpc <- BPC[[classnames[[k]][j]]]
235 # points(tic[,1]/60, tic[,2], col = cols[i], pch = pch[i], type="l") 222 # points(bpc[,1]/60, bpc[,2], col = cols[i], pch = pch[i], type="l")
236 points(tic[,1]/60, tic[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l") 223 points(bpc[,1]/60, bpc[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l")
237 colvect<-append(colvect,cols[classnames[[k]][j]]) 224 colvect<-append(colvect,cols[classnames[[k]][j]])
238 } 225 }
239 for (j in 1:length(classnames[[l]])) { 226 for (j in 1:length(classnames[[l]])) {
240 # i=class2names[j] 227 # i=class2names[j]
241 tic <- TIC[[classnames[[l]][j]]] 228 bpc <- BPC[[classnames[[l]][j]]]
242 points(tic[,1]/60, -tic[,2], col = cols[classnames[[l]][j]], pch = pch[classnames[[l]][j]], type="l") 229 points(bpc[,1]/60, -bpc[,2], col = cols[classnames[[l]][j]], pch = pch[classnames[[l]][j]], type="l")
243 colvect<-append(colvect,cols[classnames[[l]][j]]) 230 colvect<-append(colvect,cols[classnames[[l]][j]])
244 } 231 }
245 legend("topright",paste(gsub("(^.+)\\..*$","\\1",basename(files[c(classnames[[k]],classnames[[l]])]))), col = colvect, lty = lty, pch = pch) 232 legend("topright",paste(gsub("(^.+)\\..*$","\\1",basename(files[c(classnames[[k]],classnames[[l]])]))), col = colvect, lty = lty, pch = pch)
246 } 233 }
247 } 234 }
252 l=2 239 l=2
253 colvect<-NULL 240 colvect<-NULL
254 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") 241 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")
255 242
256 for (j in 1:length(classnames[[k]])) { 243 for (j in 1:length(classnames[[k]])) {
257 tic <- TIC[[classnames[[k]][j]]] 244 bpc <- BPC[[classnames[[k]][j]]]
258 # points(tic[,1]/60, tic[,2], col = cols[i], pch = pch[i], type="l") 245 # points(bpc[,1]/60, bpc[,2], col = cols[i], pch = pch[i], type="l")
259 points(tic[,1]/60, tic[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l") 246 points(bpc[,1]/60, bpc[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l")
260 colvect<-append(colvect,cols[classnames[[k]][j]]) 247 colvect<-append(colvect,cols[classnames[[k]][j]])
261 } 248 }
262 for (j in 1:length(classnames[[l]])) { 249 for (j in 1:length(classnames[[l]])) {
263 # i=class2names[j] 250 # i=class2names[j]
264 tic <- TIC[[classnames[[l]][j]]] 251 bpc <- BPC[[classnames[[l]][j]]]
265 points(tic[,1]/60, -tic[,2], col = cols[classnames[[l]][j]], pch = pch[classnames[[l]][j]], type="l") 252 points(bpc[,1]/60, -bpc[,2], col = cols[classnames[[l]][j]], pch = pch[classnames[[l]][j]], type="l")
266 colvect<-append(colvect,cols[classnames[[l]][j]]) 253 colvect<-append(colvect,cols[classnames[[l]][j]])
267 } 254 }
268 legend("topright",paste(gsub("(^.+)\\..*$","\\1",basename(files[c(classnames[[k]],classnames[[l]])]))), col = colvect, lty = lty, pch = pch) 255 legend("topright",paste(gsub("(^.+)\\..*$","\\1",basename(files[c(classnames[[k]],classnames[[l]])]))), col = colvect, lty = lty, pch = pch)
269 }#end length ==2 256 }#end length ==2
270 257
271 if (length(class)==1){ 258 if (length(class)==1){
272 k=1 259 k=1
273 ylim = range(sapply(TIC, function(x) range(x[,2]))) 260
261 ylim = range(sapply(BPC, function(x) range(x[,2])))
262
274 colvect<-NULL 263 colvect<-NULL
275 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") 264 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")
276 265
277 for (j in 1:length(classnames[[k]])) { 266 for (j in 1:length(classnames[[k]])) {
278 tic <- TIC[[classnames[[k]][j]]] 267 bpc <- BPC[[classnames[[k]][j]]]
279 # points(tic[,1]/60, tic[,2], col = cols[i], pch = pch[i], type="l") 268 # points(bpc[,1]/60, bpc[,2], col = cols[i], pch = pch[i], type="l")
280 points(tic[,1]/60, tic[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l") 269 points(bpc[,1]/60, bpc[,2], col = cols[classnames[[k]][j]], pch = pch[classnames[[k]][j]], type="l")
281 colvect<-append(colvect,cols[classnames[[k]][j]]) 270 colvect<-append(colvect,cols[classnames[[k]][j]])
282 } 271 }
283 legend("topright",paste(gsub("(^.+)\\..*$","\\1",basename(files[c(classnames[[k]])]))), col = colvect, lty = lty, pch = pch) 272 legend("topright",paste(gsub("(^.+)\\..*$","\\1",basename(files[c(classnames[[k]])]))), col = colvect, lty = lty, pch = pch)
284 }#end length ==1 273 }#end length ==1
285 dev.off() 274 dev.off()
291 } 280 }
292 281
293 ## overlay TIC from all files in current folder or from xcmsSet, create pdf 282 ## overlay TIC from all files in current folder or from xcmsSet, create pdf
294 getTIC2s <- function(files, xset=NULL, pdfname="TICs.pdf", rt=c("raw","corrected")) { 283 getTIC2s <- function(files, xset=NULL, pdfname="TICs.pdf", rt=c("raw","corrected")) {
295 require(xcms) 284 require(xcms)
296
297 #Verification for cdf files
298 stop=FALSE
299 for(i in 1:length(files)){
300 extension <- unlist(strsplit(basename(files[i]),"\\."))[length(unlist(strsplit(basename(files[i]),"\\.")))]
301 if(extension == "CDF" || extension == "cdf"){
302 stop = TRUE
303 break
304 }
305 }
306 if(stop){
307 error_message <- "You have a CDF file and there is an issue to resolve on them for chromatograms.... !"
308 print(error_message)
309 stop(error_message)
310 }
311 285
312 #create sampleMetadata, get sampleMetadata and class 286 #create sampleMetadata, get sampleMetadata and class
313 if(!is.null(xset)){ 287 if(!is.null(xset)){
314 #When files come from XCMS3 before metaMS treatment 288 #When files come from XCMS3 before metaMS treatment
315 sampleMetadata<-xset@phenoData 289 sampleMetadata<-xset@phenoData
325 299
326 N <- dim(sampleMetadata)[1] 300 N <- dim(sampleMetadata)[1]
327 TIC <- vector("list",N) 301 TIC <- vector("list",N)
328 302
329 for (i in 1:N) { 303 for (i in 1:N) {
330 cat(files[i],"\n")
331 if (!is.null(xcmsSet) && rt == "corrected") 304 if (!is.null(xcmsSet) && rt == "corrected")
332 rtcor <- xcmsSet@rt$corrected[[i]] 305 rtcor <- xcmsSet@rt$corrected[[i]]
333 else 306 else
334 rtcor <- NULL 307 rtcor <- NULL
335 TIC[[i]] <- getTIC(files[i],rtcor=rtcor) 308 TIC[[i]] <- getTIC(files[i],rtcor=rtcor)
346 319
347 ##plot start 320 ##plot start
348 if (length(class)>2){ 321 if (length(class)>2){
349 for (k in 1:(length(class)-1)){ 322 for (k in 1:(length(class)-1)){
350 for (l in (k+1):length(class)){ 323 for (l in (k+1):length(class)){
351 print(paste(class[k],"vs",class[l],sep=" ")) 324 cat(paste(class[k],"vs",class[l],"\n",sep=" "))
352 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") 325 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")
353 colvect<-NULL 326 colvect<-NULL
354 for (j in 1:length(classnames[[k]])) { 327 for (j in 1:length(classnames[[k]])) {
355 tic <- TIC[[classnames[[k]][j]]] 328 tic <- TIC[[classnames[[k]][j]]]
356 # points(tic[,1]/60, tic[,2], col = cols[i], pch = pch[i], type="l") 329 # points(tic[,1]/60, tic[,2], col = cols[i], pch = pch[i], type="l")
409 ##addition for quality control of peak picking 382 ##addition for quality control of peak picking
410 #metaMS EIC and pspectra plotting option 383 #metaMS EIC and pspectra plotting option
411 #version 20190520 384 #version 20190520
412 #only for Galaxy 385 #only for Galaxy
413 plotUnknowns<-function(resGC, unkn="", DB=NULL, fileFrom=NULL){ 386 plotUnknowns<-function(resGC, unkn="", DB=NULL, fileFrom=NULL){
414
415 #Verification for cdf files
416 stop=FALSE
417 for(i in 1:length(names(resGC$annotation))){
418 extension <- unlist(strsplit(basename(names(resGC$annotation)[i]),"\\."))[length(unlist(strsplit(basename(names(resGC$annotation)[i]),"\\.")))]
419 if(extension == "CDF" || extension == "cdf"){
420 stop = TRUE
421 break
422 }
423 }
424 if(stop){
425 error_message <- "You have a CDF file and there is an issue to resolve on them for chromatograms.... !"
426 print(error_message)
427 stop(error_message)
428 }
429 387
430 ##Annotation table each value is a pcgrp associated to the unknown 388 ##Annotation table each value is a pcgrp associated to the unknown
431 ##NOTE pcgrp index are different between xcmsSet and resGC due to filtering steps in metaMS 389 ##NOTE pcgrp index are different between xcmsSet and resGC due to filtering steps in metaMS
432 ##R. Wehrens give me some clues on that and we found a correction 390 ##R. Wehrens give me some clues on that and we found a correction
433 391