Mercurial > repos > iuc > goseq

comparison goseq.r @ 9:ef2ad746b589 draft

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 e230a8db9e090c6e0ea9577863ec6153df79e145"
author iuc
date Sun, 06 Jun 2021 22:47:36 +0000
parents 8b3e3657034e
children 602de62d995b
comparison
equal deleted inserted replaced
8:8b3e3657034e 9:ef2ad746b589
1 options( show.error.messages=F, error = function () { cat( geterrmessage(), file=stderr() ); q( "no", 1, F ) } ) 1 options(show.error.messages = F, error = function() {
2 cat(geterrmessage(), file = stderr())
3 q("no", 1, F)
4 })
2 5
3 # we need that to not crash galaxy with an UTF8 error on German LC settings. 6 # 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") 7 loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")
5 8
6 suppressPackageStartupMessages({ 9 suppressPackageStartupMessages({
7 library("goseq") 10 library("goseq")
8 library("optparse") 11 library("optparse")
9 library("dplyr") 12 library("dplyr")
10 library("ggplot2") 13 library("ggplot2")
11 }) 14 })
12 15
13 sessionInfo() 16 sessionInfo()
14 17
15 option_list <- list( 18 option_list <- list(
16 make_option(c("-d", "--dge_file"), type="character", help="Path to file with differential gene expression result"), 19 make_option("--dge_file", type = "character", help = "Path to file with differential gene expression result"),
17 make_option(c("-lf", "--length_file"), type="character", default=NULL, help="Path to tabular file mapping gene id to length"), 20 make_option("--length_file", type = "character", default = NULL, help = "Path to tabular file mapping gene id to length"),
18 make_option(c("-g", "--genome"), type="character", default=NULL, help="Genome [used for looking up correct gene length]"), 21 make_option("--genome", type = "character", default = NULL, help = "Genome [used for looking up correct gene length]"),
19 make_option(c("-i", "--gene_id"), type="character", default=NULL, help="Gene ID format of genes in DGE file"), 22 make_option("--gene_id", type = "character", default = NULL, help = "Gene ID format of genes in DGE file"),
20 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"), 23 make_option("--fetch_cats", type = "character", default = NULL, help = "Categories to get can include one or more of GO:CC, GO:BP, GO:MF, KEGG"),
21 make_option(c("-cat_file", "--category_file"), type="character", default=NULL, help="Path to tabular file with gene_id <-> category mapping"), 24 make_option("--category_file", type = "character", default = NULL, help = "Path to tabular file with gene_id <-> category mapping"),
22 make_option(c("-w","--wallenius_tab"), type="character", default=NULL, help="Path to output file with P-values estimated using wallenius distribution"), 25 make_option("--wallenius_tab", type = "character", default = NULL, help = "Path to output file with P-values estimated using wallenius distribution"),
23 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"), 26 make_option("--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"),
24 make_option(c("-r", "--repcnt"), type="integer", default=0, help="Number of repeats for sampling"), 27 make_option("--repcnt", type = "integer", default = 0, help = "Number of repeats for sampling"),
25 make_option(c("-s","--sampling_tab"), type="character", default=NULL, help="Path to output file with P-values estimated using sampling distribution"), 28 make_option("--sampling_tab", type = "character", default = NULL, help = "Path to output file with P-values estimated using sampling distribution"),
26 make_option(c("-p", "--p_adj_method"), type="character", default="BH", help="Multiple hypothesis testing correction method to use"), 29 make_option("--p_adj_method", type = "character", default = "BH", help = "Multiple hypothesis testing correction method to use"),
27 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)."), 30 make_option("--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)."),
28 make_option(c("-tp", "--top_plot"), type="character", default=NULL, help="Path to output PDF with top10 over-rep GO terms"), 31 make_option("--top_plot", type = "character", default = NULL, help = "Path to output PDF with top10 over-rep GO terms"),
29 make_option(c("-plots", "--make_plots"), default=FALSE, type="logical", help="Produce diagnostic plots?"), 32 make_option("--make_plots", default = FALSE, type = "logical", help = "Produce diagnostic plots?"),
30 make_option(c("-l","--length_bias_plot"), type="character", default=NULL, help="Path to length-bias plot"), 33 make_option("--length_bias_plot", type = "character", default = NULL, help = "Path to length-bias plot"),
31 make_option(c("-sw","--sample_vs_wallenius_plot"), type="character", default=NULL, help="Path to plot comparing sampling with wallenius p-values"), 34 make_option("--sample_vs_wallenius_plot", type = "character", default = NULL, help = "Path to plot comparing sampling with wallenius p-values"),
32 make_option(c("-rd", "--rdata"), type="character", default=NULL, help="Path to RData output file"), 35 make_option("--rdata", type = "character", default = NULL, help = "Path to RData output file"),
33 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") 36 make_option("--categories_genes_out_fp", type = "character", default = NULL, help = "Path to file with categories (GO/KEGG terms) and associated DE genes")
34 ) 37 )
35 38
36 parser <- OptionParser(usage = "%prog [options] file", option_list=option_list) 39 parser <- OptionParser(usage = "%prog [options] file", option_list = option_list)
37 args = parse_args(parser) 40 args <- parse_args(parser)
38 41
39 if (!is.null(args$fetch_cats)) { 42 if (!is.null(args$fetch_cats)) {
40 fetch_cats = unlist(strsplit(args$fetch_cats, ",")) 43 fetch_cats <- unlist(strsplit(args$fetch_cats, ","))
41 } else { 44 } else {
42 fetch_cats = "Custom" 45 fetch_cats <- "Custom"
43 } 46 }
44 47
45 # format DE genes into named vector suitable for goseq 48 # format DE genes into named vector suitable for goseq
46 # check if header is present 49 # check if header is present
47 first_line = read.delim(args$dge_file, header=FALSE, nrow=1) 50 first_line <- read.delim(args$dge_file, header = FALSE, nrow = 1)
48 second_col = toupper(first_line[, ncol(first_line)]) 51 second_col <- toupper(first_line[, ncol(first_line)])
49 if (second_col == TRUE || second_col == FALSE) { 52 if (second_col == TRUE || second_col == FALSE) {
50 dge_table = read.delim(args$dge_file, header=FALSE, sep="\t") 53 dge_table <- read.delim(args$dge_file, header = FALSE, sep = "\t")
51 } else { 54 } else {
52 dge_table = read.delim(args$dge_file, header=TRUE, sep="\t") 55 dge_table <- read.delim(args$dge_file, header = TRUE, sep = "\t")
53 } 56 }
54 genes = as.numeric(as.logical(dge_table[, ncol(dge_table)])) # Last column contains TRUE/FALSE 57 genes <- as.numeric(as.logical(dge_table[, ncol(dge_table)])) # Last column contains TRUE/FALSE
55 names(genes) = dge_table[,1] # Assuming first column contains gene names 58 names(genes) <- dge_table[, 1] # Assuming first column contains gene names
56 59
57 # gene lengths, assuming last column 60 # gene lengths, assuming last column
58 first_line = read.delim(args$length_file, header=FALSE, nrow=1) 61 first_line <- read.delim(args$length_file, header = FALSE, nrow = 1)
59 if (is.numeric(first_line[, ncol(first_line)])) { 62 if (is.numeric(first_line[, ncol(first_line)])) {
60 length_table = read.delim(args$length_file, header=FALSE, sep="\t", check.names=FALSE) 63 length_table <- read.delim(args$length_file, header = FALSE, sep = "\t", check.names = FALSE)
61 } else { 64 } else {
62 length_table = read.delim(args$length_file, header=TRUE, sep="\t", check.names=FALSE) 65 length_table <- read.delim(args$length_file, header = TRUE, sep = "\t", check.names = FALSE)
63 } 66 }
64 row.names(length_table) = length_table[,1] 67 row.names(length_table) <- length_table[, 1]
65 # get vector of gene length in same order as the genes 68 # get vector of gene length in same order as the genes
66 gene_lengths = length_table[names(genes),][, ncol(length_table)] 69 gene_lengths <- length_table[names(genes), ][, ncol(length_table)]
67 70
68 # Estimate PWF 71 # Estimate PWF
69 if (args$make_plots) { 72 if (args$make_plots) {
70 pdf(args$length_bias_plot) 73 pdf(args$length_bias_plot)
71 } 74 }
72 pwf=nullp(genes, genome=args$genome, id=args$gene_id, bias.data=gene_lengths, plot.fit=args$make_plots) 75 pwf <- nullp(genes, genome = args$genome, id = args$gene_id, bias.data = gene_lengths, plot.fit = args$make_plots)
73 if (args$make_plots) { 76 if (args$make_plots) {
74 dev.off() 77 dev.off()
75 } 78 }
76 79
77 # Fetch GO annotations if category_file hasn't been supplied: 80 # Fetch GO annotations if category_file hasn't been supplied:
78 if (is.null(args$category_file)) { 81 if (is.null(args$category_file)) {
79 go_map=getgo(genes=names(genes), genome=args$genome, id=args$gene_id, fetch.cats=fetch_cats) 82 go_map <- getgo(genes = names(genes), genome = args$genome, id = args$gene_id, fetch.cats = fetch_cats)
80 } else { 83 } else {
81 # check for header: first entry in first column must be present in genes, else it's a header 84 # check for header: first entry in first column must be present in genes, else it's a header
82 first_line = read.delim(args$category_file, header=FALSE, nrow=1) 85 first_line <- read.delim(args$category_file, header = FALSE, nrow = 1)
83 if (first_line[,1] %in% names(genes)) { 86 if (first_line[, 1] %in% names(genes)) {
84 go_map = read.delim(args$category_file, header=FALSE) 87 go_map <- read.delim(args$category_file, header = FALSE)
85 } else { 88 } else {
86 go_map = read.delim(args$category_file, header=TRUE) 89 go_map <- read.delim(args$category_file, header = TRUE)
87 } 90 }
88 } 91 }
89 92
90 results <- list() 93 results <- list()
91 94
92 runGoseq <- function(pwf, genome, gene_id, goseq_method, use_genes_without_cat, repcnt, gene2cat, p_adj_method, out_fp){ 95 run_goseq <- function(pwf, genome, gene_id, goseq_method, use_genes_without_cat, repcnt, gene2cat, p_adj_method, out_fp) {
93 out=goseq(pwf, genome=genome, id=gene_id, method=goseq_method, use_genes_without_cat=use_genes_without_cat, gene2cat=go_map) 96 out <- goseq(pwf, genome = genome, id = gene_id, method = goseq_method, use_genes_without_cat = use_genes_without_cat, gene2cat = go_map)
94 out$p.adjust.over_represented = p.adjust(out$over_represented_pvalue, method=p_adj_method) 97 out$p_adjust_over_represented <- p.adjust(out$over_represented_pvalue, method = p_adj_method)
95 out$p.adjust.under_represented = p.adjust(out$under_represented_pvalue, method=p_adj_method) 98 out$p_adjust_under_represented <- p.adjust(out$under_represented_pvalue, method = p_adj_method)
96 write.table(out, out_fp, sep="\t", row.names=FALSE, quote=FALSE) 99 write.table(out, out_fp, sep = "\t", row.names = FALSE, quote = FALSE)
97 return(out) 100 return(out)
98 } 101 }
99 102
100 # wallenius approximation of p-values 103 # wallenius approximation of p-values
101 if (!is.null(args$wallenius_tab)) results[['Wallenius']] <- runGoseq( 104 if (!is.null(args$wallenius_tab)) {
102 pwf, 105 results[["Wallenius"]] <- run_goseq(
103 genome=args$genome, 106 pwf,
104 gene_id=args$gene_id, 107 genome = args$genome,
105 goseq_method="Wallenius", 108 gene_id = args$gene_id,
106 use_genes_without_cat=args$use_genes_without_cat, 109 goseq_method = "Wallenius",
107 repcnt=args$repcnt, 110 use_genes_without_cat = args$use_genes_without_cat,
108 gene2cat=go_map, 111 repcnt = args$repcnt,
109 p_adj_method=args$p_adj_method, 112 gene2cat = go_map,
110 out_fp=args$wallenius_tab) 113 p_adj_method = args$p_adj_method,
114 out_fp = args$wallenius_tab
115 )
116 }
111 117
112 118
113 # hypergeometric (no length bias correction) 119 # hypergeometric (no length bias correction)
114 if (!is.null(args$nobias_tab)) results[['Hypergeometric']] <- runGoseq( 120 if (!is.null(args$nobias_tab)) {
115 pwf, 121 results[["Hypergeometric"]] <- run_goseq(
116 genome=args$genome, 122 pwf,
117 gene_id=args$gene_id, 123 genome = args$genome,
118 goseq_method="Hypergeometric", 124 gene_id = args$gene_id,
119 use_genes_without_cat=args$use_genes_without_cat, 125 goseq_method = "Hypergeometric",
120 repcnt=args$repcnt, 126 use_genes_without_cat = args$use_genes_without_cat,
121 gene2cat=go_map, 127 repcnt = args$repcnt,
122 p_adj_method=args$p_adj_method, 128 gene2cat = go_map,
123 out_fp=args$nobias_tab) 129 p_adj_method = args$p_adj_method,
130 out_fp = args$nobias_tab
131 )
132 }
124 133
125 # Sampling distribution 134 # Sampling distribution
126 if (args$repcnt > 0){ 135 if (args$repcnt > 0) {
127 results[['Sampling']] <- runGoseq( 136 results[["Sampling"]] <- run_goseq(
128 pwf, 137 pwf,
129 genome=args$genome, 138 genome = args$genome,
130 gene_id=args$gene_id, 139 gene_id = args$gene_id,
131 goseq_method="Sampling", 140 goseq_method = "Sampling",
132 use_genes_without_cat=args$use_genes_without_cat, 141 use_genes_without_cat = args$use_genes_without_cat,
133 repcnt=args$repcnt, 142 repcnt = args$repcnt,
134 gene2cat=go_map, 143 gene2cat = go_map,
135 p_adj_method=args$p_adj_method, 144 p_adj_method = args$p_adj_method,
136 out_fp=args$sampling_tab) 145 out_fp = args$sampling_tab
146 )
137 147
138 # Compare sampling with wallenius 148 # Compare sampling with wallenius
139 if (args$make_plots & !is.null(args$wallenius_tab)) { 149 if (args$make_plots & !is.null(args$wallenius_tab)) {
140 pdf(args$sample_vs_wallenius_plot) 150 pdf(args$sample_vs_wallenius_plot)
141 plot(log10(results[['Wallenius']][,2]), 151 plot(log10(results[["Wallenius"]][, 2]),
142 log10(results[['Sampling']][match(results[['Sampling']][,1], results[['Wallenius']][,1]), 2]), 152 log10(results[["Sampling"]][match(results[["Sampling"]][, 1], results[["Wallenius"]][, 1]), 2]),
143 xlab="log10(Wallenius p-values)", 153 xlab = "log10(Wallenius p-values)",
144 ylab="log10(Sampling p-values)", 154 ylab = "log10(Sampling p-values)",
145 xlim=c(-3,0)) 155 xlim = c(-3, 0)
146 abline(0,1,col=3,lty=2) 156 )
157 abline(0, 1, col = 3, lty = 2)
147 dev.off() 158 dev.off()
148 } 159 }
149 } 160 }
150 161
151 # Plot the top 10 162 # Plot the top 10
152 if (!is.null(args$top_plot)) { 163 if (!is.null(args$top_plot)) {
153 cats_title <- gsub("GO:","", args$fetch_cats) 164 cats_title <- gsub("GO:", "", args$fetch_cats)
154 # modified from https://bioinformatics-core-shared-training.github.io/cruk-summer-school-2018/RNASeq2018/html/06_Gene_set_testing.nb.html 165 # modified from https://bioinformatics-core-shared-training.github.io/cruk-summer-school-2018/RNASeq2018/html/06_Gene_set_testing.nb.html
155 pdf(args$top_plot) 166 pdf(args$top_plot)
156 for (m in names(results)) { 167 for (m in names(results)) {
157 p <- results[[m]] %>% 168 p <- results[[m]] %>%
158 top_n(10, wt=-over_represented_pvalue) %>% 169 top_n(10, wt = -over_represented_pvalue) %>%
159 mutate(hitsPerc=numDEInCat*100/numInCat) %>% 170 mutate(hitsPerc = numDEInCat * 100 / numInCat) %>%
160 ggplot(aes(x=hitsPerc, 171 ggplot(aes(
161 y=reorder(substr(term, 1, 40), -over_represented_pvalue), # only use 1st 40 chars of terms otherwise squashes plot 172 x = hitsPerc,
162 colour=p.adjust.over_represented, 173 y = reorder(substr(term, 1, 40), -over_represented_pvalue), # only use 1st 40 chars of terms otherwise squashes plot
163 size=numDEInCat)) + 174 colour = p_adjust_over_represented,
175 size = numDEInCat
176 )) +
164 geom_point() + 177 geom_point() +
165 expand_limits(x=0) + 178 expand_limits(x = 0) +
166 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")) + 179 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")) +
167 theme(plot.title=element_text(hjust = 0.5), plot.subtitle=element_text(hjust = 0.5)) 180 theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5))
168 print(p) 181 print(p)
169 } 182 }
170 dev.off() 183 dev.off()
171 } 184 }
172 185
173 # Extract the genes to the categories (GO/KEGG terms) 186 # Extract the genes to the categories (GO/KEGG terms)
174 if (!is.null(args$categories_genes_out_fp)) { 187 if (!is.null(args$categories_genes_out_fp)) {
175 cat2gene = split(rep(names(go_map), sapply(go_map, length)), unlist(go_map, use.names = FALSE)) 188 cat2gene <- split(rep(names(go_map), sapply(go_map, length)), unlist(go_map, use.names = FALSE))
176 # extract categories (GO/KEGG terms) for all results 189 # extract categories (GO/KEGG terms) for all results
177 categories = c() 190 categories <- c()
178 for (m in names(results)) { 191 for (m in names(results)) {
179 categories = c(categories, results[[m]]$category) 192 categories <- c(categories, results[[m]]$category)
180 } 193 }
181 categories = unique(categories) 194 categories <- unique(categories)
182 # extract the DE genes for each catge term 195 # extract the DE genes for each catge term
183 categories_genes = data.frame(Categories=categories, DEgenes=rep('', length(categories))) 196 categories_genes <- data.frame(category = categories, de_genes = rep("", length(categories)))
184 categories_genes$DEgenes = as.character(categories_genes$DEgenes) 197 categories_genes$de_genes <- as.character(categories_genes$de_genes)
185 rownames(categories_genes) = categories 198 rownames(categories_genes) <- categories
186 for (cat in categories){ 199 for (cat in categories) {
187 tmp = pwf[cat2gene[[cat]],] 200 tmp <- pwf[cat2gene[[cat]], ]
188 tmp = rownames(tmp[tmp$DEgenes > 0, ]) 201 tmp <- rownames(tmp[tmp$DEgenes > 0, ])
189 categories_genes[cat, 'DEgenes'] = paste(tmp, collapse=',') 202 categories_genes[cat, "de_genes"] <- paste(tmp, collapse = ",")
190 } 203 }
191 # output 204 # output
192 write.table(categories_genes, args$categories_genes_out_fp, sep = "\t", row.names=FALSE, quote=FALSE) 205 write.table(categories_genes, args$categories_genes_out_fp, sep = "\t", row.names = FALSE, quote = FALSE)
193 } 206 }
194 207
195 # Output RData file 208 # Output RData file
196 if (!is.null(args$rdata)) { 209 if (!is.null(args$rdata)) {
197 save.image(file=args$rdata) 210 save.image(file = args$rdata)
198 } 211 }