changeset 3:783e8b70b047 draft

planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/goseq commit 46c4278d292ab4d76dc5f3f74c3109c3179be7ef

--- a/goseq.r	Mon Oct 23 11:19:12 2017 -0400
+++ b/goseq.r	Mon Sep 24 06:29:03 2018 -0400
@@ -6,6 +6,8 @@
 suppressPackageStartupMessages({
     library("goseq")
     library("optparse")
+    library("dplyr")
+    library("ggplot2")
 })
 
 option_list <- list(
@@ -25,7 +27,8 @@
                 help="A large number of gene may have no GO term annotated. If this option is set to FALSE, genes without category will be ignored in the calculation of p-values(default behaviour). If TRUE these genes will count towards the total number of genes outside the tested category (default behaviour prior to version 1.15.2)."),
     make_option(c("-plots", "--make_plots"), default=FALSE, type="logical", help="produce diagnostic plots?"),
     make_option(c("-fc", "--fetch_cats"), default=NULL, type="character", help="Categories to get can include one or more of GO:CC, GO:BP, GO:MF, KEGG"),
-    make_option(c("-rd", "--rdata"), default=NULL, type="character", help="Path to RData output file.")
+    make_option(c("-rd", "--rdata"), default=NULL, type="character", help="Path to RData output file."),
+    make_option(c("-tp", "--top_plot"), default=NULL, type="logical", help="Output PDF with top10 over-rep GO terms?")
     )
 
 parser <- OptionParser(usage = "%prog [options] file", option_list=option_list)
@@ -37,9 +40,7 @@
 length_file = args$length_file
 genome = args$genome
 gene_id = args$gene_id
-wallenius_tab = args$wallenius_tab
 sampling_tab = args$sampling_tab
-nobias_tab = args$nobias_tab
 length_bias_plot = args$length_bias_plot
 sample_vs_wallenius_plot = args$sample_vs_wallenius_plot
 repcnt = args$repcnt
@@ -50,6 +51,8 @@
 
 if (!is.null(args$fetch_cats)) {
   fetch_cats = unlist(strsplit(args$fetch_cats, ","))
+} else {
+  fetch_cats = "Custom"
 }
 
 # format DE genes into named vector suitable for goseq
@@ -84,7 +87,9 @@
   pdf(length_bias_plot)
 }
 pwf=nullp(genes, genome = genome, id = gene_id, bias.data = gene_lengths, plot.fit=make_plots)
-graphics.off()
+if (make_plots != 'false') {
+  dev.off()
+}
 
 # Fetch GO annotations if category_file hasn't been supplied:
 if (category_file == "FALSE") {
@@ -99,20 +104,24 @@
     }
 }
 
+results <- list()
+
 # wallenius approximation of p-values
-if (wallenius_tab != "" && wallenius_tab!="None") {
+if (!is.null(args$wallenius_tab)) {
   GO.wall=goseq(pwf, genome = genome, id = gene_id, use_genes_without_cat = use_genes_without_cat, gene2cat=go_map)
   GO.wall$p.adjust.over_represented = p.adjust(GO.wall$over_represented_pvalue, method=p_adj_method)
   GO.wall$p.adjust.under_represented = p.adjust(GO.wall$under_represented_pvalue, method=p_adj_method)
-  write.table(GO.wall, wallenius_tab, sep="\t", row.names = FALSE, quote = FALSE)
+  write.table(GO.wall, args$wallenius_tab, sep="\t", row.names = FALSE, quote = FALSE)
+  results[['Wallenius']] <- GO.wall
 }
 
 # hypergeometric (no length bias correction)
-if (nobias_tab != "" && nobias_tab != "None") {
+if (!is.null(args$nobias_tab)) {
   GO.nobias=goseq(pwf, genome = genome, id = gene_id, method="Hypergeometric", use_genes_without_cat = use_genes_without_cat, gene2cat=go_map)
   GO.nobias$p.adjust.over_represented = p.adjust(GO.nobias$over_represented_pvalue, method=p_adj_method)
   GO.nobias$p.adjust.under_represented = p.adjust(GO.nobias$under_represented_pvalue, method=p_adj_method)
-  write.table(GO.nobias, nobias_tab, sep="\t", row.names = FALSE, quote = FALSE)
+  write.table(GO.nobias, args$nobias_tab, sep="\t", row.names = FALSE, quote = FALSE)
+  results[['Hypergeometric']] <- GO.nobias
 }
 
 # Sampling distribution
@@ -134,8 +143,29 @@
      xlab="log10(Wallenius p-values)",ylab="log10(Sampling p-values)",
      xlim=c(-3,0))
      abline(0,1,col=3,lty=2)
-  graphics.off()
+  dev.off()
   }
+  results[['Sampling']] <- GO.samp
+}
+
+if (!is.null(args$top_plot)) {
+  # modified from https://bioinformatics-core-shared-training.github.io/cruk-summer-school-2018/RNASeq2018/html/06_Gene_set_testing.nb.html
+  pdf("top10.pdf")
+  for (m in names(results)) {
+    p <- results[[m]] %>%
+      top_n(10, wt=-p.adjust.over_represented)  %>%
+      mutate(hitsPerc=numDEInCat*100/numInCat) %>%
+      ggplot(aes(x=hitsPerc,
+                   y=term,
+                   colour=p.adjust.over_represented,
+                   size=numDEInCat)) +
+      geom_point() +
+      expand_limits(x=0) +
+      labs(x="% DE in category", y="Category", colour="adj. P value", size="Count", title=paste("Top over-represented categories in", fetch_cats), subtitle=paste(m, " method")) +
+      theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5))
+    print(p)
+  }
+  dev.off()
 }
 
 # Output RData file

--- a/goseq.xml	Mon Oct 23 11:19:12 2017 -0400
+++ b/goseq.xml	Mon Sep 24 06:29:03 2018 -0400
@@ -1,12 +1,14 @@
-<tool id="goseq" name="goseq" version="1.26.0">
+<tool id="goseq" name="goseq" version="1.30.0">
     <description>tests for overrepresented gene categories</description>
     <requirements>
-        <requirement type="package" version="1.3.2">r-optparse</requirement>
-        <requirement type="package" version="1.26.0">bioconductor-goseq</requirement>
-        <requirement type="package" version="3.3.0">bioconductor-org.hs.eg.db</requirement>
-        <requirement type="package" version="3.4.0">bioconductor-org.dm.eg.db</requirement>
-        <requirement type="package" version="3.4.1">bioconductor-org.dr.eg.db</requirement>
-        <requirement type="package" version="3.4.0">bioconductor-org.mm.eg.db</requirement>
+        <requirement type="package" version="1.6.0">r-optparse</requirement>
+        <requirement type="package" version="1.30.0">bioconductor-goseq</requirement>
+        <requirement type="package" version="3.5.0">bioconductor-org.hs.eg.db</requirement>
+        <requirement type="package" version="3.5.0">bioconductor-org.dm.eg.db</requirement>
+        <requirement type="package" version="3.5.0">bioconductor-org.dr.eg.db</requirement>
+        <requirement type="package" version="3.5.0">bioconductor-org.mm.eg.db</requirement>
+        <requirement type="package" version="0.7.6">r-dplyr</requirement>
+        <requirement type="package" version="3.0.0">r-ggplot2</requirement>
     </requirements>
     <stdio>
         <regex match="Execution halted"
@@ -23,7 +25,7 @@
                description="An undefined error occured, please check your input carefully and contact your administrator." />
     </stdio>
     <version_command><![CDATA[
-echo $(R --version | grep version | grep -v GNU)", goseq version" $(R --vanilla --slave -e "library(goseq); cat(sessionInfo()\$otherPkgs\$goseq\$Version)" 2> /dev/null | grep -v -i "WARNING: ")", optparse version" $(R --vanilla --slave -e "library(optparse); cat(sessionInfo()\$otherPkgs\$optparse\$Version)" 2> /dev/null | grep -v -i "WARNING: ")", org.Hs.eg.db version" $(R --vanilla --slave -e "library(org.Hs.eg.db); cat(sessionInfo()\$otherPkgs\$org.Hs.eg.db\$Version)" 2> /dev/null | grep -v -i "WARNING: ")", org.Dr.eg.db version" $(R --vanilla --slave -e "library(org.Dr.eg.db); cat(sessionInfo()\$otherPkgs\$org.Dr.eg.db\$Version)" 2> /dev/null | grep -v -i "WARNING: ")", org.Dm.eg.db version" $(R --vanilla --slave -e "library(org.Dm.eg.db); cat(sessionInfo()\$otherPkgs\$org.Dm.eg.db\$Version)" 2> /dev/null | grep -v -i "WARNING: ")", org.Mm.eg.db version" $(R --vanilla --slave -e "library(org.Mm.eg.db); cat(sessionInfo()\$otherPkgs\$org.Mm.eg.db\$Version)" 2> /dev/null | grep -v -i "WARNING: ")
+echo $(R --version | grep version | grep -v GNU)", goseq version" $(R --vanilla --slave -e "library(goseq); cat(sessionInfo()\$otherPkgs\$goseq\$Version)" 2> /dev/null | grep -v -i "WARNING: ")", optparse version" $(R --vanilla --slave -e "library(optparse); cat(sessionInfo()\$otherPkgs\$optparse\$Version)" 2> /dev/null | grep -v -i "WARNING: ")", org.Hs.eg.db version" $(R --vanilla --slave -e "library(org.Hs.eg.db); cat(sessionInfo()\$otherPkgs\$org.Hs.eg.db\$Version)" 2> /dev/null | grep -v -i "WARNING: ")", org.Dr.eg.db version" $(R --vanilla --slave -e "library(org.Dr.eg.db); cat(sessionInfo()\$otherPkgs\$org.Dr.eg.db\$Version)" 2> /dev/null | grep -v -i "WARNING: ")", org.Dm.eg.db version" $(R --vanilla --slave -e "library(org.Dm.eg.db); cat(sessionInfo()\$otherPkgs\$org.Dm.eg.db\$Version)" 2> /dev/null | grep -v -i "WARNING: ")", org.Mm.eg.db version" $(R --vanilla --slave -e "library(org.Mm.eg.db); cat(sessionInfo()\$otherPkgs\$org.Mm.eg.db\$Version)" 2> /dev/null | grep -v -i "WARNING: ")", dplyr version" $(R --vanilla --slave -e "library(dplyr); cat(sessionInfo()\$otherPkgs\$dplyr\$Version)" 2> /dev/null | grep -v -i "WARNING: ")", ggplot2 version" $(R --vanilla --slave -e "library(ggplot2); cat(sessionInfo()\$otherPkgs\$ggplot2\$Version)" 2> /dev/null | grep -v -i "WARNING: ")
     ]]></version_command>
     <command><![CDATA[
 Rscript '$__tool_directory__/goseq.r'
@@ -51,11 +53,15 @@
 --make_plots '$out.make_plots'
 --length_bias_plot '$length_bias_plot'
 --sample_vs_wallenius_plot '$sample_vs_wallenius_plot'
+
 --rdata '$out.rdata_out'
-
 --p_adj_method '$adv.p_adj_method'
 --use_genes_without_cat '$adv.use_genes_without_cat'
 
+#if $out.topgo_plot:
+    --top_plot '$out.topgo_plot'
+#end if
+
     ]]></command>
 
     <!-- Input Files-->
@@ -100,6 +106,7 @@
 
         <!-- Output Options -->
         <section name="out" title="Output Options">
+            <param name="topgo_plot" type="boolean" checked="false" label="Output Top GO terms plot?" help="Output a PDF plot of the Top 10 over-represented GO terms. Default: No" />
             <param name="make_plots" type="boolean" checked="false" label="Produce diagnostic plots?" help="This will produce the length bias (PWF) plot. If both sampling and wallenius methods are selected, it will also produce a plot comparing their p-values. These plots may help you compare the different p-value estimation methods that goseq can use. Default: No" />
             <param name="rdata_out" type="boolean" checked="false" label="Output RData file?" help="Output all the data used by R to construct the tables and plots, can be loaded into R. Default: No" />
         </section>
@@ -120,7 +127,7 @@
 
     <outputs>
         <data name="wallenius_tab" format="tabular" label="${tool.name} on ${on_string}: Ranked category list - Wallenius method">
-            <filter>methods['wallenius]'</filter>
+            <filter>methods['wallenius']</filter>
         </data>
         <data name="sampling_tab" format="tabular" label="${tool.name} on ${on_string}: Ranked category list - Sampling method">
             <filter>methods['repcnt'] != 0</filter>
@@ -139,9 +146,23 @@
         <data name="rdata" format="rdata" from_work_dir="goseq_analysis.RData" label="${tool.name} on ${on_string}: RData file">
             <filter>out['rdata_out']</filter>
         </data>
+        <data name="top_plot" format="pdf" from_work_dir="top10.pdf" label="${tool.name} on ${on_string}: Top over-represented GO terms plot">
+            <filter>methods['wallenius']</filter>
+            <filter>out['topgo_plot']</filter>
+        </data>
     </outputs>
 
     <tests>
+        <!-- Ensure top plot is output -->
+        <test expect_num_outputs="2">
+            <param name="dge_file" value="dge_list.tab" ftype="tabular" />
+            <param name="length_file" value="gene_length.tab" ftype="tabular" />
+            <param name="catSource" value="history" />
+            <param name="category_file" value="category.tab" ftype="tabular" />
+            <param name="use_genes_without_cat" value="true" />
+            <param name="topgo_plot" value="true" />
+            <output name="top_plot" ftype="pdf" file="topgo.pdf" compare="sim_size"/>
+        </test>
         <!-- Ensure Wallenius table is output -->
         <test expect_num_outputs="1">
             <param name="dge_file" value="dge_list.tab" ftype="tabular" />
@@ -149,7 +170,12 @@
             <param name="catSource" value="history" />
             <param name="category_file" value="category.tab" ftype="tabular" />
             <param name="use_genes_without_cat" value="true" />
-            <output name="wallenius_tab" file="wal.tab" compare="contains" />
+            <output name="wallenius_tab">
+                <assert_contents>
+                    <has_text_matching expression="category.*over_represented_pvalue.*under_represented_pvalue.*numDEInCat.*numInCat.*term.*ontology.*p.adjust.over_represented.*p.adjust.under_represented" />
+                    <has_text_matching expression="GO:0000278.*0.01" />
+                </assert_contents>
+            </output>
         </test>
         <!-- Ensure getting GO categories works -->
         <test expect_num_outputs="1">
@@ -159,7 +185,12 @@
             <param name="genome" value="hg38" />
             <param name="gene_id" value="ensGene" />
             <param name="use_genes_without_cat" value="true" />
-            <output name="wallenius_tab" ftype="tabular" file="getgo.hg38.tab" compare="contains"/>
+            <output name="wallenius_tab">
+                <assert_contents>
+                    <has_text_matching expression="category.*over_represented_pvalue.*under_represented_pvalue.*numDEInCat.*numInCat.*term.*ontology.*p.adjust.over_represented.*p.adjust.under_represented" />
+                    <has_text_matching expression="GO:0005576.*5.2" />
+                </assert_contents>
+            </output>
         </test>
         <!-- Ensure getting GO categories for another genome (zebrafish) works -->
         <test expect_num_outputs="1">
@@ -169,7 +200,12 @@
             <param name="genome" value="danRer10"/>
             <param name="gene_id" value="ensGene" />
             <param name="use_genes_without_cat" value="true" />
-            <output name="wallenius_tab" ftype="tabular" file="getgo.danRer10.tab" compare="contains"/>
+            <output name="wallenius_tab">
+                <assert_contents>
+                    <has_text_matching expression="category.*over_represented_pvalue.*under_represented_pvalue.*numDEInCat.*numInCat.*term.*ontology.*p.adjust.over_represented.*p.adjust.under_represented" />
+                    <has_text_matching expression="GO:0031324.*0.08" />
+                </assert_contents>
+            </output>
         </test>
         <!-- Ensure length bias plot works -->
         <test expect_num_outputs="2">
@@ -277,7 +313,7 @@
 **Outputs**
 
 * This tool outputs a tabular file containing a ranked list of gene categories, similar to below. The default output is the Wallenius method table. If the Sampling and/or Hypergeometric methods are also selected, additional tables are produced.
-* Optionally, this tool can also output some diagnostic plots and an RData file, see **Output Options** above.
+* Optionally, this tool can also output a plot of the top 10 over-represented categories, some diagnostic plots and an RData file, see **Output Options** above.
 
 Example:
 

--- a/test-data/getgo.danRer10.tab	Mon Oct 23 11:19:12 2017 -0400
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,10 +0,0 @@
-category	over_represented_pvalue	under_represented_pvalue	numDEInCat	numInCat	term	ontology	p.adjust.over_represented	p.adjust.under_represented
-GO:0031324	0.019289727841568	0.997004018824821	6	9	negative regulation of cellular metabolic process	BP	1	1
-GO:0040011	0.0219399815699082	0.993554925323586	10	19	locomotion	BP	1	1
-GO:0048738	0.0232122438335162	1	3	3	cardiac muscle tissue development	BP	1	1
-GO:0031101	0.0232122438335199	1	3	3	fin regeneration	BP	1	1
-GO:0042246	0.0232122438335199	1	3	3	tissue regeneration	BP	1	1
-GO:0007050	0.023212243833521	1	3	3	cell cycle arrest	BP	1	1
-GO:0019783	0.0254384360641003	0.998148600664743	4	5	ubiquitin-like protein-specific protease activity	MF	1	1
-GO:0036459	0.0254384360641003	0.998148600664743	4	5	thiol-dependent ubiquitinyl hydrolase activity	MF	1	1
-GO:0101005	0.0254384360641003	0.998148600664743	4	5	ubiquitinyl hydrolase activity	MF	1	1

--- a/test-data/getgo.hg38.tab	Mon Oct 23 11:19:12 2017 -0400
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,10 +0,0 @@
-category	over_represented_pvalue	under_represented_pvalue	numDEInCat	numInCat	term	ontology	p.adjust.over_represented	p.adjust.under_represented
-GO:0005576	4.72734295222294e-05	0.999979271555286	56	142	extracellular region	CC	0.329456825863645	1
-GO:0005840	0.000150633625443482	0.999987765310632	9	12	ribosome	CC	0.329456825863645	1
-GO:0044763	0.000210237360853053	0.999883100939053	148	473	single-organism cellular process	BP	0.329456825863645	1
-GO:0044699	0.000229197548055812	0.999873090122854	158	513	single-organism process	BP	0.329456825863645	1
-GO:0065010	0.000394294879818402	0.999824474827037	43	108	extracellular membrane-bounded organelle	CC	0.329456825863645	1
-GO:0070062	0.000394294879818402	0.999824474827037	43	108	extracellular exosome	CC	0.329456825863645	1
-GO:0008150	0.000409074003076654	0.999785179807024	191	656	biological_process	BP	0.329456825863645	1
-GO:0005488	0.000447980265756431	0.99975072864471	175	589	binding	MF	0.329456825863645	1
-GO:0005198	0.000511195682086445	0.999905085898726	13	21	structural molecule activity	MF	0.329456825863645	1

Binary file test-data/topgo.pdf has changed

--- a/test-data/wal.tab	Mon Oct 23 11:19:12 2017 -0400
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,3 +0,0 @@
-category	over_represented_pvalue	under_represented_pvalue	numDEInCat	numInCat	term	ontology	p.adjust.over_represented	p.adjust.under_represented
-GO:0000278	0.0112350612534339	0.999376653834006	4	5	mitotic cell cycle	BP	0.0224701225068678	0.999376653834006
-GO:0000003	1	0.805913166914892	0	1	reproduction	BP	1	0.999376653834006