comparison lib.r @ 36:b6a8143c397c draft

planemo upload commit aa538ae477bc2f901c95d916e11d70eab75b5e00

35:b02797fbead1

    system(paste0("gm convert ",filebase,"_box/*.png ",pdfBoxOutput))

}

#@author G. Le Corguille
#This function convert if it is required the Retention Time in minutes
RTSecondToMinute <- function(variableMetadata, convertRTMinute) {
    if (convertRTMinute){
        #converting the retention times (seconds) into minutes
        print("converting the retention times into minutes in the variableMetadata")

    variableMetadata=cbind(name=variableMetadata$name, namecustom=namecustom, variableMetadata[,!(colnames(variableMetadata) %in% c("name"))])
    return(variableMetadata)
}

#The function annotateDiffreport without the corr function which bugs
annotatediff <- function(xset=xset, listArguments=listArguments, variableMetadataOutput="variableMetadata.tsv", dataMatrixOutput="dataMatrix.tsv") {
    # Resolve the bug with x11, with the function png
    options(bitmapType='cairo')

    #Check if the fillpeaks step has been done previously, if it hasn't, there is an error message and the execution is stopped.
    res=try(is.null(xset@filled))

    # launch annotate
    xa = do.call("annotate", listArguments4annotate)
    peakList=getPeaklist(xa,intval=listArguments[["intval"]])
    peakList=cbind(groupnames(xa@xcmsSet),peakList); colnames(peakList)[1] = c("name");

    # --- dataMatrix ---
    dataMatrix = peakList[,(make.names(colnames(peakList)) %in% c("name", make.names(sampnames(xa@xcmsSet))))]
    write.table(dataMatrix, sep="\t", quote=FALSE, row.names=FALSE, file=dataMatrixOutput)


    # --- Multi condition : diffreport ---
    diffrepOri=NULL
    if (!is.null(listArguments[["runDiffreport"]]) & nlevels(sampclass(xset))>=2) {
        #Check if the fillpeaks step has been done previously, if it hasn't, there is an error message and the execution is stopped.

                    dir.create("tabular")
                    write.table(diffrep, sep="\t", quote=FALSE, row.names=FALSE, file=paste("tabular/",filebase,"_tsv.tabular",sep=""))

                    if (listArguments[["eicmax"]] != 0) {
                        diffreport_png2pdf(filebase)
                    }
                }
            }
        }
    }

    # --- variableMetadata ---
    variableMetadata=peakList[,!(make.names(colnames(peakList)) %in% c(make.names(sampnames(xa@xcmsSet))))]
    variableMetadata = RTSecondToMinute(variableMetadata, listArguments[["convertRTMinute"]])

############################################################
## 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(

    } 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))))

    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)

        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

                      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]

    ## 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.")

        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)

                    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)) {

                    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)
}

36:b6a8143c397c

    system(paste0("gm convert ",filebase,"_box/*.png ",pdfBoxOutput))

}

#@author G. Le Corguille
#The function create a zip archive from the different png generated by diffreport
diffreport_png2zip <- function() {
    zip("eic.zip",dir(pattern="_eic"))
    zip("box.zip",dir(pattern="_box"))
}

#The function create a zip archive from the different tabular generated by diffreport
diffreport_tabular2zip <- function() {
    zip("tabular.zip",dir(pattern="tabular/*"))
}

#@author G. Le Corguille
#This function convert if it is required the Retention Time in minutes
RTSecondToMinute <- function(variableMetadata, convertRTMinute) {
    if (convertRTMinute){
        #converting the retention times (seconds) into minutes
        print("converting the retention times into minutes in the variableMetadata")

    variableMetadata=cbind(name=variableMetadata$name, namecustom=namecustom, variableMetadata[,!(colnames(variableMetadata) %in% c("name"))])
    return(variableMetadata)
}

#The function annotateDiffreport without the corr function which bugs
annotatediff <- function(xset=xset, listArguments=listArguments, variableMetadataOutput="variableMetadata.tsv") {
    # Resolve the bug with x11, with the function png
    options(bitmapType='cairo')

    #Check if the fillpeaks step has been done previously, if it hasn't, there is an error message and the execution is stopped.
    res=try(is.null(xset@filled))

    # launch annotate
    xa = do.call("annotate", listArguments4annotate)
    peakList=getPeaklist(xa,intval=listArguments[["intval"]])
    peakList=cbind(groupnames(xa@xcmsSet),peakList); colnames(peakList)[1] = c("name");

    # --- Multi condition : diffreport ---
    diffrepOri=NULL
    if (!is.null(listArguments[["runDiffreport"]]) & nlevels(sampclass(xset))>=2) {
        #Check if the fillpeaks step has been done previously, if it hasn't, there is an error message and the execution is stopped.

                    dir.create("tabular")
                    write.table(diffrep, sep="\t", quote=FALSE, row.names=FALSE, file=paste("tabular/",filebase,"_tsv.tabular",sep=""))

                    if (listArguments[["eicmax"]] != 0) {
                        if (listArguments[["png2"]] == "pdf")
                            diffreport_png2pdf(filebase)
                    }
                }
            }
        }
        if (listArguments[["png2"]] == "zip")
            diffreport_png2zip()
        if (listArguments[["tabular2"]] == "zip")
            diffreport_tabular2zip()
    }

    # --- variableMetadata ---
    variableMetadata=peakList[,!(make.names(colnames(peakList)) %in% c(make.names(sampnames(xa@xcmsSet))))]
    variableMetadata = RTSecondToMinute(variableMetadata, listArguments[["convertRTMinute"]])

############################################################
## 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(

    } 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))))

    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)

        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

                      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]

    ## 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.")

        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)

                    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)) {

                    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)
}