Mercurial > repos > proteore > proteore_topgo

comparison topGO_enrichment.R @ 11:fa2e27165d5d draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit 4efc56eb769fbceb66c64181441ff8781d523454-dirty
author proteore
date Mon, 04 Mar 2019 08:37:49 -0500
parents e3430084c996
children 8eaa43ba1bfc
comparison
equal deleted inserted replaced
10:e3430084c996 11:fa2e27165d5d
153 return(values) 153 return(values)
154 } 154 }
155 155
156 createDotPlot = function(data, onto){ 156 createDotPlot = function(data, onto){
157 157
158 values = deleteInfChar(data$pvalues) 158 values = deleteInfChar(data$pvalues)
159 values = roundValues(values) 159 values = roundValues(values)
160 values = as.numeric(values) 160 values = as.numeric(values)
161 161
162 geneRatio = data$Significant/data$Annotated 162 geneRatio = data$Significant/data$Annotated
163 goTerms = data$Term 163 goTerms = data$Term
164 count = data$Significant 164 count = data$Significant
165 165
166 labely = paste("GO terms",onto,sep=" ") 166 labely = paste("GO terms",onto,sep=" ")
167 ggplot(data,aes(x=geneRatio,y=goTerms, color=values,size=count)) +geom_point( ) + scale_colour_gradientn(colours=c("red","violet","blue")) + xlab("Gene Ratio") + ylab(labely) + labs(color="p-values\n" ) 167 ggplot(data,aes(x=geneRatio,y=goTerms, color=values,size=count)) +geom_point( ) + scale_colour_gradientn(colours=c("red","violet","blue")) + xlab("Gene Ratio") + ylab(labely) + labs(color="p-values\n" )
168 ggsave("dotplot.png", device = "png", dpi = 320, limitsize = TRUE, width = 15, height = 15, units="cm") 168 ggsave("dotplot.png", device = "png", dpi = 320, limitsize = TRUE, width = 15, height = 15, units="cm")
169 } 169 }
170 170
171 createBarPlot = function(data, onto){ 171 createBarPlot = function(data, onto){
172 172
173 173 values = deleteInfChar(data$pvalues)
174 values = deleteInfChar(data$pvalues) 174 values = roundValues(values)
175 values = roundValues(values) 175 values = as.numeric(values)
176 values = as.numeric(values) 176
177 177 goTerms = data$Term
178 goTerms = data$Term 178 count = data$Significant
179 count = data$Significant 179
180 180 labely = paste("GO terms",onto,sep=" ")
181 labely = paste("GO terms",onto,sep=" ") 181 ggplot(data, aes(x=goTerms, y=count,fill=values,scale(scale = 0.5))) + ylab("Gene count") + xlab(labely) +geom_bar(stat="identity") + scale_fill_gradientn(colours=c("red","violet","blue")) + coord_flip() + labs(fill="p-values\n")
182 ggplot(data, aes(x=goTerms, y=count,fill=values,scale(scale = 0.5))) + ylab("Gene count") + xlab(labely) +geom_bar(stat="identity") + scale_fill_gradientn(colours=c("red","violet","blue")) + coord_flip() + labs(fill="p-values\n") 182 ggsave("barplot.png", device = "png", dpi = 320, limitsize = TRUE, width = 15, height = 15, units="cm")
183 ggsave("barplot.png", device = "png", dpi = 320, limitsize = TRUE, width = 15, height = 15, units="cm")
184 } 183 }
185 184
186 # Produce the different outputs 185 # Produce the different outputs
187 createOutputs = function(result, cut_result,text, barplot, dotplot, onto){ 186 createOutputs = function(result, cut_result,text, barplot, dotplot, onto){
188 187
189 if (is.null(result)){ 188 if (is.null(result)){
190 if (text){ 189 err_msg = "None of the input ids can be found in the org package data, enrichment analysis cannot be realized. \n The inputs ids probably either have no associated GO terms or are not ENSG identifiers (e.g : ENSG00000012048)."
191 err_msg = "None of the input ids can be found in the org package data, enrichment analysis cannot be realized. \n The inputs ids probably either have no associated GO terms or are not ENSG identifiers (e.g : ENSG00000012048)." 190 write.table(err_msg, file='result', quote=FALSE, sep='\t', col.names = F, row.names = F)
192 write.table(err_msg, file='result', quote=FALSE, sep='\t', col.names = T, row.names = F) 191 }else if (is.null(cut_result)){
193 } 192 err_msg = "Threshold was too stringent, no GO term found with pvalue equal or lesser than the threshold value."
194 if (barplot){ 193 write.table(err_msg, file='result.tsv', quote=FALSE, sep='\t', col.names = F, row.names = F)
195 png(filename="barplot.png") 194 }else {
196 plot.new() 195 write.table(cut_result, file='result.tsv', quote=FALSE, sep='\t', col.names = T, row.names = F)
197 #text(0,0,err_msg) 196
198 dev.off() 197 if (barplot){createBarPlot(cut_result, onto)}
199 } 198 if (dotplot){createDotPlot(cut_result, onto)}
200 if (dotplot){
201 png(filename="dotplot.png")
202 plot.new()
203 #text(0,0,err_msg)
204 dev.off()
205 }
206 opt <- options(show.error.messages=FALSE)
207 on.exit(options(opt))
208 stop("null result")
209 }
210
211 if (is.null(cut_result)){
212 if (text){
213 err_msg = "Threshold was too stringent, no GO term found with pvalue equal or lesser than the threshold value."
214 write.table(err_msg, file='result', quote=FALSE, sep='\t', col.names = T, row.names = F)
215 }
216 if (barplot){
217 png(filename="barplot.png")
218 plot.new()
219 text(0,0,err_msg)
220 dev.off()
221 }
222 if (dotplot){
223 png(filename="dotplot.png")
224 plot.new()
225 text(0,0,err_msg)
226 dev.off()
227 }
228 opt <- options(show.error.messages=FALSE)
229 on.exit(options(opt))
230 stop("null cut_result")
231 }
232
233 if (text){
234 write.table(cut_result, file='result', quote=FALSE, sep='\t', col.names = T, row.names = F)
235 }
236
237 if (barplot){
238 createBarPlot(cut_result, onto)
239 }
240
241 if (dotplot){
242 createDotPlot(cut_result, onto)
243 } 199 }
244 } 200 }
245 201
246 # Launch enrichment analysis and return result data from the analysis or the null 202 # Launch enrichment analysis and return result data from the analysis or the null
247 # object if the enrichment could not be done. 203 # object if the enrichment could not be done.
316 sample = trimws(unlist(strsplit(tab[,column],";"))) 272 sample = trimws(unlist(strsplit(tab[,column],";")))
317 } 273 }
318 274
319 #check of ENS ids 275 #check of ENS ids
320 if (! any(check_ens_ids(sample))){ 276 if (! any(check_ens_ids(sample))){
321 print("no ensembl gene ids found in your ids list, please check your IDs in input or the selected column of your input file") 277 stop("no ensembl gene ids found in your ids list, please check your IDs in input or the selected column of your input file")
322 stop()
323 } 278 }
324 279
325 #get input if background genes 280 #get input if background genes
326 if (background){ 281 if (background){
327 if (background_input_type=="copy_paste"){ 282 if (background_input_type=="copy_paste"){
330 background_tab=read_file(background_genes,background_header) 285 background_tab=read_file(background_genes,background_header)
331 background_sample = unique(trimws(unlist(strsplit(background_tab[,background_column],";")))) 286 background_sample = unique(trimws(unlist(strsplit(background_tab[,background_column],";"))))
332 } 287 }
333 #check of ENS ids 288 #check of ENS ids
334 if (! any(check_ens_ids(background_sample))){ 289 if (! any(check_ens_ids(background_sample))){
335 print("no ensembl gene ids found in your background ids list, please check your IDs in input or the selected column of your input file") 290 stop("no ensembl gene ids found in your background ids list, please check your IDs in input or the selected column of your input file")
336 stop()
337 } 291 }
338 } else { 292 } else {
339 background_sample=NULL 293 background_sample=NULL
340 } 294 }
341 295