comparison goseq.r @ 4:ae39895af5fe draft

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

3:783e8b70b047

dge_file = args$dge_file
category_file = args$category_file
length_file = args$length_file
genome = args$genome
gene_id = args$gene_id
sampling_tab = args$sampling_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

}

results <- list()

# wallenius approximation of p-values
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, args$wallenius_tab, sep="\t", row.names = FALSE, quote = FALSE)
  results[['Wallenius']] <- GO.wall
}

# hypergeometric (no length bias correction)
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, args$nobias_tab, sep="\t", row.names = FALSE, quote = FALSE)
  results[['Hypergeometric']] <- GO.nobias

  }
  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()
}

4:ae39895af5fe

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

}

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
}

# 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

  }
  results[['Sampling']] <- GO.samp
}

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")
  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=substr(term, 1, 40), # only use 1st 40 chars of terms otherwise squashes plot
                   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", cats_title), subtitle=paste(m, " method")) +
      theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5))
    print(p)
  }
  dev.off()
}