# HG changeset patch # User iuc # Date 1567770646 14400 # Node ID 8b3e3657034e58a6ab724c6941c72f259e6ce899 # Parent 67c29afac85fe0f1a6f9642f7da1436426e1aeec "planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/goseq commit 8e19f8bcaea6f607a1eaa14bb88f2d625ed63df0" diff -r 67c29afac85f -r 8b3e3657034e goseq.r --- a/goseq.r Sun Mar 17 10:27:17 2019 -0400 +++ b/goseq.r Fri Sep 06 07:50:46 2019 -0400 @@ -10,47 +10,32 @@ library("ggplot2") }) +sessionInfo() + option_list <- list( make_option(c("-d", "--dge_file"), type="character", help="Path to file with differential gene expression result"), - make_option(c("-w","--wallenius_tab"), type="character", help="Path to output file with P-values estimated using wallenius distribution."), - make_option(c("-s","--sampling_tab"), type="character", default=FALSE, help="Path to output file with P-values estimated using sampling distribution."), - make_option(c("-n","--nobias_tab"), type="character", default=FALSE, help="Path to output file with P-values estimated using hypergeometric distribution and no correction for gene length bias."), - make_option(c("-l","--length_bias_plot"), type="character", default=FALSE, help="Path to length-bias plot."), - make_option(c("-sw","--sample_vs_wallenius_plot"), type="character", default=FALSE, help="Path to plot comparing sampling with wallenius p-values."), - make_option(c("-r", "--repcnt"), type="integer", default=100, help="Number of repeats for sampling"), - make_option(c("-lf", "--length_file"), type="character", default="FALSE", help = "Path to tabular file mapping gene id to length"), - make_option(c("-cat_file", "--category_file"), default="FALSE", type="character", help = "Path to tabular file with gene_id <-> category mapping."), - make_option(c("-g", "--genome"), default=NULL, type="character", help = "Genome [used for looking up correct gene length]"), - make_option(c("-i", "--gene_id"), default=NULL, type="character", help = "Gene ID format of genes in DGE file"), - make_option(c("-p", "--p_adj_method"), default="BH", type="character", help="Multiple hypothesis testing correction method to use"), - make_option(c("-cat", "--use_genes_without_cat"), default=FALSE, type="logical", - 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("-tp", "--top_plot"), default=NULL, type="logical", help="Output PDF with top10 over-rep GO terms?") + make_option(c("-lf", "--length_file"), type="character", default=NULL, help="Path to tabular file mapping gene id to length"), + make_option(c("-g", "--genome"), type="character", default=NULL, help="Genome [used for looking up correct gene length]"), + make_option(c("-i", "--gene_id"), type="character", default=NULL, help="Gene ID format of genes in DGE file"), + make_option(c("-fc", "--fetch_cats"), type="character", default=NULL, help="Categories to get can include one or more of GO:CC, GO:BP, GO:MF, KEGG"), + make_option(c("-cat_file", "--category_file"), type="character", default=NULL, help="Path to tabular file with gene_id <-> category mapping"), + make_option(c("-w","--wallenius_tab"), type="character", default=NULL, help="Path to output file with P-values estimated using wallenius distribution"), + make_option(c("-n","--nobias_tab"), type="character", default=NULL, help="Path to output file with P-values estimated using hypergeometric distribution and no correction for gene length bias"), + make_option(c("-r", "--repcnt"), type="integer", default=0, help="Number of repeats for sampling"), + make_option(c("-s","--sampling_tab"), type="character", default=NULL, help="Path to output file with P-values estimated using sampling distribution"), + make_option(c("-p", "--p_adj_method"), type="character", default="BH", help="Multiple hypothesis testing correction method to use"), + make_option(c("-cat", "--use_genes_without_cat"), type="logical", default=FALSE, 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("-tp", "--top_plot"), type="character", default=NULL, help="Path to output PDF with top10 over-rep GO terms"), + make_option(c("-plots", "--make_plots"), default=FALSE, type="logical", help="Produce diagnostic plots?"), + make_option(c("-l","--length_bias_plot"), type="character", default=NULL, help="Path to length-bias plot"), + make_option(c("-sw","--sample_vs_wallenius_plot"), type="character", default=NULL, help="Path to plot comparing sampling with wallenius p-values"), + make_option(c("-rd", "--rdata"), type="character", default=NULL, help="Path to RData output file"), + make_option(c("-g2g", "--categories_genes_out_fp"), type="character", default=NULL, help="Path to file with categories (GO/KEGG terms) and associated DE genes") ) parser <- OptionParser(usage = "%prog [options] file", option_list=option_list) args = parse_args(parser) -# Vars: -dge_file = args$dge_file -category_file = args$category_file -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 -p_adj_method = args$p_adj_method -use_genes_without_cat = args$use_genes_without_cat -make_plots = args$make_plots -rdata = args$rdata - if (!is.null(args$fetch_cats)) { fetch_cats = unlist(strsplit(args$fetch_cats, ",")) } else { @@ -59,101 +44,115 @@ # format DE genes into named vector suitable for goseq # check if header is present -first_line = read.delim(dge_file, header = FALSE, nrow=1) +first_line = read.delim(args$dge_file, header=FALSE, nrow=1) second_col = toupper(first_line[, ncol(first_line)]) if (second_col == TRUE || second_col == FALSE) { - dge_table = read.delim(dge_file, header = FALSE, sep="\t") + dge_table = read.delim(args$dge_file, header=FALSE, sep="\t") } else { - dge_table = read.delim(dge_file, header = TRUE, sep="\t") + dge_table = read.delim(args$dge_file, header=TRUE, sep="\t") } -genes = as.numeric(as.logical(dge_table[,ncol(dge_table)])) # Last column contains TRUE/FALSE +genes = as.numeric(as.logical(dge_table[, ncol(dge_table)])) # Last column contains TRUE/FALSE names(genes) = dge_table[,1] # Assuming first column contains gene names # gene lengths, assuming last column -if (length_file != "FALSE" ) { - first_line = read.delim(length_file, header = FALSE, nrow=1) - if (is.numeric(first_line[, ncol(first_line)])) { - length_table = read.delim(length_file, header=FALSE, sep="\t", check.names=FALSE) - } else { - length_table = read.delim(length_file, header=TRUE, sep="\t", check.names=FALSE) - } - row.names(length_table) = length_table[,1] - gene_lengths = length_table[names(genes),][,ncol(length_table)] - } else { - gene_lengths = getlength(names(genes), genome, gene_id) - } +first_line = read.delim(args$length_file, header=FALSE, nrow=1) +if (is.numeric(first_line[, ncol(first_line)])) { + length_table = read.delim(args$length_file, header=FALSE, sep="\t", check.names=FALSE) +} else { + length_table = read.delim(args$length_file, header=TRUE, sep="\t", check.names=FALSE) +} +row.names(length_table) = length_table[,1] +# get vector of gene length in same order as the genes +gene_lengths = length_table[names(genes),][, ncol(length_table)] # Estimate PWF - -if (make_plots != 'false') { - pdf(length_bias_plot) +if (args$make_plots) { + pdf(args$length_bias_plot) } -pwf=nullp(genes, genome = genome, id = gene_id, bias.data = gene_lengths, plot.fit=make_plots) -if (make_plots != 'false') { +pwf=nullp(genes, genome=args$genome, id=args$gene_id, bias.data=gene_lengths, plot.fit=args$make_plots) +if (args$make_plots) { dev.off() } # Fetch GO annotations if category_file hasn't been supplied: -if (category_file == "FALSE") { - go_map=getgo(genes = names(genes), genome=genome, id=gene_id, fetch.cats=fetch_cats) - } else { +if (is.null(args$category_file)) { + go_map=getgo(genes=names(genes), genome=args$genome, id=args$gene_id, fetch.cats=fetch_cats) +} else { # check for header: first entry in first column must be present in genes, else it's a header - first_line = read.delim(category_file, header = FALSE, nrow=1) + first_line = read.delim(args$category_file, header=FALSE, nrow=1) if (first_line[,1] %in% names(genes)) { - go_map = read.delim(category_file, header = FALSE) - } else { - go_map = read.delim(category_file, header= TRUE) - } + go_map = read.delim(args$category_file, header=FALSE) + } else { + go_map = read.delim(args$category_file, header=TRUE) + } } results <- list() -# wallenius approximation of p-values -if (wallenius_tab != FALSE) { - 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, args$wallenius_tab, sep="\t", row.names = FALSE, quote = FALSE) - results[['Wallenius']] <- GO.wall +runGoseq <- function(pwf, genome, gene_id, goseq_method, use_genes_without_cat, repcnt, gene2cat, p_adj_method, out_fp){ + out=goseq(pwf, genome=genome, id=gene_id, method=goseq_method, use_genes_without_cat=use_genes_without_cat, gene2cat=go_map) + out$p.adjust.over_represented = p.adjust(out$over_represented_pvalue, method=p_adj_method) + out$p.adjust.under_represented = p.adjust(out$under_represented_pvalue, method=p_adj_method) + write.table(out, out_fp, sep="\t", row.names=FALSE, quote=FALSE) + return(out) } +# wallenius approximation of p-values +if (!is.null(args$wallenius_tab)) results[['Wallenius']] <- runGoseq( + pwf, + genome=args$genome, + gene_id=args$gene_id, + goseq_method="Wallenius", + use_genes_without_cat=args$use_genes_without_cat, + repcnt=args$repcnt, + gene2cat=go_map, + p_adj_method=args$p_adj_method, + out_fp=args$wallenius_tab) + + # hypergeometric (no length bias correction) -if (nobias_tab != FALSE) { - 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, args$nobias_tab, sep="\t", row.names = FALSE, quote = FALSE) - results[['Hypergeometric']] <- GO.nobias -} +if (!is.null(args$nobias_tab)) results[['Hypergeometric']] <- runGoseq( + pwf, + genome=args$genome, + gene_id=args$gene_id, + goseq_method="Hypergeometric", + use_genes_without_cat=args$use_genes_without_cat, + repcnt=args$repcnt, + gene2cat=go_map, + p_adj_method=args$p_adj_method, + out_fp=args$nobias_tab) # Sampling distribution -if (repcnt > 0) { - - # capture the sampling progress so it doesn't fill stdout - zz <- file("/dev/null", open = "wt") - sink(zz) - GO.samp=goseq(pwf, genome = genome, id = gene_id, method="Sampling", repcnt=repcnt, use_genes_without_cat = use_genes_without_cat, gene2cat=go_map) - sink() +if (args$repcnt > 0){ + results[['Sampling']] <- runGoseq( + pwf, + genome=args$genome, + gene_id=args$gene_id, + goseq_method="Sampling", + use_genes_without_cat=args$use_genes_without_cat, + repcnt=args$repcnt, + gene2cat=go_map, + p_adj_method=args$p_adj_method, + out_fp=args$sampling_tab) - GO.samp$p.adjust.over_represented = p.adjust(GO.samp$over_represented_pvalue, method=p_adj_method) - GO.samp$p.adjust.under_represented = p.adjust(GO.samp$under_represented_pvalue, method=p_adj_method) - write.table(GO.samp, sampling_tab, sep="\t", row.names = FALSE, quote = FALSE) # Compare sampling with wallenius - if (make_plots == TRUE) { - pdf(sample_vs_wallenius_plot) - plot(log10(GO.wall[,2]), log10(GO.samp[match(GO.samp[,1],GO.wall[,1]),2]), - xlab="log10(Wallenius p-values)",ylab="log10(Sampling p-values)", - xlim=c(-3,0)) - abline(0,1,col=3,lty=2) - dev.off() + if (args$make_plots & !is.null(args$wallenius_tab)) { + pdf(args$sample_vs_wallenius_plot) + plot(log10(results[['Wallenius']][,2]), + log10(results[['Sampling']][match(results[['Sampling']][,1], results[['Wallenius']][,1]), 2]), + xlab="log10(Wallenius p-values)", + ylab="log10(Sampling p-values)", + xlim=c(-3,0)) + abline(0,1,col=3,lty=2) + dev.off() } - results[['Sampling']] <- GO.samp } +# Plot the top 10 if (!is.null(args$top_plot)) { cats_title <- gsub("GO:","", args$fetch_cats) # modified from https://bioinformatics-core-shared-training.github.io/cruk-summer-school-2018/RNASeq2018/html/06_Gene_set_testing.nb.html - pdf("top10.pdf") + pdf(args$top_plot) for (m in names(results)) { p <- results[[m]] %>% top_n(10, wt=-over_represented_pvalue) %>% @@ -165,16 +164,35 @@ 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", cats_title), subtitle=paste(m, " method")) + - theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5)) + theme(plot.title=element_text(hjust = 0.5), plot.subtitle=element_text(hjust = 0.5)) print(p) } dev.off() } +# Extract the genes to the categories (GO/KEGG terms) +if (!is.null(args$categories_genes_out_fp)) { + cat2gene = split(rep(names(go_map), sapply(go_map, length)), unlist(go_map, use.names = FALSE)) + # extract categories (GO/KEGG terms) for all results + categories = c() + for (m in names(results)) { + categories = c(categories, results[[m]]$category) + } + categories = unique(categories) + # extract the DE genes for each catge term + categories_genes = data.frame(Categories=categories, DEgenes=rep('', length(categories))) + categories_genes$DEgenes = as.character(categories_genes$DEgenes) + rownames(categories_genes) = categories + for (cat in categories){ + tmp = pwf[cat2gene[[cat]],] + tmp = rownames(tmp[tmp$DEgenes > 0, ]) + categories_genes[cat, 'DEgenes'] = paste(tmp, collapse=',') + } + # output + write.table(categories_genes, args$categories_genes_out_fp, sep = "\t", row.names=FALSE, quote=FALSE) +} + # Output RData file if (!is.null(args$rdata)) { - save.image(file = "goseq_analysis.RData") + save.image(file=args$rdata) } - - -sessionInfo() diff -r 67c29afac85f -r 8b3e3657034e goseq.xml --- a/goseq.xml Sun Mar 17 10:27:17 2019 -0400 +++ b/goseq.xml Fri Sep 06 07:50:46 2019 -0400 @@ -1,14 +1,18 @@ - + tests for overrepresented gene categories + + 1.36.0 + galaxy0 + - bioconductor-goseq - bioconductor-org.hs.eg.db - bioconductor-org.dm.eg.db - bioconductor-org.dr.eg.db - bioconductor-org.mm.eg.db - r-dplyr - r-ggplot2 - r-optparse + bioconductor-goseq + bioconductor-org.hs.eg.db + bioconductor-org.dm.eg.db + bioconductor-org.dr.eg.db + bioconductor-org.mm.eg.db + r-dplyr + r-ggplot2 + r-optparse +#if str($out.make_plots) == 'TRUE' + --make_plots '$out.make_plots' + --length_bias_plot '$length_bias_plot' + #if $methods.repcnt != 0 and $methods.wallenius + --sample_vs_wallenius_plot '$sample_vs_wallenius_plot' + #end if +#end if - +#if $out.cat_genes + --categories_genes_out_fp '$cat_genes_tab' +#end if + +#if $out.rdata_out + --rdata '$rdata' +#end if + ]]> + @@ -96,21 +112,12 @@ -
- - - + + +
- - -
- - - -
-
@@ -121,10 +128,16 @@ - + +
+ +
+ + + +
- methods['wallenius'] @@ -143,33 +156,42 @@ methods['wallenius'] out['make_plots'] - - out['rdata_out'] - - + methods['wallenius'] out['topgo_plot'] + + out['cat_genes'] + + + out['rdata_out'] + - - + - - - - + + + + +
+ + + +
+
+ + +
+
+ + + + +
-
- - - - - - - @@ -177,18 +199,41 @@ - - + + - - - - + + + + + + +
+ + + +
+
+ + +
+
+ + + + +
- + + + + + + +
@@ -196,14 +241,31 @@ - - - - + + + + + + +
+ + + +
+
+ + +
+
+ + + + +
- +
@@ -211,31 +273,76 @@ - - + + + + +
+ + + +
+
+ + +
+
+ + + + +
-
- - - - + + + + +
+ + + +
+
+ + +
+
+ + + + +
- - - - - + + + + +
+ + + +
+
+ + +
+
+ + + + +
@@ -244,14 +351,28 @@ - - - - + + + + +
+ + + +
+
+ + +
+
+ + + + +
-