Mercurial > repos > workflow4metabolomics > camera_combinexsannos

comparison lib.r @ 0:139ff66b0b5d draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit f69695e76674862ed9c77c1c127f459b4df42464
author workflow4metabolomics
date Fri, 26 Jul 2019 16:49:18 -0400
parents
children ea15115a5b3f
comparison
equal deleted inserted replaced
-1:000000000000 0:139ff66b0b5d
1 # lib.r
2
3 #@author G. Le Corguille
4 # solve an issue with batch if arguments are logical TRUE/FALSE
5 parseCommandArgs <- function(...) {
6 args <- batch::parseCommandArgs(...)
7 for (key in names(args)) {
8 if (args[key] %in% c("TRUE","FALSE"))
9 args[key] = as.logical(args[key])
10 }
11 return(args)
12 }
13
14 #@author G. Le Corguille
15 # This function will
16 # - load the packages
17 # - display the sessionInfo
18 loadAndDisplayPackages <- function(pkgs) {
19 for(pkg in pkgs) suppressPackageStartupMessages( stopifnot( library(pkg, quietly=TRUE, logical.return=TRUE, character.only=TRUE)))
20
21 sessioninfo = sessionInfo()
22 cat(sessioninfo$R.version$version.string,"\n")
23 cat("Main packages:\n")
24 for (pkg in names(sessioninfo$otherPkgs)) { cat(paste(pkg,packageVersion(pkg)),"\t") }; cat("\n")
25 cat("Other loaded packages:\n")
26 for (pkg in names(sessioninfo$loadedOnly)) { cat(paste(pkg,packageVersion(pkg)),"\t") }; cat("\n")
27 }
28
29 # This function retrieve a xset like object
30 #@author Gildas Le Corguille lecorguille@sb-roscoff.fr
31 getxcmsSetObject <- function(xobject) {
32 # XCMS 1.x
33 if (class(xobject) == "xcmsSet")
34 return (xobject)
35 # XCMS 3.x
36 if (class(xobject) == "XCMSnExp") {
37 # Get the legacy xcmsSet object
38 suppressWarnings(xset <- as(xobject, 'xcmsSet'))
39 if (is.null(xset@phenoData$sample_group))
40 sampclass(xset) = "."
41 else
42 sampclass(xset) <- xset@phenoData$sample_group
43 if (!is.null(xset@phenoData$sample_name))
44 rownames(xset@phenoData) = xset@phenoData$sample_name
45 return (xset)
46 }
47 }
48
49 #@author G. Le Corguille
50 #The function create a pdf from the different png generated by diffreport
51 diffreport_png2pdf <- function(filebase) {
52 dir.create("pdf")
53
54 pdfEicOutput = paste0("pdf/",filebase,"-eic_pdf.pdf")
55 pdfBoxOutput = paste0("pdf/",filebase,"-box_pdf.pdf")
56
57 system(paste0("gm convert ",filebase,"_eic/*.png ",pdfEicOutput))
58 system(paste0("gm convert ",filebase,"_box/*.png ",pdfBoxOutput))
59
60 }
61
62 #@author G. Le Corguille
63 #The function create a zip archive from the different png generated by diffreport
64 diffreport_png2zip <- function() {
65 zip("eic.zip", dir(pattern="_eic"), zip=Sys.which("zip"))
66 zip("box.zip", dir(pattern="_box"), zip=Sys.which("zip"))
67 }
68
69 #The function create a zip archive from the different tabular generated by diffreport
70 diffreport_tabular2zip <- function() {
71 zip("tabular.zip", dir(pattern="tabular/*"), zip=Sys.which("zip"))
72 }
73
74 #@author G. Le Corguille
75 #This function convert if it is required the Retention Time in minutes
76 RTSecondToMinute <- function(variableMetadata, convertRTMinute) {
77 if (convertRTMinute){
78 #converting the retention times (seconds) into minutes
79 print("converting the retention times into minutes in the variableMetadata")
80 variableMetadata[,"rt"]=variableMetadata[,"rt"]/60
81 variableMetadata[,"rtmin"]=variableMetadata[,"rtmin"]/60
82 variableMetadata[,"rtmax"]=variableMetadata[,"rtmax"]/60
83 }
84 return (variableMetadata)
85 }
86
87 #@author G. Le Corguille
88 #This function format ions identifiers
89 formatIonIdentifiers <- function(variableMetadata, numDigitsRT=0, numDigitsMZ=0) {
90 splitDeco = strsplit(as.character(variableMetadata$name),"_")
91 idsDeco = sapply(splitDeco, function(x) { deco=unlist(x)[2]; if (is.na(deco)) return ("") else return(paste0("_",deco)) })
92 namecustom = make.unique(paste0("M",round(variableMetadata[,"mz"],numDigitsMZ),"T",round(variableMetadata[,"rt"],numDigitsRT),idsDeco))
93 variableMetadata=cbind(name=variableMetadata$name, namecustom=namecustom, variableMetadata[,!(colnames(variableMetadata) %in% c("name"))])
94 return(variableMetadata)
95 }
96
97 #The function annotateDiffreport without the corr function which bugs
98 annotatediff <- function(xset=xset, args=args, variableMetadataOutput="variableMetadata.tsv") {
99 # Resolve the bug with x11, with the function png
100 options(bitmapType='cairo')
101
102 #Check if the fillpeaks step has been done previously, if it hasn't, there is an error message and the execution is stopped.
103 res=try(is.null(xset@filled))
104
105 # ------ annot -------
106 args$calcCiS=as.logical(args$calcCiS)
107 args$calcIso=as.logical(args$calcIso)
108 args$calcCaS=as.logical(args$calcCaS)
109
110 # common parameters
111 args4annotate = list(object=xset,
112 nSlaves=args$nSlaves,sigma=args$sigma,perfwhm=args$perfwhm,
113 maxcharge=args$maxcharge,maxiso=args$maxiso,minfrac=args$minfrac,
114 ppm=args$ppm,mzabs=args$mzabs,quick=args$quick,
115 polarity=args$polarity,max_peaks=args$max_peaks,intval=args$intval)
116
117 # quick == FALSE
118 if(args$quick==FALSE) {
119 args4annotate = append(args4annotate,
120 list(graphMethod=args$graphMethod,cor_eic_th=args$cor_eic_th,pval=args$pval,
121 calcCiS=args$calcCiS,calcIso=args$calcIso,calcCaS=args$calcCaS))
122 # no ruleset
123 if (!is.null(args$multiplier)) {
124 args4annotate = append(args4annotate,
125 list(multiplier=args$multiplier))
126 }
127 # ruleset
128 else {
129 rulset=read.table(args$rules, h=T, sep=";")
130 if (ncol(rulset) < 4) rulset=read.table(args$rules, h=T, sep="\t")
131 if (ncol(rulset) < 4) rulset=read.table(args$rules, h=T, sep=",")
132 if (ncol(rulset) < 4) {
133 error_message="Your ruleset file seems not well formatted. The column separators accepted are ; , and tabulation"
134 print(error_message)
135 stop(error_message)
136 }
137
138 args4annotate = append(args4annotate,
139 list(rules=rulset))
140 }
141 }
142
143
144 # launch annotate
145 xa = do.call("annotate", args4annotate)
146 peakList=getPeaklist(xa,intval=args$intval)
147 peakList=cbind(groupnames(xa@xcmsSet),peakList); colnames(peakList)[1] = c("name");
148
149 # --- Multi condition : diffreport ---
150 diffrepOri=NULL
151 if (!is.null(args$runDiffreport) & nlevels(sampclass(xset))>=2) {
152 #Check if the fillpeaks step has been done previously, if it hasn't, there is an error message and the execution is stopped.
153 res=try(is.null(xset@filled))
154 classes=levels(sampclass(xset))
155 x=1:(length(classes)-1)
156 for (i in seq(along=x) ) {
157 y=1:(length(classes))
158 for (n in seq(along=y)){
159 if(i+n <= length(classes)){
160 filebase=paste(classes[i],class2=classes[i+n],sep="-vs-")
161
162 diffrep=diffreport(
163 object=xset,class1=classes[i],class2=classes[i+n],
164 filebase=filebase,eicmax=args$eicmax,eicwidth=args$eicwidth,
165 sortpval=TRUE,value=args$value,h=args$h,w=args$w,mzdec=args$mzdec,missing=0)
166
167 diffrepOri = diffrep
168
169 # renamming of the column rtmed to rt to fit with camera peaklist function output
170 colnames(diffrep)[colnames(diffrep)=="rtmed"] <- "rt"
171 colnames(diffrep)[colnames(diffrep)=="mzmed"] <- "mz"
172
173 # combines results and reorder columns
174 diffrep = merge(peakList, diffrep[,c("name","fold","tstat","pvalue")], by.x="name", by.y="name", sort=F)
175 diffrep = cbind(diffrep[,!(colnames(diffrep) %in% c(sampnames(xa@xcmsSet)))],diffrep[,(colnames(diffrep) %in% c(sampnames(xa@xcmsSet)))])
176
177 diffrep = RTSecondToMinute(diffrep, args$convertRTMinute)
178 diffrep = formatIonIdentifiers(diffrep, numDigitsRT=args$numDigitsRT, numDigitsMZ=args$numDigitsMZ)
179
180 if(args$sortpval){
181 diffrep=diffrep[order(diffrep$pvalue), ]
182 }
183
184 dir.create("tabular", showWarnings = FALSE)
185 write.table(diffrep, sep="\t", quote=FALSE, row.names=FALSE, file=paste("tabular/",filebase,"_tsv.tabular",sep=""))
186
187 if (args$eicmax != 0) {
188 if (args$png2 == "pdf")
189 diffreport_png2pdf(filebase)
190 }
191 }
192 }
193 }
194 if (args$png2 == "zip")
195 diffreport_png2zip()
196 if (args$tabular2 == "zip")
197 diffreport_tabular2zip()
198 }
199
200 # --- variableMetadata ---
201 variableMetadata=peakList[,!(make.names(colnames(peakList)) %in% c(make.names(sampnames(xa@xcmsSet))))]
202 variableMetadata = RTSecondToMinute(variableMetadata, args$convertRTMinute)
203 variableMetadata = formatIonIdentifiers(variableMetadata, numDigitsRT=args$numDigitsRT, numDigitsMZ=args$numDigitsMZ)
204 # if we have 2 conditions, we keep stat of diffrep
205 if (!is.null(args$runDiffreport) & nlevels(sampclass(xset))==2) {
206 variableMetadata = merge(variableMetadata, diffrep[,c("name","fold","tstat","pvalue")],by.x="name", by.y="name", sort=F)
207 if(exists("args[[\"sortpval\"]]")){
208 variableMetadata=variableMetadata[order(variableMetadata$pvalue), ]
209 }
210 }
211
212 variableMetadataOri=variableMetadata
213 write.table(variableMetadata, sep="\t", quote=FALSE, row.names=FALSE, file=variableMetadataOutput)
214
215 return(list("xa"=xa,"diffrep"=diffrepOri,"variableMetadata"=variableMetadataOri));
216
217 }
218
219
220 combinexsAnnos_function <- function(xaP, xaN, diffrepP=NULL,diffrepN=NULL,
221 pos=TRUE,tol=2,ruleset=NULL,keep_meta=TRUE, convertRTMinute=F, numDigitsMZ=0,
222 numDigitsRT=0, variableMetadataOutput="variableMetadata.tsv"){
223
224 #Load the two Rdata to extract the xset objects from positive and negative mode
225 cat("\tObject xset from positive mode\n")
226 print(xaP)
227 cat("\n")
228
229 cat("\tObject xset from negative mode\n")
230 print(xaN)
231 cat("\n")
232
233 cat("\n")
234 cat("\tCombining...\n")
235 #Convert the string to numeric for creating matrix
236 row=as.numeric(strsplit(ruleset,",")[[1]][1])
237 column=as.numeric(strsplit(ruleset,",")[[1]][2])
238 ruleset=cbind(row,column)
239 #Test if the file comes from an older version tool
240 if ((!is.null(xaP)) & (!is.null(xaN))) {
241 #Launch the combinexsannos function from CAMERA
242 cAnnot=combinexsAnnos(xaP, xaN,pos=pos,tol=tol,ruleset=ruleset)
243 } else {
244 stop("You must relauch the CAMERA.annotate step with the lastest version.")
245 }
246
247 if(pos){
248 xa=xaP
249 mode="neg. Mode"
250 } else {
251 xa=xaN
252 mode="pos. Mode"
253 }
254
255 peakList=getPeaklist(xa)
256 peakList=cbind(groupnames(xa@xcmsSet),peakList); colnames(peakList)[1] = c("name");
257 variableMetadata=cbind(peakList, cAnnot[, c("isotopes", "adduct", "pcgroup",mode)]);
258 variableMetadata=variableMetadata[,!(colnames(variableMetadata) %in% c(sampnames(xa@xcmsSet)))]
259
260 #Test if there are more than two classes (conditions)
261 if ( nlevels(sampclass(xaP@xcmsSet))==2 & (!is.null(diffrepN)) & (!is.null(diffrepP))) {
262 diffrepP = diffrepP[,c("name","fold","tstat","pvalue")]; colnames(diffrepP) = paste("P.",colnames(diffrepP),sep="")
263 diffrepN = diffrepN[,c("name","fold","tstat","pvalue")]; colnames(diffrepN) = paste("N.",colnames(diffrepN),sep="")
264
265 variableMetadata = merge(variableMetadata, diffrepP, by.x="name", by.y="P.name")
266 variableMetadata = merge(variableMetadata, diffrepN, by.x="name", by.y="N.name")
267 }
268
269 rownames(variableMetadata) = NULL
270 #TODO: checker
271 #colnames(variableMetadata)[1:2] = c("name","mz/rt");
272
273 variableMetadata = RTSecondToMinute(variableMetadata, convertRTMinute)
274 variableMetadata = formatIonIdentifiers(variableMetadata, numDigitsRT=numDigitsRT, numDigitsMZ=numDigitsMZ)
275
276 #If the user want to keep only the metabolites which match a difference
277 if(keep_meta){
278 variableMetadata=variableMetadata[variableMetadata[,c(mode)]!="",]
279 }
280
281 #Write the output into a tsv file
282 write.table(variableMetadata, sep="\t", quote=FALSE, row.names=FALSE, file=variableMetadataOutput)
283 return(variableMetadata);
284
285 }
286
287 # This function get the raw file path from the arguments
288 getRawfilePathFromArguments <- function(singlefile, zipfile, args) {
289 if (!is.null(args$zipfile)) zipfile = args$zipfile
290 if (!is.null(args$zipfilePositive)) zipfile = args$zipfilePositive
291 if (!is.null(args$zipfileNegative)) zipfile = args$zipfileNegative
292
293 if (!is.null(args$singlefile_galaxyPath)) {
294 singlefile_galaxyPaths = args$singlefile_galaxyPath;
295 singlefile_sampleNames = args$singlefile_sampleName
296 }
297 if (!is.null(args$singlefile_galaxyPathPositive)) {
298 singlefile_galaxyPaths = args$singlefile_galaxyPathPositive;
299 singlefile_sampleNames = args$singlefile_sampleNamePositive
300 }
301 if (!is.null(args$singlefile_galaxyPathNegative)) {
302 singlefile_galaxyPaths = args$singlefile_galaxyPathNegative;
303 singlefile_sampleNames = args$singlefile_sampleNameNegative
304 }
305 if (exists("singlefile_galaxyPaths")){
306 singlefile_galaxyPaths = unlist(strsplit(singlefile_galaxyPaths,","))
307 singlefile_sampleNames = unlist(strsplit(singlefile_sampleNames,","))
308
309 singlefile=NULL
310 for (singlefile_galaxyPath_i in seq(1:length(singlefile_galaxyPaths))) {
311 singlefile_galaxyPath=singlefile_galaxyPaths[singlefile_galaxyPath_i]
312 singlefile_sampleName=singlefile_sampleNames[singlefile_galaxyPath_i]
313 singlefile[[singlefile_sampleName]] = singlefile_galaxyPath
314 }
315 }
316 for (argument in c("zipfile", "zipfilePositive", "zipfileNegative",
317 "singlefile_galaxyPath", "singlefile_sampleName",
318 "singlefile_galaxyPathPositive", "singlefile_sampleNamePositive",
319 "singlefile_galaxyPathNegative","singlefile_sampleNameNegative")) {
320 args[[argument]]=NULL
321 }
322 return(list(zipfile=zipfile, singlefile=singlefile, args=args))
323 }
324
325
326 # This function retrieve the raw file in the working directory
327 # - if zipfile: unzip the file with its directory tree
328 # - if singlefiles: set symlink with the good filename
329 retrieveRawfileInTheWorkingDirectory <- function(singlefile, zipfile) {
330 if(!is.null(singlefile) && (length("singlefile")>0)) {
331 for (singlefile_sampleName in names(singlefile)) {
332 singlefile_galaxyPath = singlefile[[singlefile_sampleName]]
333 if(!file.exists(singlefile_galaxyPath)){
334 error_message=paste("Cannot access the sample:",singlefile_sampleName,"located:",singlefile_galaxyPath,". Please, contact your administrator ... if you have one!")
335 print(error_message); stop(error_message)
336 }
337
338 file.symlink(singlefile_galaxyPath,singlefile_sampleName)
339 }
340 directory = "."
341
342 }
343 if(!is.null(zipfile) && (zipfile!="")) {
344 if(!file.exists(zipfile)){
345 error_message=paste("Cannot access the Zip file:",zipfile,". Please, contact your administrator ... if you have one!")
346 print(error_message)
347 stop(error_message)
348 }
349
350 #list all file in the zip file
351 #zip_files=unzip(zipfile,list=T)[,"Name"]
352
353 #unzip
354 suppressWarnings(unzip(zipfile, unzip="unzip"))
355
356 #get the directory name
357 filesInZip=unzip(zipfile, list=T);
358 directories=unique(unlist(lapply(strsplit(filesInZip$Name,"/"), function(x) x[1])));
359 directories=directories[!(directories %in% c("__MACOSX")) & file.info(directories)$isdir]
360 directory = "."
361 if (length(directories) == 1) directory = directories
362
363 cat("files_root_directory\t",directory,"\n")
364
365 }
366 return (directory)
367 }
368
369 #@TODO: remove this function as soon as we can use xcms 3.x.x from Bioconductor 3.7
370 # https://github.com/sneumann/CAMERA/issues/33#issuecomment-405168524
371 # https://github.com/sneumann/xcms/commit/950a3fe794cdb6b0fda88696e31aab3d97a3b7dd
372 ############################################################
373 ## getEIC
374 getEIC <- function(object, mzrange, rtrange = 200,
375 groupidx, sampleidx = sampnames(object),
376 rt = c("corrected", "raw")) {
377
378 files <- filepaths(object)
379 grp <- groups(object)
380 samp <- sampnames(object)
381 prof <- profinfo(object)
382
383 rt <- match.arg(rt)
384
385 if (is.numeric(sampleidx))
386 sampleidx <- sampnames(object)[sampleidx]
387 sampidx <- match(sampleidx, sampnames(object))
388
389 if (!missing(groupidx)) {
390 if (is.numeric(groupidx))
391 groupidx <- groupnames(object)[unique(as.integer(groupidx))]
392 grpidx <- match(groupidx, groupnames(object, template = groupidx))
393 }
394
395 if (missing(mzrange)) {
396 if (missing(groupidx))
397 stop("No m/z range or groups specified")
398 if (any(is.na(groupval(object, value = "mz"))))
399 warning(
400 "`NA` values in xcmsSet. Use fillPeaks() on the object to fill",
401 "-in missing peak values. Note however that this will also ",
402 "insert intensities of 0 for peaks that can not be filled in.")
403 mzmin <- apply(groupval(object, value = "mzmin"), 1, min, na.rm = TRUE)
404 mzmax <- apply(groupval(object, value = "mzmax"), 1, max, na.rm = TRUE)
405 mzrange <- matrix(c(mzmin[grpidx], mzmax[grpidx]), ncol = 2)
406 ## if (any(is.na(groupval(object, value = "mz"))))
407 ## stop('Please use fillPeaks() to fill up NA values !')
408 ## mzmin <- -rowMax(-groupval(object, value = "mzmin"))
409 ## mzmax <- rowMax(groupval(object, value = "mzmax"))
410 ## mzrange <- matrix(c(mzmin[grpidx], mzmax[grpidx]), ncol = 2)
411 } else if (all(c("mzmin","mzmax") %in% colnames(mzrange)))
412 mzrange <- mzrange[,c("mzmin", "mzmax"),drop=FALSE]
413 else if (is.null(dim(mzrange)))
414 stop("mzrange must be a matrix")
415 colnames(mzrange) <- c("mzmin", "mzmax")
416
417 if (length(rtrange) == 1) {
418 if (missing(groupidx))
419 rtrange <- matrix(rep(range(object@rt[[rt]][sampidx]), nrow(mzrange)),
420 ncol = 2, byrow = TRUE)
421 else {
422 rtrange <- retexp(grp[grpidx,c("rtmin","rtmax"),drop=FALSE], rtrange)
423 }
424 } else if (is.null(dim(rtrange)))
425 stop("rtrange must be a matrix or single number")
426 colnames(rtrange) <- c("rtmin", "rtmax")
427
428 ## Ensure that we've got corrected retention time if requested.
429 if (is.null(object@rt[[rt]]))
430 stop(rt, " retention times not present in 'object'!")
431
432 ## Ensure that the defined retention time range is within the rtrange of the
433 ## object: we're using the max minimal rt of all files and the min maximal rt
434 rtrs <- lapply(object@rt[[rt]], range)
435 rtr <- c(max(unlist(lapply(rtrs, "[", 1))),
436 min(unlist(lapply(rtrs, "[", 2))))
437 ## Check if we've got a range which is completely off:
438 if (any(rtrange[, "rtmin"] >= rtr[2] | rtrange[, "rtmax"] <= rtr[1])) {
439 outs <- which(rtrange[, "rtmin"] >= rtr[2] |
440 rtrange[, "rtmax"] <= rtr[1])
441 stop(length(outs), " of the specified 'rtrange' are completely outside ",
442 "of the retention time range of 'object' which is (", rtr[1], ", ",
443 rtr[2], "). The first was: (", rtrange[outs[1], "rtmin"], ", ",
444 rtrange[outs[1], "rtmax"], "!")
445 }
446 lower_rt_outside <- rtrange[, "rtmin"] < rtr[1]
447 upper_rt_outside <- rtrange[, "rtmax"] > rtr[2]
448 if (any(lower_rt_outside) | any(upper_rt_outside)) {
449 ## Silently fix these ranges.
450 rtrange[lower_rt_outside, "rtmin"] <- rtr[1]
451 rtrange[upper_rt_outside, "rtmax"] <- rtr[2]
452 }
453
454 if (missing(groupidx))
455 gnames <- character(0)
456 else
457 gnames <- groupidx
458
459 eic <- vector("list", length(sampleidx))
460 names(eic) <- sampleidx
461
462 for (i in seq(along = sampidx)) {
463
464 ## cat(sampleidx[i], "")
465 flush.console()
466 ## getXcmsRaw takes care of rt correction, susetting to scanrage and other
467 ## stuff.
468 lcraw <- getXcmsRaw(object, sampleidx = sampidx[i], rt=rt)
469 currenteic <- xcms::getEIC(lcraw, mzrange, rtrange, step = prof$step)
470 eic[[i]] <- currenteic@eic[[1]]
471 rm(lcraw)
472 gc()
473 }
474 ## cat("\n")
475
476 invisible(new("xcmsEIC", eic = eic, mzrange = mzrange, rtrange = rtrange,
477 rt = rt, groupnames = gnames))
478 }
479
480 #@TODO: remove this function as soon as we can use xcms 3.x.x from Bioconductor 3.7
481 # https://github.com/sneumann/CAMERA/issues/33#issuecomment-405168524
482 # https://github.com/sneumann/xcms/commit/950a3fe794cdb6b0fda88696e31aab3d97a3b7dd
483 ############################################################
484 ## diffreport
485 diffreport = function(object,
486 class1 = levels(sampclass(object))[1],
487 class2 = levels(sampclass(object))[2],
488 filebase = character(),
489 eicmax = 0, eicwidth = 200,
490 sortpval = TRUE,
491 classeic = c(class1,class2),
492 value = c("into","maxo","intb"),
493 metlin = FALSE,
494 h = 480, w = 640, mzdec=2,
495 missing = numeric(), ...) {
496
497 if ( nrow(object@groups)<1 || length(object@groupidx) <1) {
498 stop("No group information. Use group().")
499 }
500
501 if (!is.numeric(w) || !is.numeric(h))
502 stop("'h' and 'w' have to be numeric")
503 ## require(multtest) || stop("Couldn't load multtest")
504
505 value <- match.arg(value)
506 groupmat <- groups(object)
507 if (length(groupmat) == 0)
508 stop("No group information found")
509 samples <- sampnames(object)
510 n <- length(samples)
511 classlabel <- sampclass(object)
512 classlabel <- levels(classlabel)[as.vector(unclass(classlabel))]
513
514 values <- groupval(object, "medret", value=value)
515 indecies <- groupval(object, "medret", value = "index")
516
517 if (!all(c(class1,class2) %in% classlabel))
518 stop("Incorrect Class Labels")
519
520 ## c1 and c2 are column indices of class1 and class2 resp.
521 c1 <- which(classlabel %in% class1)
522 c2 <- which(classlabel %in% class2)
523 ceic <- which(classlabel %in% classeic)
524 if (length(intersect(c1, c2)) > 0)
525 stop("Intersecting Classes")
526
527 ## Optionally replace NA values with the value provided with missing
528 if (length(missing)) {
529 if (is.numeric(missing)) {
530 ## handles NA, Inf and -Inf
531 values[, c(c1, c2)][!is.finite(values[, c(c1, c2)])] <- missing[1]
532 } else
533 stop("'missing' should be numeric")
534 }
535 ## Check against missing Values
536 if (any(is.na(values[, c(c1, c2)])))
537 warning("`NA` values in xcmsSet. Use fillPeaks() on the object to fill",
538 "-in missing peak values. Note however that this will also ",
539 "insert intensities of 0 for peaks that can not be filled in.")
540
541 mean1 <- rowMeans(values[,c1,drop=FALSE], na.rm = TRUE)
542 mean2 <- rowMeans(values[,c2,drop=FALSE], na.rm = TRUE)
543
544 ## Calculate fold change.
545 ## For foldchange <1 set fold to 1/fold
546 ## See tstat to check which was higher
547 fold <- mean2 / mean1
548 fold[!is.na(fold) & fold < 1] <- 1/fold[!is.na(fold) & fold < 1]
549
550 testval <- values[,c(c1,c2)]
551 ## Replace eventual infinite values with NA (CAMERA issue #33)
552 testval[is.infinite(testval)] <- NA
553 testclab <- c(rep(0,length(c1)),rep(1,length(c2)))
554
555 if (min(length(c1), length(c2)) >= 2) {
556 tstat <- mt.teststat(testval, testclab, ...)
557 pvalue <- xcms:::pval(testval, testclab, tstat)
558 } else {
559 message("Too few samples per class, skipping t-test.")
560 tstat <- pvalue <- rep(NA,nrow(testval))
561 }
562 stat <- data.frame(fold = fold, tstat = tstat, pvalue = pvalue)
563 if (length(levels(sampclass(object))) >2) {
564 pvalAnova<-c()
565 for(i in 1:nrow(values)){
566 var<-as.numeric(values[i,])
567 ano<-summary(aov(var ~ sampclass(object)) )
568 pvalAnova<-append(pvalAnova, unlist(ano)["Pr(>F)1"])
569 }
570 stat<-cbind(stat, anova= pvalAnova)
571 }
572 if (metlin) {
573 neutralmass <- groupmat[,"mzmed"] + ifelse(metlin < 0, 1, -1)
574 metlin <- abs(metlin)
575 digits <- ceiling(-log10(metlin))+1
576 metlinurl <-
577 paste("http://metlin.scripps.edu/simple_search_result.php?mass_min=",
578 round(neutralmass - metlin, digits), "&mass_max=",
579 round(neutralmass + metlin, digits), sep="")
580 values <- cbind(metlin = metlinurl, values)
581 }
582 twosamp <- cbind(name = groupnames(object), stat, groupmat, values)
583 if (sortpval) {
584 tsidx <- order(twosamp[,"pvalue"])
585 twosamp <- twosamp[tsidx,]
586 rownames(twosamp) <- 1:nrow(twosamp)
587 values<-values[tsidx,]
588 } else
589 tsidx <- 1:nrow(values)
590
591 if (length(filebase))
592 write.table(twosamp, paste(filebase, ".tsv", sep = ""), quote = FALSE, sep = "\t", col.names = NA)
593
594 if (eicmax > 0) {
595 if (length(unique(peaks(object)[,"rt"])) > 1) {
596 ## This looks like "normal" LC data
597
598 eicmax <- min(eicmax, length(tsidx))
599 eics <- getEIC(object, rtrange = eicwidth*1.1, sampleidx = ceic,
600 groupidx = tsidx[seq(length = eicmax)])
601
602 if (length(filebase)) {
603 eicdir <- paste(filebase, "_eic", sep="")
604 boxdir <- paste(filebase, "_box", sep="")
605 dir.create(eicdir)
606 dir.create(boxdir)
607 if (capabilities("png")){
608 xcms:::xcmsBoxPlot(values[seq(length = eicmax),],
609 sampclass(object), dirpath=boxdir, pic="png", width=w, height=h)
610 png(file.path(eicdir, "%003d.png"), width = w, height = h)
611 } else {
612 xcms:::xcmsBoxPlot(values[seq(length = eicmax),],
613 sampclass(object), dirpath=boxdir, pic="pdf", width=w, height=h)
614 pdf(file.path(eicdir, "%003d.pdf"), width = w/72,
615 height = h/72, onefile = FALSE)
616 }
617 }
618 plot(eics, object, rtrange = eicwidth, mzdec=mzdec)
619
620 if (length(filebase))
621 dev.off()
622 } else {
623 ## This looks like a direct-infusion single spectrum
624 if (length(filebase)) {
625 eicdir <- paste(filebase, "_eic", sep="")
626 boxdir <- paste(filebase, "_box", sep="")
627 dir.create(eicdir)
628 dir.create(boxdir)
629 if (capabilities("png")){
630 xcmsBoxPlot(values[seq(length = eicmax),],
631 sampclass(object), dirpath=boxdir, pic="png",
632 width=w, height=h)
633 png(file.path(eicdir, "%003d.png"), width = w, height = h,
634 units = "px")
635 } else {
636 xcmsBoxPlot(values[seq(length = eicmax),],
637 sampclass(object), dirpath=boxdir, pic="pdf",
638 width=w, height=h)
639 pdf(file.path(eicdir, "%003d.pdf"), width = w/72,
640 height = h/72, onefile = FALSE)
641 }
642 }
643
644 plotSpecWindow(object, gidxs = tsidx[seq(length = eicmax)], borderwidth=1)
645
646 if (length(filebase))
647 dev.off()
648 }
649 }
650
651 invisible(twosamp)
652 }