Mercurial > repos > proteore > proteore_topgo

comparison topGO_enrichment.R @ 16:7f1ce70f0f09 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload commit 29a04769c0546be759c38bbcc157c8bc09ec1115-dirty"
author proteore
date Mon, 17 May 2021 14:40:03 +0000
parents 8eaa43ba1bfc
children
comparison
equal deleted inserted replaced
15:159596ec807d 16:7f1ce70f0f09
1 options(warn=-1) #TURN OFF WARNINGS !!!!!! 1 options(warn = -1) #TURN OFF WARNINGS !!!!!!
2 2
3 suppressMessages(library(ggplot2)) 3 suppressMessages(library(ggplot2))
4 suppressMessages(library(topGO)) 4 suppressMessages(library(topGO))
5 5
6 get_args <- function(){ 6 get_args <- function() {
7 7
8 ## Collect arguments 8 ## Collect arguments
9 args <- commandArgs(TRUE) 9 args <- commandArgs(TRUE)
10 10
11 ## Default setting when no arguments passed 11 ## Default setting when no arguments passed
12 if(length(args) < 1) { 12 if (length(args) < 1) {
13 args <- c("--help") 13 args <- c("--help")
14 } 14 }
15 15
16 ## Help section 16 ## Help section
17 if("--help" %in% args) { 17 if ("--help" %in% args) {
18 cat("Pathview R script 18 cat("Pathview R script
19 Arguments: 19 Arguments:
20 --help Print this test 20 --help Print this test
21 --input_type 21 --input_type
22 --onto 22 --onto
23 --option 23 --option
24 --correction 24 --correction
25 --threshold 25 --threshold
26 --text 26 --text
28 --column 28 --column
29 --geneuniverse 29 --geneuniverse
30 --header 30 --header
31 31
32 Example: 32 Example:
33 Rscript --vanilla enrichment_v3.R --inputtype=tabfile (or copypaste) --input=file.txt --ontology='BP/CC/MF' --option=option 33 Rscript --vanilla enrichment_v3.R --inputtype=tabfile (or copypaste)
34 (e.g : classic/elim...) --threshold=threshold --correction=correction --textoutput=text --barplotoutput=barplot --dotplotoutput=dotplot 34 --input=file.txt --ontology='BP/CC/MF' --option=option
35 (e.g : classic/elim...) --threshold=threshold --correction=correction
36 --textoutput=text --barplotoutput=barplot --dotplotoutput=dotplot
35 --column=column --geneuniver=human \n\n") 37 --column=column --geneuniver=human \n\n")
36 38
37 q(save="no") 39 q(save = "no")
38 } 40 }
39 41
40 parseArgs <- function(x) strsplit(sub("^--", "", x), "=") 42 parseargs <- function(x) strsplit(sub("^--", "", x), "=")
41 argsDF <- as.data.frame(do.call("rbind", parseArgs(args))) 43 argsdf <- as.data.frame(do.call("rbind", parseargs(args)))
42 args <- as.list(as.character(argsDF$V2)) 44 args <- as.list(as.character(argsdf$V2))
43 names(args) <- argsDF$V1 45 names(args) <- argsdf$V1
44 46
45 return(args) 47 return(args)
46 } 48 }
47 49
48 read_file <- function(path,header){ 50 read_file <- function(path, header) {
49 file <- try(read.csv(path,header=header, sep="\t",stringsAsFactors = FALSE, quote="\"", check.names = F),silent=TRUE) 51 file <- try(read.csv(path, header = header, sep = "\t",
50 if (inherits(file,"try-error")){ 52 stringsAsFactors = FALSE, quote = "\"", check.names = F),
53 silent = TRUE)
54 if (inherits(file, "try-error")) {
51 stop("File not found !") 55 stop("File not found !")
52 }else{ 56 }else {
53 return(file) 57 return(file)
54 } 58 }
55 } 59 }
56 60
57 get_list_from_cp <-function(list){ 61 get_list_from_cp <- function(list) {
58 list = gsub(";"," ",list) 62 list <- gsub(";", " ", list)
59 list = strsplit(list, "[ \t\n]+")[[1]] 63 list <- strsplit(list, "[ \t\n]+")[[1]]
60 list = list[list != ""] #remove empty entry 64 list <- list[list != ""] #remove empty entry
61 list = gsub("-.+", "", list) #Remove isoform accession number (e.g. "-2") 65 list <- gsub("-.+", "", list) #Remove isoform accession number (e.g. "-2")
62 return(list) 66 return(list)
63 } 67 }
64 68
65 check_ens_ids <- function(vector) { 69 check_ens_ids <- function(vector) {
66 ens_pattern = "^(ENS[A-Z]+[0-9]{11}|[A-Z]{3}[0-9]{3}[A-Za-z](-[A-Za-z])?|CG[0-9]+|[A-Z0-9]+\\.[0-9]+|YM[A-Z][0-9]{3}[a-z][0-9])$" 70 ens_pattern <-
67 return(grepl(ens_pattern,vector)) 71 "^(ENS[A-Z]+[0-9]{11}|[A-Z]{3}[0-9]{3}[A-Za-z](-[A-Za-z])?
68 } 72 |CG[0-9]+|[A-Z0-9]+\\.[0-9]+|YM[A-Z][0-9]{3}[a-z][0-9])$"
69 73 return(grepl(ens_pattern, vector))
70 str2bool <- function(x){ 74 }
71 if (any(is.element(c("t","true"),tolower(x)))){ 75
72 return (TRUE) 76 str2bool <- function(x) {
73 }else if (any(is.element(c("f","false"),tolower(x)))){ 77 if (any(is.element(c("t", "true"), tolower(x)))) {
74 return (FALSE) 78 return(TRUE)
75 }else{ 79 }else if (any(is.element(c("f", "false"), tolower(x)))) {
80 return(FALSE)
81 }else {
76 return(NULL) 82 return(NULL)
77 } 83 }
78 } 84 }
79 85
80 # Some libraries such as GOsummaries won't be able to treat the values such as 86 # Some libraries such as GOsummaries won't be able to treat the values such as
81 # "< 1e-30" produced by topGO. As such it is important to delete the < char 87 # "< 1e-30" produced by topGO. As such it is important to delete the < char
82 # with the deleteInfChar function. Nevertheless the user will have access to the original results in the text output. 88 # with the deleteinfchar function. Nevertheless the user will have access to
83 deleteInfChar = function(values){ 89 #the original results in the text output.
84 90 deleteinfchar <- function(values) {
85 lines = grep("<",values) 91
86 if (length(lines)!=0){ 92 lines <- grep("<", values)
87 for (line in lines){ 93 if (length(lines) != 0) {
88 values[line]=gsub("<","",values[line]) 94 for (line in lines) {
95 values[line] <- gsub("<", "", values[line])
89 } 96 }
90 } 97 }
91 return(values) 98 return(values)
92 } 99 }
93 100
94 corrMultipleTesting = function(result, myGOdata,correction,threshold){ 101 #nolint start
102 corrMultipleTesting = function(result, mygodata, correction, threshold){
95 103
96 # adjust for multiple testing 104 # adjust for multiple testing
97 if (correction!="none"){ 105 if (correction != "none"){
98 # GenTable : transforms the result object into a list. Filters can be applied 106 # GenTable : transforms the result object into a list. Filters can be applied
99 # (e.g : with the topNodes argument, to get for instance only the n first 107 # (e.g : with the topNodes argument, to get for instance only the n first
100 # GO terms with the lowest pvalues), but as we want to apply a correction we 108 # GO terms with the lowest pvalues), but as we want to apply a correction we
101 # take all the GO terms, no matter their pvalues 109 # take all the GO terms, no matter their pvalues
102 allRes <- GenTable(myGOdata, test = result, orderBy = "result", ranksOf = "result",topNodes=length(attributes(result)$score)) 110 allRes <- GenTable(mygodata, test = result, orderBy = "result",
111 ranksOf = "result", topNodes = length(attributes(result)$score))
103 # Some pvalues given by topGO are not numeric (e.g : "<1e-30). As such, these 112 # Some pvalues given by topGO are not numeric (e.g : "<1e-30). As such, these
104 # values are converted to 1e-30 to be able to correct the pvalues 113 # values are converted to 1e-30 to be able to correct the pvalues
105 pvaluestmp = deleteInfChar(allRes$test) 114 pvaluestmp = deleteinfchar(allRes$test)
106 115
107 # the correction is done from the modified pvalues 116 # the correction is done from the modified pvalues
108 allRes$qvalues = p.adjust(pvaluestmp, method = as.character(correction), n = length(pvaluestmp)) 117 allRes$qvalues = p.adjust(pvaluestmp, method = as.character(correction),
118 n = length(pvaluestmp))
109 allRes = as.data.frame(allRes) 119 allRes = as.data.frame(allRes)
110 120
111 # Rename the test column by pvalues, so that is more explicit 121 # Rename the test column by pvalues, so that is more explicit
112 nb = which(names(allRes) %in% c("test")) 122 nb = which(names(allRes) %in% c("test"))
113 123
114 names(allRes)[nb] = "pvalues" 124 names(allRes)[nb] = "pvalues"
115 125
116 allRes = allRes[which(as.numeric(allRes$pvalues) <= threshold),] 126 allRes = allRes[which(as.numeric(allRes$pvalues) <= threshold), ]
117 if (length(allRes$pvalues)==0){ 127 if (length(allRes$pvalues) == 0) {
118 print("Threshold was too stringent, no GO term found with pvalue equal or lesser than the threshold value") 128 print("Threshold was too stringent, no GO term found with pvalue
129 equal or lesser than the threshold value")
119 return(NULL) 130 return(NULL)
120 } 131 }
121 allRes = allRes[order(allRes$qvalues),] 132 allRes = allRes[order(allRes$qvalues), ]
122 } 133 }
123 134
124 if (correction=="none"){ 135 if (correction == "none"){
125 # get all the go terms under user threshold 136 # get all the go terms under user threshold
126 mysummary <- summary(attributes(result)$score <= threshold) 137 mysummary <- summary(attributes(result)$score <= threshold)
127 numsignif <- as.integer(mysummary[[3]]) 138 numsignif <- as.integer(mysummary[[3]])
128 # get all significant nodes 139 # get all significant nodes
129 allRes <- GenTable(myGOdata, test = result, orderBy = "result", ranksOf = "result",topNodes=numsignif) 140 allRes <- GenTable(mygodata, test = result, orderBy = "result",
141 ranksOf = "result", topNodes = numsignif)
130 142
131 143
132 allRes = as.data.frame(allRes) 144 allRes = as.data.frame(allRes)
133 # Rename the test column by pvalues, so that is more explicit 145 # Rename the test column by pvalues, so that is more explicit
134 nb = which(names(allRes) %in% c("test")) 146 nb = which(names(allRes) %in% c("test"))
135 names(allRes)[nb] = "pvalues" 147 names(allRes)[nb] = "pvalues"
136 if (numsignif==0){ 148 if (numsignif == 0) {
137 149
138 print("Threshold was too stringent, no GO term found with pvalue equal or lesser than the threshold value") 150 print("Threshold was too stringent, no GO term found with pvalue
151 equal or lesser than the threshold value")
139 return(NULL) 152 return(NULL)
140 } 153 }
141 154
142 allRes = allRes[order(allRes$pvalues),] 155 allRes = allRes[order(allRes$pvalues), ]
143 } 156 }
144 157
145 return(allRes) 158 return(allRes)
146 } 159 }
147 160 #nolint end
148 # roundValues will simplify the results by rounding down the values. For instance 1.1e-17 becomes 1e-17 161
149 roundValues = function(values){ 162 #roundvalues will simplify the results by rounding down the values.
150 for (line in 1:length(values)){ 163 #For instance 1.1e-17 becomes 1e-17
151 values[line]=as.numeric(gsub(".*e","1e",as.character(values[line]))) 164 roundvalues <- function(values) {
165 for (line in seq_len(length(values))) {
166 values[line] <- as.numeric(gsub(".*e", "1e", as.character(values[line])))
152 } 167 }
153 return(values) 168 return(values)
154 } 169 }
155 170
156 createDotPlot = function(data, onto){ 171 #nolint start
157 172 createDotPlot = function(data, onto) {
158 values = deleteInfChar(data$pvalues) 173
159 values = roundValues(values) 174 values = deleteinfchar(data$pvalues)
175 values = roundvalues(values)
160 values = as.numeric(values) 176 values = as.numeric(values)
161 177
162 geneRatio = data$Significant/data$Annotated 178 geneRatio = data$Significant / data$Annotated
163 goTerms = data$Term 179 goTerms = data$Term
164 count = data$Significant 180 count = data$Significant
165 181
166 labely = paste("GO terms",onto,sep=" ") 182 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" ) 183 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") 184 ggsave("dotplot.png", device = "png", dpi = 320, limitsize = TRUE,
169 } 185 width = 15, height = 15, units = "cm")
170 186 }
171 createBarPlot = function(data, onto){ 187
172 188 createBarPlot = function(data, onto) {
173 values = deleteInfChar(data$pvalues) 189
174 values = roundValues(values) 190 values = deleteinfchar(data$pvalues)
191 values = roundvalues(values)
175 values = as.numeric(values) 192 values = as.numeric(values)
176 193
177 goTerms = data$Term 194 goTerms = data$Term
178 count = data$Significant 195 count = data$Significant
179 196
180 labely = paste("GO terms",onto,sep=" ") 197 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") 198 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") 199 ggsave("barplot.png", device = "png", dpi = 320, limitsize = TRUE,
183 } 200 width = 15, height = 15, units = "cm")
201 }
202
203 #nolint end
184 204
185 # Produce the different outputs 205 # Produce the different outputs
186 createOutputs = function(result, cut_result,text, barplot, dotplot, onto){ 206 createoutputs <- function(result, cut_result, text, barplot, dotplot, onto) {
187 207
188 if (is.null(result)){ 208 if (is.null(result)) {
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)." 209 err_msg <- "None of the input ids can be found in the org package data,
190 write.table(err_msg, file='result', quote=FALSE, sep='\t', col.names = F, row.names = F) 210 enrichment analysis cannot be realized. \n The inputs ids probably
191 }else if (is.null(cut_result)){ 211 either have no associated GO terms or are not ENSG identifiers
192 err_msg = "Threshold was too stringent, no GO term found with pvalue equal or lesser than the threshold value." 212 (e.g : ENSG00000012048)."
193 write.table(err_msg, file='result.tsv', quote=FALSE, sep='\t', col.names = F, row.names = F) 213 write.table(err_msg, file = "result", quote = FALSE, sep = "\t",
214 col.names = F, row.names = F)
215 }else if (is.null(cut_result)) {
216 err_msg <- "Threshold was too stringent, no GO term found with pvalue equal
217 or lesser than the threshold value."
218 write.table(err_msg, file = "result.tsv", quote = FALSE, sep = "\t",
219 col.names = F, row.names = F)
194 }else { 220 }else {
195 write.table(cut_result, file='result.tsv', quote=FALSE, sep='\t', col.names = T, row.names = F) 221 write.table(cut_result, file = "result.tsv", quote = FALSE, sep = "\t",
196 222 col.names = T, row.names = F)
197 if (barplot){createBarPlot(cut_result, onto)} 223
198 if (dotplot){createDotPlot(cut_result, onto)} 224 if (barplot) {
199 } 225 createBarPlot(cut_result, onto) #nolint
200 } 226 }
201 227 if (dotplot) {
202 # Launch enrichment analysis and return result data from the analysis or the null 228 createDotPlot(cut_result, onto) #nolint
203 # object if the enrichment could not be done. 229 }
204 goEnrichment = function(geneuniverse,sample,background_sample,onto){ 230 }
205 231 }
206 if (is.null(background_sample)){ 232
207 xx = annFUN.org(onto,mapping=geneuniverse,ID="ensembl") # get all the GO terms of the corresponding ontology (BP/CC/MF) and all their associated ensembl ids according to the org package 233 # Launch enrichment analysis and return result data from the analysis or the
208 allGenes = unique(unlist(xx)) # check if the genes given by the user can be found in the org package (gene universe), that is in allGenes 234 # null object if the enrichment could not be done.
235 goenrichment <- function(geneuniverse, sample, background_sample, onto) {
236
237 if (is.null(background_sample)) {
238 xx <- annFUN.org(onto, mapping = geneuniverse, ID = "ensembl") #nolint
239 #get all the GO terms of the corresponding ontology (BP/CC/MF)
240 #and all their associated ensembl ids according to the org package
241
242 #nolint start
243
244 allGenes <- unique(unlist(xx))
245 #check if the genes given by the user can be found in the org package
246 #(gene universe), that is in allGenes
209 } else { 247 } else {
210 allGenes = background_sample 248 allGenes <- background_sample
211 } 249 }
212 250
213 if (length(intersect(sample,allGenes))==0){ 251 if (length(intersect(sample,allGenes)) == 0) {
214 print("None of the input ids can be found in the org package data, enrichment analysis cannot be realized. \n The inputs ids probably have no associated GO terms.") 252 print("None of the input ids can be found in the org package data,
215 return(c(NULL,NULL)) 253 enrichment analysis cannot be realized. \n The inputs ids probably
216 } 254 have no associated GO terms.")
217 255 return(c(NULL, NULL))
218 geneList = factor(as.integer(allGenes %in% sample)) #duplicated ids in sample count only for one 256 }
257
258 geneList <- factor(as.integer(allGenes %in% sample))
259 #duplicated ids in sample count only for one
219 if (length(levels(geneList)) == 1 ){ 260 if (length(levels(geneList)) == 1 ){
220 stop("All or none of the background genes are found in tested genes dataset, enrichment analysis can't be done") 261 stop("All or none of the background genes are found in tested genes dataset,
262 enrichment analysis can't be done")
221 } 263 }
222 names(geneList) <- allGenes 264 names(geneList) <- allGenes
223 265
224 #topGO enrichment 266 #nolint end
225 267
268 #topGO enrichment
269
226 # Creation of a topGOdata object 270 # Creation of a topGOdata object
227 # It will contain : the list of genes of interest, the GO annotations and the GO hierarchy 271 # It will contain : the list of genes of interest, the GO annotations and
228 # Parameters : 272 # the GO hierarchy
273 # Parameters :
229 # ontology : character string specifying the ontology of interest (BP, CC, MF) 274 # ontology : character string specifying the ontology of interest (BP, CC, MF)
230 # allGenes : named vector of type numeric or factor 275 # allGenes : named vector of type numeric or factor
231 # annot : tells topGO how to map genes to GO annotations. 276 # annot : tells topGO how to map genes to GO annotations.
232 # argument not used here : nodeSize : at which minimal number of GO annotations 277 # argument not used here : nodeSize : at which minimal number of GO
233 # do we consider a gene 278 # annotations do we consider a gene
234 279
235 myGOdata = new("topGOdata", description="SEA with TopGO", ontology=onto, allGenes=geneList, annot = annFUN.org, mapping=geneuniverse,ID="ensembl") 280 mygodata <- new("topGOdata", description = "SEA with TopGO", ontology = onto,
236 281 allGenes = geneList, annot = annFUN.org,
282 mapping = geneuniverse, ID = "ensembl")
283
237 # Performing enrichment tests 284 # Performing enrichment tests
238 result <- runTest(myGOdata, algorithm=option, statistic="fisher") 285 result <- runTest(mygodata, algorithm = option, statistic = "fisher") #nolint
239 return(c(result,myGOdata)) 286 return(c(result, mygodata))
240 } 287 }
241 288
242 args <- get_args() 289 args <- get_args()
243 290
244 #save(args,file="/home/dchristiany/proteore_project/ProteoRE/tools/topGO/args.rda") 291
245 #load("/home/dchristiany/proteore_project/ProteoRE/tools/topGO/args.rda") 292 input_type <- args$inputtype
246 293 input <- args$input
247 input_type = args$inputtype 294 onto <- args$ontology
248 input = args$input 295 option <- args$option
249 onto = args$ontology 296 correction <- args$correction
250 option = args$option 297 threshold <- as.numeric(args$threshold)
251 correction = args$correction 298 text <- str2bool(args$textoutput)
252 threshold = as.numeric(args$threshold) 299 barplot <- "barplot" %in% unlist(strsplit(args$plot, ","))
253 text = str2bool(args$textoutput) 300 dotplot <- "dotplot" %in% unlist(strsplit(args$plot, ","))
254 barplot = "barplot" %in% unlist(strsplit(args$plot,",")) 301 column <- as.numeric(gsub("c", "", args$column))
255 dotplot = "dotplot" %in% unlist(strsplit(args$plot,",")) 302 geneuniverse <- args$geneuniverse
256 column = as.numeric(gsub("c","",args$column)) 303 header <- str2bool(args$header)
257 geneuniverse = args$geneuniverse 304 background <- str2bool(args$background)
258 header = str2bool(args$header) 305 if (background) {
259 background = str2bool(args$background) 306 background_genes <- args$background_genes
260 if (background){ 307 background_input_type <- args$background_input_type
261 background_genes = args$background_genes 308 background_header <- str2bool(args$background_header)
262 background_input_type = args$background_input_type 309 background_column <- as.numeric(gsub("c", "", args$background_column))
263 background_header = str2bool(args$background_header)
264 background_column = as.numeric(gsub("c","",args$background_column))
265 } 310 }
266 311
267 #get input 312 #get input
268 if (input_type=="copy_paste"){ 313 if (input_type == "copy_paste") {
269 sample <- get_list_from_cp(input) 314 sample <- get_list_from_cp(input)
270 } else if (input_type=="file"){ 315 } else if (input_type == "file") {
271 tab=read_file(input,header) 316 tab <- read_file(input, header)
272 sample = trimws(unlist(strsplit(tab[,column],";"))) 317 sample <- trimws(unlist(strsplit(tab[, column], ";")))
273 } 318 }
274 319
275 #check of ENS ids 320 #check of ENS ids
276 if (! any(check_ens_ids(sample))){ 321 if (! any(check_ens_ids(sample))) {
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("no ensembl gene ids found in your ids list,
323 please check your IDs in input or the selected column of your input file")
278 } 324 }
279 325
280 #get input if background genes 326 #get input if background genes
281 if (background){ 327 if (background) {
282 if (background_input_type=="copy_paste"){ 328 if (background_input_type == "copy_paste") {
283 background_sample <- get_list_from_cp(background_genes) 329 background_sample <- get_list_from_cp(background_genes)
284 } else if (background_input_type=="file"){ 330 } else if (background_input_type == "file") {
285 background_tab=read_file(background_genes,background_header) 331 background_tab <- read_file(background_genes, background_header)
286 background_sample = unique(trimws(unlist(strsplit(background_tab[,background_column],";")))) 332 background_sample <- unique(trimws(unlist(
333 strsplit(background_tab[, background_column], ";"))))
287 } 334 }
288 #check of ENS ids 335 #check of ENS ids
289 if (! any(check_ens_ids(background_sample))){ 336 if (! any(check_ens_ids(background_sample))) {
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") 337 stop("no ensembl gene ids found in your background ids list,
338 please check your IDs in input or the selected column of your input file")
291 } 339 }
292 } else { 340 } else {
293 background_sample=NULL 341 background_sample <- NULL
294 } 342 }
295 343
296 # Launch enrichment analysis 344 # Launch enrichment analysis
297 allresult = suppressMessages(goEnrichment(geneuniverse,sample,background_sample,onto)) 345 allresult <- suppressMessages(goenrichment(geneuniverse, sample,
298 result = allresult[1][[1]] 346 background_sample, onto))
299 myGOdata = allresult[2][[1]] 347 result <- allresult[1][[1]]
300 if (!is.null(result)){ 348 mygodata <- allresult[2][[1]]
301 cut_result = corrMultipleTesting(result,myGOdata, correction,threshold) # Adjust the result with a multiple testing correction or not and with the user, p-value cutoff 349 if (!is.null(result)) {
302 }else{ 350 cut_result <- corrMultipleTesting(result, mygodata, correction, threshold)
303 cut_result=NULL 351 #Adjust the result with a multiple testing correction or not and with the
304 } 352 #user, p-value cutoff
305 353 }else {
306 createOutputs(result, cut_result,text, barplot, dotplot, onto) 354 cut_result <- NULL
355 }
356
357 createoutputs(result, cut_result, text, barplot, dotplot, onto)