Mercurial > repos > mvdbeek > r_goseq_1_22_0
diff goseq.r @ 18:5fb82111ec62 draft
planemo upload for repository https://github.com/ARTbio/tools-artbio/tree/master/tools/goseq_1_22_0 commit fdd0811efc61c31f88ff17096fbe8ee8cfacd766-dirty
author | mvdbeek |
---|---|
date | Fri, 26 Feb 2016 11:42:43 -0500 |
parents | 1b03f6232900 |
children | 9442d1bf6d93 |
line wrap: on
line diff
--- a/goseq.r Fri Feb 26 08:49:59 2016 -0500 +++ b/goseq.r Fri Feb 26 11:42:43 2016 -0500 @@ -1,6 +1,6 @@ sink(stdout(), type = "message") -library(goseq) -library(optparse) +suppressWarnings(suppressMessages(library(goseq))) +suppressWarnings(suppressMessages(library(optparse))) option_list <- list( make_option(c("-d", "--dge_file"), type="character", help="Path to file with differential gene expression result"), @@ -34,15 +34,13 @@ sample_vs_wallenius_plot = args$sample_vs_wallenius_plot repcnt = args$repcnt - # format DE genes into vector suitable for use with goseq dge_table = read.delim(dge_file, header = TRUE, sep="\t", check.names = FALSE) genes = as.integer(dge_table[,p_adj_column]<p_adj_cutoff) names(genes) = dge_table[,1] # Assuming first row contains gene names # Get gene lengths - -if (length_file) { +if (length_file !=FALSE) { length_table = read.delim(length_file, header=TRUE, sep="\t", check.names=FALSE, row.names=TRUE) gene_lengths = length_table[names(genes),]$length } else { @@ -53,7 +51,7 @@ pdf(length_bias_plot) pwf=nullp(genes, genome, gene_id, gene_lengths) -dev.off() +message = dev.off() # wallenius approximation of p-values GO.wall=goseq(pwf, genome, gene_id) @@ -68,7 +66,7 @@ xlab="log10(Wallenius p-values)",ylab="log10(Sampling p-values)", xlim=c(-3,0)) abline(0,1,col=3,lty=2) - dev.off() + message = dev.off() write.table(GO.samp, sampling_tab, sep="\t", row.names = FALSE, quote = FALSE) }