comparison limma_voom.R @ 11:7e8af58c8052 draft

planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/limma_voom commit 4bcbd83574ecf7194f3370aa883a9573563afdbc

10:e26047c4562d

    "robOpt", "b", 0, "logical",
    "trend", "t", 1, "double",
    "weightOpt", "w", 0, "logical",
    "topgenes", "G", 1, "integer",
    "treatOpt", "T", 0, "logical",
    "plots", "P", 1, "character"),
    byrow=TRUE, ncol=4)
opt <- getopt(spec)


if (is.null(opt$matrixPath) & is.null(opt$filesPath)) {

if (is.null(opt$treatOpt)) {
    wantTreat <- FALSE
} else {
    wantTreat <- TRUE
}


if (!is.null(opt$filesPath)) {
    # Process the separate count files (adapted from DESeq2 wrapper)

    countfiles <- lapply(sampleTable$filename, function(x){read.delim(x, row.names=1)})
    counts <- do.call("cbind", countfiles)

} else {
    # Process the single count matrix
    counts <- read.table(opt$matrixPath, header=TRUE, sep="\t", stringsAsFactors=FALSE)
    row.names(counts) <- counts[, 1]
    counts <- counts[ , -1]
    countsRows <- nrow(counts)

    # Process factors
    if (is.null(opt$factInput)) {
            factorData <- read.table(opt$factFile, header=TRUE, sep="\t")
            factors <- factorData[, -1, drop=FALSE]
    }  else {
            factors <- unlist(strsplit(opt$factInput, "|", fixed=TRUE))
            factorData <- list()
            for (fact in factors) {

    }
}

 # if annotation file provided
if (haveAnno) {
    geneanno <- read.table(opt$annoPath, header=TRUE, sep="\t", stringsAsFactors=FALSE)
}

#Create output directory
dir.create(opt$outPath, showWarnings=FALSE)

    colnames(design) <- gsub(factorList[i], "", colnames(design), fixed=TRUE)
}

# Calculating normalising factors
print("Calculating Normalisation Factors")
y <- calcNormFactors(y, method=opt$normOpt)

# Generate contrasts information
print("Generating Contrasts")
contrasts <- makeContrasts(contrasts=contrastData, levels=design)

    # Plot MD plots for individual samples
    print("Generating MD plots for samples")
    pdf(mdsamOutPdf, width=6.5, height=10)
    par(mfrow=c(3, 2))
    for (i in 1:nsamples) {
        plotMD(y, column = i)
        abline(h=0, col="red", lty=2, lwd=2)
    }
    linkName <- "MDPlots_Samples.pdf"
    linkAddr <- "mdplots_samples.pdf"
    linkData <- rbind(linkData, c(linkName, linkAddr))

        fit <- eBayes(voomFit, robust=FALSE)
    }
    plotData <- vData
}


print("Generating DE results")

if (wantTreat) {
    print("Applying TREAT method")

11:7e8af58c8052

    "robOpt", "b", 0, "logical",
    "trend", "t", 1, "double",
    "weightOpt", "w", 0, "logical",
    "topgenes", "G", 1, "integer",
    "treatOpt", "T", 0, "logical",
    "plots", "P", 1, "character",
    "libinfoOpt", "L", 0, "logical"),
    byrow=TRUE, ncol=4)
opt <- getopt(spec)


if (is.null(opt$matrixPath) & is.null(opt$filesPath)) {

if (is.null(opt$treatOpt)) {
    wantTreat <- FALSE
} else {
    wantTreat <- TRUE
}

if (is.null(opt$libinfoOpt)) {
    wantLibinfo <- FALSE
} else {
    wantLibinfo <- TRUE
}


if (!is.null(opt$filesPath)) {
    # Process the separate count files (adapted from DESeq2 wrapper)

    countfiles <- lapply(sampleTable$filename, function(x){read.delim(x, row.names=1)})
    counts <- do.call("cbind", countfiles)

} else {
    # Process the single count matrix
    counts <- read.table(opt$matrixPath, header=TRUE, sep="\t", strip.white=TRUE, stringsAsFactors=FALSE)
    row.names(counts) <- counts[, 1]
    counts <- counts[ , -1]
    countsRows <- nrow(counts)

    # Process factors
    if (is.null(opt$factInput)) {
            factorData <- read.table(opt$factFile, header=TRUE, sep="\t", strip.white=TRUE)
            factors <- factorData[, -1, drop=FALSE]
    }  else {
            factors <- unlist(strsplit(opt$factInput, "|", fixed=TRUE))
            factorData <- list()
            for (fact in factors) {

    }
}

 # if annotation file provided
if (haveAnno) {
    geneanno <- read.table(opt$annoPath, header=TRUE, sep="\t", quote= "", strip.white=TRUE, stringsAsFactors=FALSE)
}

#Create output directory
dir.create(opt$outPath, showWarnings=FALSE)

    colnames(design) <- gsub(factorList[i], "", colnames(design), fixed=TRUE)
}

# Calculating normalising factors
print("Calculating Normalisation Factors")
logcounts <- y #store for plots
y <- calcNormFactors(y, method=opt$normOpt)

# Generate contrasts information
print("Generating Contrasts")
contrasts <- makeContrasts(contrasts=contrastData, levels=design)

    # Plot MD plots for individual samples
    print("Generating MD plots for samples")
    pdf(mdsamOutPdf, width=6.5, height=10)
    par(mfrow=c(3, 2))
    for (i in 1:nsamples) {
        if (opt$normOpt != "none") {
            plotMD(logcounts, column=i, main=paste(colnames(logcounts)[i], "(before)"))
            abline(h=0, col="red", lty=2, lwd=2)
        }
        plotMD(y, column=i)
        abline(h=0, col="red", lty=2, lwd=2)
    }
    linkName <- "MDPlots_Samples.pdf"
    linkAddr <- "mdplots_samples.pdf"
    linkData <- rbind(linkData, c(linkName, linkAddr))

        fit <- eBayes(voomFit, robust=FALSE)
    }
    plotData <- vData
}

 # Save library size info
if (wantLibinfo) {
    efflibsize <- round(y$samples$lib.size * y$samples$norm.factors)
    libsizeinfo <- cbind(y$samples, EffectiveLibrarySize=efflibsize)
    libsizeinfo$lib.size <- round(libsizeinfo$lib.size)
    names(libsizeinfo)[names(libsizeinfo)=="sampleID"] <- "SampleID"
    names(libsizeinfo)[names(libsizeinfo)=="lib.size"] <- "LibrarySize"
    names(libsizeinfo)[names(libsizeinfo)=="norm.factors"] <- "NormalisationFactor"
    write.table(libsizeinfo, file="libsizeinfo", row.names=FALSE, sep="\t", quote=FALSE)
}

print("Generating DE results")

if (wantTreat) {
    print("Applying TREAT method")