Mercurial > repos > mvdbeek > donotinstall

comparison goseq.r @ 0:8c4a74852e63 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/goseq commit 16ab3a91da29ad4fae86d2c0b1b8b862a3dedcb5-dirty
author mvdbeek
date Thu, 01 Jun 2017 11:23:54 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:8c4a74852e63
1 options( show.error.messages=F, error = function () { cat( geterrmessage(), file=stderr() ); q( "no", 1, F ) } )
2
3 # we need that to not crash galaxy with an UTF8 error on German LC settings.
4 loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")
5
6 suppressPackageStartupMessages({
7 library("goseq")
8 library("optparse")
9 })
10
11 option_list <- list(
12 make_option(c("-d", "--dge_file"), type="character", help="Path to file with differential gene expression result"),
13 make_option(c("-w","--wallenius_tab"), type="character", help="Path to output file with P-values estimated using wallenius distribution."),
14 make_option(c("-s","--sampling_tab"), type="character", default=FALSE, help="Path to output file with P-values estimated using wallenius distribution."),
15 make_option(c("-n","--nobias_tab"), type="character", default=FALSE, help="Path to output file with P-values estimated using wallenius distribution and no correction for gene length bias."),
16 make_option(c("-l","--length_bias_plot"), type="character", default=FALSE, help="Path to length-bias plot."),
17 make_option(c("-sw","--sample_vs_wallenius_plot"), type="character", default=FALSE, help="Path to plot comparing sampling with wallenius p-values."),
18 make_option(c("-r", "--repcnt"), type="integer", default=100, help="Number of repeats for sampling"),
19 make_option(c("-lf", "--length_file"), type="character", default="FALSE", help = "Path to tabular file mapping gene id to length"),
20 make_option(c("-cat_file", "--category_file"), default="FALSE", type="character", help = "Path to tabular file with gene_id <-> category mapping."),
21 make_option(c("-g", "--genome"), default=NULL, type="character", help = "Genome [used for looking up correct gene length]"),
22 make_option(c("-i", "--gene_id"), default=NULL, type="character", help = "Gene ID format of genes in DGE file"),
23 make_option(c("-p", "--p_adj_method"), default="BH", type="character", help="Multiple hypothesis testing correction method to use"),
24 make_option(c("-cat", "--use_genes_without_cat"), default=FALSE, type="logical",
25 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)."),
26 make_option(c("-plots", "--make_plots"), default=FALSE, type="logical", help="produce diagnostic plots?")
27 )
28
29 parser <- OptionParser(usage = "%prog [options] file", option_list=option_list)
30 args = parse_args(parser)
31
32 # Vars:
33 dge_file = args$dge_file
34 category_file = args$category_file
35 length_file = args$length_file
36 genome = args$genome
37 gene_id = args$gene_id
38 wallenius_tab = args$wallenius_tab
39 sampling_tab = args$sampling_tab
40 nobias_tab = args$nobias_tab
41 length_bias_plot = args$length_bias_plot
42 sample_vs_wallenius_plot = args$sample_vs_wallenius_plot
43 repcnt = args$repcnt
44 p_adj_method = args$p_adj_method
45 use_genes_without_cat = args$use_genes_without_cat
46 make_plots = args$make_plots
47
48 # format DE genes into named vector suitable for goseq
49 dge_table = read.delim(dge_file, header = FALSE, sep="\t")
50 genes = as.numeric(as.logical(dge_table[,ncol(dge_table)])) # Last column contains TRUE/FALSE
51 names(genes) = dge_table[,1] # Assuming first column contains gene names
52
53 # gene lengths, assuming last column
54 if (length_file != "FALSE" ) {
55 first_line = read.delim(dge_file, header = FALSE, nrow=1)
56 if (is.numeric(first_line[, ncol(first_line)])) {
57 length_table = read.delim(length_file, header=FALSE, sep="\t", check.names=FALSE)
58 } else {
59 length_table = read.delim(length_file, header=TRUE, sep="\t", check.names=FALSE)
60 }
61 row.names(length_table) = length_table[,1]
62 gene_lengths = length_table[names(genes),][,ncol(length_table)]
63 } else {
64 gene_lengths = getlength(names(genes), genome, gene_id)
65 }
66
67 # Estimate PWF
68
69 if (make_plots == TRUE) {
70 pdf(length_bias_plot)
71 }
72 pwf=nullp(genes, genome = genome, id = gene_id, bias.data = gene_lengths, plot.fit=make_plots)
73 graphics.off()
74
75 # Fetch GO annotations if category_file hasn't been supplied:
76 if (category_file == "FALSE") {
77 go_map=getgo(genes = names(genes), genome = genome, id = gene_id, fetch.cats=c("GO:CC", "GO:BP", "GO:MF", "KEGG"))
78 } else {
79 # check for header: first entry in first column must be present in genes, else it's a header
80 first_line = read.delim(category_file, header = FALSE, nrow=1)
81 if (first_line[,1] %in% names(genes)) {
82 go_map = read.delim(category_file, header = FALSE)
83 } else {
84 go_map = read.delim(category_file, header= TRUE)
85 }
86 }
87
88 # wallenius approximation of p-values
89 if (wallenius_tab != "" && wallenius_tab!="None") {
90 GO.wall=goseq(pwf, genome = genome, id = gene_id, use_genes_without_cat = use_genes_without_cat, gene2cat=go_map)
91 GO.wall$p.adjust.over_represented = p.adjust(GO.wall$over_represented_pvalue, method=p_adj_method)
92 GO.wall$p.adjust.under_represented = p.adjust(GO.wall$under_represented_pvalue, method=p_adj_method)
93 write.table(GO.wall, wallenius_tab, sep="\t", row.names = FALSE, quote = FALSE)
94 }
95
96 # hypergeometric (no length bias correction)
97 if (nobias_tab != "" && nobias_tab != "None") {
98 GO.nobias=goseq(pwf, genome = genome, id = gene_id, method="Hypergeometric", use_genes_without_cat = use_genes_without_cat, gene2cat=go_map)
99 GO.nobias$p.adjust.over_represented = p.adjust(GO.nobias$over_represented_pvalue, method=p_adj_method)
100 GO.nobias$p.adjust.under_represented = p.adjust(GO.nobias$under_represented_pvalue, method=p_adj_method)
101 write.table(GO.nobias, nobias_tab, sep="\t", row.names = FALSE, quote = FALSE)
102 }
103
104 # Sampling distribution
105 if (repcnt > 0) {
106 GO.samp=goseq(pwf, genome = genome, id = gene_id, method="Sampling", repcnt=repcnt, use_genes_without_cat = use_genes_without_cat, gene2cat=go_map)
107 GO.samp$p.adjust.over_represented = p.adjust(GO.samp$over_represented_pvalue, method=p_adj_method)
108 GO.samp$p.adjust.under_represented = p.adjust(GO.samp$under_represented_pvalue, method=p_adj_method)
109 write.table(GO.samp, sampling_tab, sep="\t", row.names = FALSE, quote = FALSE)
110 # Compare sampling with wallenius
111 if (make_plots == TRUE) {
112 pdf(sample_vs_wallenius_plot)
113 plot(log10(GO.wall[,2]), log10(GO.samp[match(GO.samp[,1],GO.wall[,1]),2]),
114 xlab="log10(Wallenius p-values)",ylab="log10(Sampling p-values)",
115 xlim=c(-3,0))
116 abline(0,1,col=3,lty=2)
117 graphics.off()
118 }
119 }
120
121 sessionInfo()