diff lib.r @ 1:ea15115a5b3f draft

"planemo upload commit 4fcbbcbc6d6b0a59e801870d31fe886a920ef429"
author workflow4metabolomics
date Thu, 13 Feb 2020 17:23:27 -0500
parents 139ff66b0b5d
children c4c13745e797
line wrap: on
line diff
--- a/lib.r	Fri Jul 26 16:49:18 2019 -0400
+++ b/lib.r	Thu Feb 13 17:23:27 2020 -0500
@@ -366,287 +366,3 @@
     return (directory)
 }
 
-#@TODO: remove this function as soon as we can use xcms 3.x.x from Bioconductor 3.7
-# https://github.com/sneumann/CAMERA/issues/33#issuecomment-405168524
-# https://github.com/sneumann/xcms/commit/950a3fe794cdb6b0fda88696e31aab3d97a3b7dd
-############################################################
-## getEIC
-getEIC <- function(object, mzrange, rtrange = 200,
-                   groupidx, sampleidx = sampnames(object),
-                   rt = c("corrected", "raw")) {
-
-    files <- filepaths(object)
-    grp <- groups(object)
-    samp <- sampnames(object)
-    prof <- profinfo(object)
-
-    rt <- match.arg(rt)
-
-    if (is.numeric(sampleidx))
-    sampleidx <- sampnames(object)[sampleidx]
-    sampidx <- match(sampleidx, sampnames(object))
-
-    if (!missing(groupidx)) {
-        if (is.numeric(groupidx))
-        groupidx <- groupnames(object)[unique(as.integer(groupidx))]
-        grpidx <- match(groupidx, groupnames(object, template = groupidx))
-    }
-
-    if (missing(mzrange)) {
-        if (missing(groupidx))
-        stop("No m/z range or groups specified")
-        if (any(is.na(groupval(object, value = "mz"))))
-        warning(
-        "`NA` values in xcmsSet. Use fillPeaks() on the object to fill",
-        "-in missing peak values. Note however that this will also ",
-        "insert intensities of 0 for peaks that can not be filled in.")
-        mzmin <- apply(groupval(object, value = "mzmin"), 1, min, na.rm = TRUE)
-        mzmax <- apply(groupval(object, value = "mzmax"), 1, max, na.rm = TRUE)
-        mzrange <- matrix(c(mzmin[grpidx], mzmax[grpidx]), ncol = 2)
-        ## if (any(is.na(groupval(object, value = "mz"))))
-        ##     stop('Please use fillPeaks() to fill up NA values !')
-        ## mzmin <- -rowMax(-groupval(object, value = "mzmin"))
-        ## mzmax <- rowMax(groupval(object, value = "mzmax"))
-        ## mzrange <- matrix(c(mzmin[grpidx], mzmax[grpidx]), ncol = 2)
-    } else if (all(c("mzmin","mzmax") %in% colnames(mzrange)))
-    mzrange <- mzrange[,c("mzmin", "mzmax"),drop=FALSE]
-    else if (is.null(dim(mzrange)))
-    stop("mzrange must be a matrix")
-    colnames(mzrange) <- c("mzmin", "mzmax")
-
-    if (length(rtrange) == 1) {
-        if (missing(groupidx))
-        rtrange <- matrix(rep(range(object@rt[[rt]][sampidx]), nrow(mzrange)),
-        ncol = 2, byrow = TRUE)
-        else {
-            rtrange <- retexp(grp[grpidx,c("rtmin","rtmax"),drop=FALSE], rtrange)
-        }
-    } else if (is.null(dim(rtrange)))
-    stop("rtrange must be a matrix or single number")
-    colnames(rtrange) <- c("rtmin", "rtmax")
-
-    ## Ensure that we've got corrected retention time if requested.
-    if (is.null(object@rt[[rt]]))
-    stop(rt, " retention times not present in 'object'!")
-
-    ## Ensure that the defined retention time range is within the rtrange of the
-    ## object: we're using the max minimal rt of all files and the min maximal rt
-    rtrs <- lapply(object@rt[[rt]], range)
-    rtr <- c(max(unlist(lapply(rtrs, "[", 1))),
-    min(unlist(lapply(rtrs, "[", 2))))
-    ## Check if we've got a range which is completely off:
-    if (any(rtrange[, "rtmin"] >= rtr[2] | rtrange[, "rtmax"] <= rtr[1])) {
-        outs <- which(rtrange[, "rtmin"] >= rtr[2] |
-        rtrange[, "rtmax"] <= rtr[1])
-        stop(length(outs), " of the specified 'rtrange' are completely outside ",
-        "of the retention time range of 'object' which is (", rtr[1], ", ",
-        rtr[2], "). The first was: (", rtrange[outs[1], "rtmin"], ", ",
-        rtrange[outs[1], "rtmax"], "!")
-    }
-    lower_rt_outside <- rtrange[, "rtmin"] < rtr[1]
-    upper_rt_outside <- rtrange[, "rtmax"] > rtr[2]
-    if (any(lower_rt_outside) | any(upper_rt_outside)) {
-        ## Silently fix these ranges.
-        rtrange[lower_rt_outside, "rtmin"] <- rtr[1]
-        rtrange[upper_rt_outside, "rtmax"] <- rtr[2]
-    }
-
-    if (missing(groupidx))
-    gnames <- character(0)
-    else
-    gnames <- groupidx
-
-    eic <- vector("list", length(sampleidx))
-    names(eic) <- sampleidx
-
-    for (i in seq(along = sampidx)) {
-
-        ## cat(sampleidx[i], "")
-        flush.console()
-        ## getXcmsRaw takes care of rt correction, susetting to scanrage and other
-        ## stuff.
-        lcraw <- getXcmsRaw(object, sampleidx = sampidx[i], rt=rt)
-        currenteic <- xcms::getEIC(lcraw, mzrange, rtrange, step = prof$step)
-        eic[[i]] <- currenteic@eic[[1]]
-        rm(lcraw)
-        gc()
-    }
-    ## cat("\n")
-
-    invisible(new("xcmsEIC", eic = eic, mzrange = mzrange, rtrange = rtrange,
-    rt = rt, groupnames = gnames))
-}
-
-#@TODO: remove this function as soon as we can use xcms 3.x.x from Bioconductor 3.7
-# https://github.com/sneumann/CAMERA/issues/33#issuecomment-405168524
-# https://github.com/sneumann/xcms/commit/950a3fe794cdb6b0fda88696e31aab3d97a3b7dd
-############################################################
-## diffreport
-diffreport = function(object,
-                      class1 = levels(sampclass(object))[1],
-                      class2 = levels(sampclass(object))[2],
-                      filebase = character(),
-                      eicmax = 0, eicwidth = 200,
-                      sortpval = TRUE,
-                      classeic = c(class1,class2),
-                      value = c("into","maxo","intb"),
-                      metlin = FALSE,
-                      h = 480, w = 640, mzdec=2,
-                      missing = numeric(), ...) {
-
-    if ( nrow(object@groups)<1 || length(object@groupidx) <1) {
-        stop("No group information. Use group().")
-    }
-
-    if (!is.numeric(w) || !is.numeric(h))
-        stop("'h' and 'w' have to be numeric")
-    ## require(multtest) || stop("Couldn't load multtest")
-
-    value <- match.arg(value)
-    groupmat <- groups(object)
-    if (length(groupmat) == 0)
-    stop("No group information found")
-    samples <- sampnames(object)
-    n <- length(samples)
-    classlabel <- sampclass(object)
-    classlabel <- levels(classlabel)[as.vector(unclass(classlabel))]
-
-    values <- groupval(object, "medret", value=value)
-    indecies <- groupval(object, "medret", value = "index")
-
-    if (!all(c(class1,class2) %in% classlabel))
-    stop("Incorrect Class Labels")
-
-    ## c1 and c2 are column indices of class1 and class2 resp.
-    c1 <- which(classlabel %in% class1)
-    c2 <- which(classlabel %in% class2)
-    ceic <- which(classlabel %in% classeic)
-    if (length(intersect(c1, c2)) > 0)
-    stop("Intersecting Classes")
-
-    ## Optionally replace NA values with the value provided with missing
-    if (length(missing)) {
-        if (is.numeric(missing)) {
-            ## handles NA, Inf and -Inf
-            values[, c(c1, c2)][!is.finite(values[, c(c1, c2)])] <- missing[1]
-        } else
-        stop("'missing' should be numeric")
-    }
-    ## Check against missing Values
-    if (any(is.na(values[, c(c1, c2)])))
-    warning("`NA` values in xcmsSet. Use fillPeaks() on the object to fill",
-    "-in missing peak values. Note however that this will also ",
-    "insert intensities of 0 for peaks that can not be filled in.")
-
-    mean1 <- rowMeans(values[,c1,drop=FALSE], na.rm = TRUE)
-    mean2 <- rowMeans(values[,c2,drop=FALSE], na.rm = TRUE)
-
-    ## Calculate fold change.
-    ## For foldchange <1 set fold to 1/fold
-    ## See tstat to check which was higher
-    fold <- mean2 / mean1
-    fold[!is.na(fold) & fold < 1] <- 1/fold[!is.na(fold) & fold < 1]
-
-    testval <- values[,c(c1,c2)]
-    ## Replace eventual infinite values with NA (CAMERA issue #33)
-    testval[is.infinite(testval)] <- NA
-    testclab <- c(rep(0,length(c1)),rep(1,length(c2)))
-
-    if (min(length(c1), length(c2)) >= 2) {
-        tstat <- mt.teststat(testval, testclab, ...)
-        pvalue <- xcms:::pval(testval, testclab, tstat)
-    } else {
-        message("Too few samples per class, skipping t-test.")
-        tstat <- pvalue <- rep(NA,nrow(testval))
-    }
-    stat <- data.frame(fold = fold, tstat = tstat, pvalue = pvalue)
-    if (length(levels(sampclass(object))) >2) {
-        pvalAnova<-c()
-        for(i in 1:nrow(values)){
-            var<-as.numeric(values[i,])
-            ano<-summary(aov(var ~ sampclass(object)) )
-            pvalAnova<-append(pvalAnova, unlist(ano)["Pr(>F)1"])
-        }
-        stat<-cbind(stat, anova= pvalAnova)
-    }
-    if (metlin) {
-        neutralmass <- groupmat[,"mzmed"] + ifelse(metlin < 0, 1, -1)
-        metlin <- abs(metlin)
-        digits <- ceiling(-log10(metlin))+1
-        metlinurl <-
-        paste("http://metlin.scripps.edu/simple_search_result.php?mass_min=",
-        round(neutralmass - metlin, digits), "&mass_max=",
-        round(neutralmass + metlin, digits), sep="")
-        values <- cbind(metlin = metlinurl, values)
-    }
-    twosamp <- cbind(name = groupnames(object), stat, groupmat, values)
-    if (sortpval) {
-        tsidx <- order(twosamp[,"pvalue"])
-        twosamp <- twosamp[tsidx,]
-        rownames(twosamp) <- 1:nrow(twosamp)
-        values<-values[tsidx,]
-    } else
-    tsidx <- 1:nrow(values)
-
-    if (length(filebase))
-    write.table(twosamp, paste(filebase, ".tsv", sep = ""), quote = FALSE, sep = "\t", col.names = NA)
-
-    if (eicmax > 0) {
-        if (length(unique(peaks(object)[,"rt"])) > 1) {
-            ## This looks like "normal" LC data
-
-            eicmax <- min(eicmax, length(tsidx))
-            eics <- getEIC(object, rtrange = eicwidth*1.1, sampleidx = ceic,
-            groupidx = tsidx[seq(length = eicmax)])
-
-            if (length(filebase)) {
-                eicdir <- paste(filebase, "_eic", sep="")
-                boxdir <- paste(filebase, "_box", sep="")
-                dir.create(eicdir)
-                dir.create(boxdir)
-                if (capabilities("png")){
-                    xcms:::xcmsBoxPlot(values[seq(length = eicmax),],
-                    sampclass(object), dirpath=boxdir, pic="png",  width=w, height=h)
-                    png(file.path(eicdir, "%003d.png"), width = w, height = h)
-                } else {
-                    xcms:::xcmsBoxPlot(values[seq(length = eicmax),],
-                    sampclass(object), dirpath=boxdir, pic="pdf", width=w, height=h)
-                    pdf(file.path(eicdir, "%003d.pdf"), width = w/72,
-                    height = h/72, onefile = FALSE)
-                }
-            }
-            plot(eics, object, rtrange = eicwidth, mzdec=mzdec)
-
-            if (length(filebase))
-            dev.off()
-        } else {
-            ## This looks like a direct-infusion single spectrum
-            if (length(filebase)) {
-                eicdir <- paste(filebase, "_eic", sep="")
-                boxdir <- paste(filebase, "_box", sep="")
-                dir.create(eicdir)
-                dir.create(boxdir)
-                if (capabilities("png")){
-                    xcmsBoxPlot(values[seq(length = eicmax),],
-                    sampclass(object), dirpath=boxdir, pic="png",
-                    width=w, height=h)
-                    png(file.path(eicdir, "%003d.png"), width = w, height = h,
-                    units = "px")
-                } else {
-                    xcmsBoxPlot(values[seq(length = eicmax),],
-                    sampclass(object), dirpath=boxdir, pic="pdf",
-                    width=w, height=h)
-                    pdf(file.path(eicdir, "%003d.pdf"), width = w/72,
-                    height = h/72, onefile = FALSE)
-                }
-            }
-
-            plotSpecWindow(object, gidxs = tsidx[seq(length = eicmax)], borderwidth=1)
-
-            if (length(filebase))
-            dev.off()
-        }
-    }
-
-    invisible(twosamp)
-}