# HG changeset patch # User iuc # Date 1441442265 14400 # Node ID 982bd8bfa44c4f86860cf9a36c020ae572fcf6c2 # Parent 4ba2afa87fce84e6f37d40420d6b98765eb14440 planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/deseq2 commit 470804fe1d65ce5bf49804c7f249c51f6a567f9f diff -r 4ba2afa87fce -r 982bd8bfa44c deseq2.xml --- a/deseq2.xml Fri Sep 04 08:26:02 2015 -0400 +++ b/deseq2.xml Sat Sep 05 04:37:45 2015 -0400 @@ -6,16 +6,29 @@ atlas deseq2 + + + + + /dev/null | grep -v -i "WARNING: ") + echo $(R --version | grep version | grep -v GNU)", DESeq2 version" $(R --vanilla --slave -e "library(DESeq2); cat(sessionInfo()\$otherPkgs\$DESeq2\$Version)" 2> /dev/null | grep -v -i "WARNING: ") ]]> - - - - - + threshold - if (sum(!use) == 0) { - h <- hist(res$pvalue, breaks=0:50/50, plot=FALSE) - barplot(height = h$counts, - col = "powderblue", space = 0, xlab="p-values", ylab="frequency", - main=paste("Histogram of p-values for",title_suffix)) - text(x = c(0, length(h1$counts)), y = 0, label=paste(c(0,1)), adj=c(0.5,1.7), xpd=NA) - } else { - h1 <- hist(res$pvalue[!use], breaks=0:50/50, plot=FALSE) - h2 <- hist(res$pvalue[use], breaks=0:50/50, plot=FALSE) - colori <- c("filtered (low count)"="khaki", "not filtered"="powderblue") - barplot(height = rbind(h1$counts, h2$counts), beside = FALSE, - col = colori, space = 0, xlab="p-values", ylab="frequency", - main=paste("Histogram of p-values for",title_suffix)) - text(x = c(0, length(h1$counts)), y = 0, label=paste(c(0,1)), adj=c(0.5,1.7), xpd=NA) - legend("topright", fill=rev(colori), legend=rev(names(colori)), bg="white") - } - plotMA(res, main= paste("MA-plot for",title_suffix), ylim=range(res$log2FoldChange, na.rm=TRUE)) -} - -if (verbose) { - cat(paste("primary factor:",primaryFactor,"\n")) - if (length(factors) > 1) { - cat(paste("other factors in design:",paste(factors[-length(factors)],collapse=","),"\n")) - } - cat("\n---------------------\n") -} - -# if JSON input from Galaxy, path is absolute -# otherwise, from sample_table, assume it is relative -dir <- if (is.null(opt$factors)) { - "." -} else { - "" -} - -# construct the object -dds <- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable, - directory = dir, - design = designFormula) - -if (verbose) cat(paste(ncol(dds), "samples with counts over", nrow(dds), "genes\n")) - -# optional outlier behavior -if (is.null(opt$outlier_replace_off)) { - minRep <- 7 -} else { - minRep <- Inf - if (verbose) cat("outlier replacement off\n") -} -if (is.null(opt$outlier_filter_off)) { - cooksCutoff <- TRUE -} else { - cooksCutoff <- FALSE - if (verbose) cat("outlier filtering off\n") -} - -# optional automatic mean filtering -if (is.null(opt$auto_mean_filter_off)) { - independentFiltering <- TRUE -} else { - independentFiltering <- FALSE - if (verbose) cat("automatic filtering on the mean off\n") -} - -# shrinkage of LFCs -if (is.null(opt$beta_prior_off)) { - betaPrior <- TRUE -} else { - betaPrior <- FALSE - if (verbose) cat("beta prior off\n") -} - -# dispersion fit type -if (is.null(opt$fit_type)) { - fitType <- "parametric" -} else { - fitType <- c("parametric","local","mean")[opt$fit_type] -} - -if (verbose) cat(paste("using disperion fit type:",fitType,"\n")) - -# run the analysis -dds <- DESeq(dds, fitType=fitType, betaPrior=betaPrior, minReplicatesForReplace=minRep) - -# create the generic plots and leave the device open -if (!is.null(opt$plots)) { - if (verbose) cat("creating plots\n") - pdf(opt$plots) - generateGenericPlots(dds, factors) -} - -n <- nlevels(colData(dds)[[primaryFactor]]) -allLevels <- levels(colData(dds)[[primaryFactor]]) - -if (is.null(opt$many_contrasts)) { - # only contrast the first and second level of the primary factor - ref <- allLevels[1] - lvl <- allLevels[2] - res <- results(dds, contrast=c(primaryFactor, lvl, ref), - cooksCutoff=cooksCutoff, - independentFiltering=independentFiltering) - if (verbose) { - cat("summary of results\n") - cat(paste0(primaryFactor,": ",lvl," vs ",ref,"\n")) - print(summary(res)) - } - resSorted <- res[order(res$padj),] - outDF <- as.data.frame(resSorted) - outDF$geneID <- rownames(outDF) - outDF <- outDF[,c("geneID", "baseMean", "log2FoldChange", "lfcSE", "stat", "pvalue", "padj")] - filename <- opt$outfile - write.table(outDF, file=filename, sep="\t", quote=FALSE, row.names=FALSE, col.names=FALSE) - if (independentFiltering) { - threshold <- unname(attr(res, "filterThreshold")) - } else { - threshold <- 0 - } - title_suffix <- paste0(primaryFactor,": ",lvl," vs ",ref) - if (!is.null(opt$plots)) { - generateSpecificPlots(res, threshold, title_suffix) - } -} else { - # rotate through the possible contrasts of the primary factor - # write out a sorted table of results with the contrast as a suffix - # add contrast specific plots to the device - for (i in seq_len(n-1)) { - ref <- allLevels[i] - contrastLevels <- allLevels[(i+1):n] - for (lvl in contrastLevels) { - res <- results(dds, contrast=c(primaryFactor, lvl, ref), - cooksCutoff=cooksCutoff, - independentFiltering=independentFiltering) - resSorted <- res[order(res$padj),] - outDF <- as.data.frame(resSorted) - outDF$geneID <- rownames(outDF) - outDF <- outDF[,c("geneID", "baseMean", "log2FoldChange", "lfcSE", "stat", "pvalue", "padj")] - filename <- paste0(opt$outfile,".",primaryFactor,"_",lvl,"_vs_",ref) - write.table(outDF, file=filename, sep="\t", quote=FALSE, row.names=FALSE, col.names=FALSE) - if (independentFiltering) { - threshold <- unname(attr(res, "filterThreshold")) - } else { - threshold <- 0 - } - title_suffix <- paste0(primaryFactor,": ",lvl," vs ",ref) - if (!is.null(opt$plots)) { - generateSpecificPlots(res, threshold, title_suffix) - } - } - } -} - -# close the plot device -if (!is.null(opt$plots)) { - cat("closing plot device\n") - dev.off() -} - -cat("Session information:\n\n") - -sessionInfo() - diff -r 4ba2afa87fce -r 982bd8bfa44c tool_dependencies.xml --- a/tool_dependencies.xml Fri Sep 04 08:26:02 2015 -0400 +++ b/tool_dependencies.xml Sat Sep 05 04:37:45 2015 -0400 @@ -4,6 +4,6 @@ - +