| Previous changeset 2:ab492df30cdf (2017-10-23) Next changeset 4:ae39895af5fe (2018-09-30) |
|
Commit message:
planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/goseq commit 46c4278d292ab4d76dc5f3f74c3109c3179be7ef |
|
modified:
goseq.r goseq.xml |
|
added:
test-data/topgo.pdf |
|
removed:
test-data/getgo.danRer10.tab test-data/getgo.hg38.tab test-data/wal.tab |
| b |
| diff -r ab492df30cdf -r 783e8b70b047 goseq.r --- 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 |
| b |
| diff -r ab492df30cdf -r 783e8b70b047 goseq.xml --- a/goseq.xml Mon Oct 23 11:19:12 2017 -0400 +++ b/goseq.xml Mon Sep 24 06:29:03 2018 -0400 |
| [ |
| b'@@ -1,12 +1,14 @@\n-<tool id="goseq" name="goseq" version="1.26.0">\n+<tool id="goseq" name="goseq" version="1.30.0">\n <description>tests for overrepresented gene categories</description>\n <requirements>\n- <requirement type="package" version="1.3.2">r-optparse</requirement>\n- <requirement type="package" version="1.26.0">bioconductor-goseq</requirement>\n- <requirement type="package" version="3.3.0">bioconductor-org.hs.eg.db</requirement>\n- <requirement type="package" version="3.4.0">bioconductor-org.dm.eg.db</requirement>\n- <requirement type="package" version="3.4.1">bioconductor-org.dr.eg.db</requirement>\n- <requirement type="package" version="3.4.0">bioconductor-org.mm.eg.db</requirement>\n+ <requirement type="package" version="1.6.0">r-optparse</requirement>\n+ <requirement type="package" version="1.30.0">bioconductor-goseq</requirement>\n+ <requirement type="package" version="3.5.0">bioconductor-org.hs.eg.db</requirement>\n+ <requirement type="package" version="3.5.0">bioconductor-org.dm.eg.db</requirement>\n+ <requirement type="package" version="3.5.0">bioconductor-org.dr.eg.db</requirement>\n+ <requirement type="package" version="3.5.0">bioconductor-org.mm.eg.db</requirement>\n+ <requirement type="package" version="0.7.6">r-dplyr</requirement>\n+ <requirement type="package" version="3.0.0">r-ggplot2</requirement>\n </requirements>\n <stdio>\n <regex match="Execution halted"\n@@ -23,7 +25,7 @@\n description="An undefined error occured, please check your input carefully and contact your administrator." />\n </stdio>\n <version_command><![CDATA[\n-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: ")\n+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 |'..b'GO terms plot">\n+ <filter>methods[\'wallenius\']</filter>\n+ <filter>out[\'topgo_plot\']</filter>\n+ </data>\n </outputs>\n \n <tests>\n+ <!-- Ensure top plot is output -->\n+ <test expect_num_outputs="2">\n+ <param name="dge_file" value="dge_list.tab" ftype="tabular" />\n+ <param name="length_file" value="gene_length.tab" ftype="tabular" />\n+ <param name="catSource" value="history" />\n+ <param name="category_file" value="category.tab" ftype="tabular" />\n+ <param name="use_genes_without_cat" value="true" />\n+ <param name="topgo_plot" value="true" />\n+ <output name="top_plot" ftype="pdf" file="topgo.pdf" compare="sim_size"/>\n+ </test>\n <!-- Ensure Wallenius table is output -->\n <test expect_num_outputs="1">\n <param name="dge_file" value="dge_list.tab" ftype="tabular" />\n@@ -149,7 +170,12 @@\n <param name="catSource" value="history" />\n <param name="category_file" value="category.tab" ftype="tabular" />\n <param name="use_genes_without_cat" value="true" />\n- <output name="wallenius_tab" file="wal.tab" compare="contains" />\n+ <output name="wallenius_tab">\n+ <assert_contents>\n+ <has_text_matching expression="category.*over_represented_pvalue.*under_represented_pvalue.*numDEInCat.*numInCat.*term.*ontology.*p.adjust.over_represented.*p.adjust.under_represented" />\n+ <has_text_matching expression="GO:0000278.*0.01" />\n+ </assert_contents>\n+ </output>\n </test>\n <!-- Ensure getting GO categories works -->\n <test expect_num_outputs="1">\n@@ -159,7 +185,12 @@\n <param name="genome" value="hg38" />\n <param name="gene_id" value="ensGene" />\n <param name="use_genes_without_cat" value="true" />\n- <output name="wallenius_tab" ftype="tabular" file="getgo.hg38.tab" compare="contains"/>\n+ <output name="wallenius_tab">\n+ <assert_contents>\n+ <has_text_matching expression="category.*over_represented_pvalue.*under_represented_pvalue.*numDEInCat.*numInCat.*term.*ontology.*p.adjust.over_represented.*p.adjust.under_represented" />\n+ <has_text_matching expression="GO:0005576.*5.2" />\n+ </assert_contents>\n+ </output>\n </test>\n <!-- Ensure getting GO categories for another genome (zebrafish) works -->\n <test expect_num_outputs="1">\n@@ -169,7 +200,12 @@\n <param name="genome" value="danRer10"/>\n <param name="gene_id" value="ensGene" />\n <param name="use_genes_without_cat" value="true" />\n- <output name="wallenius_tab" ftype="tabular" file="getgo.danRer10.tab" compare="contains"/>\n+ <output name="wallenius_tab">\n+ <assert_contents>\n+ <has_text_matching expression="category.*over_represented_pvalue.*under_represented_pvalue.*numDEInCat.*numInCat.*term.*ontology.*p.adjust.over_represented.*p.adjust.under_represented" />\n+ <has_text_matching expression="GO:0031324.*0.08" />\n+ </assert_contents>\n+ </output>\n </test>\n <!-- Ensure length bias plot works -->\n <test expect_num_outputs="2">\n@@ -277,7 +313,7 @@\n **Outputs**\n \n * 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.\n-* Optionally, this tool can also output some diagnostic plots and an RData file, see **Output Options** above.\n+* 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.\n \n Example:\n \n' |
| b |
| diff -r ab492df30cdf -r 783e8b70b047 test-data/getgo.danRer10.tab --- a/test-data/getgo.danRer10.tab Mon Oct 23 11:19:12 2017 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| b |
| @@ -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 |
| b |
| diff -r ab492df30cdf -r 783e8b70b047 test-data/getgo.hg38.tab --- a/test-data/getgo.hg38.tab Mon Oct 23 11:19:12 2017 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| b |
| @@ -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 |
| b |
| diff -r ab492df30cdf -r 783e8b70b047 test-data/topgo.pdf |
| b |
| Binary file test-data/topgo.pdf has changed |
| b |
| diff -r ab492df30cdf -r 783e8b70b047 test-data/wal.tab --- a/test-data/wal.tab Mon Oct 23 11:19:12 2017 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| b |
| @@ -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 |