Mercurial > repos > proteore > proteore_clusterprofiler
view GO-enrich.R @ 12:85f039f53414 draft
planemo upload commit b8671ffe2e12dc6612b971a3e6e1dc71496aefd0
author | proteore |
---|---|
date | Fri, 24 Jan 2020 05:12:09 -0500 |
parents | d951677a50d4 |
children | f6107b8ae8f8 |
line wrap: on
line source
options(warn=-1) #TURN OFF WARNINGS !!!!!! suppressMessages(library(clusterProfiler,quietly = TRUE)) # Read file and return file content as data.frame read_file <- function(path,header){ file <- try(read.csv(path,header=header, sep="\t",stringsAsFactors = FALSE, quote="",check.names = F),silent=TRUE) if (inherits(file,"try-error")){ stop("File not found !") }else{ file <- file[!apply(is.na(file) | file == "", 1, all), , drop=FALSE] return(file) } } #return the number of character from the longest description found (from the 10 first) max_str_length_10_first <- function(vector){ vector <- as.vector(vector) nb_description = length(vector) if (nb_description >= 10){nb_description=10} return(max(nchar(vector[1:nb_description]))) } str2bool <- function(x){ if (any(is.element(c("t","true"),tolower(x)))){ return (TRUE) }else if (any(is.element(c("f","false"),tolower(x)))){ return (FALSE) }else{ return(NULL) } } #used before the limit was set to 50 characters width_by_max_char <- function (nb_max_char) { if (nb_max_char < 50 ){ width=600 } else if (nb_max_char < 75) { width=800 } else if (nb_max_char < 100) { width=900 } else { width=1000 } return (width) } repartition_GO <- function(geneid, orgdb, ontology, level=3, readable=TRUE) { ggo<-groupGO(gene=geneid, OrgDb = orgdb, ont=ontology, level=level, readable=TRUE) if (length(ggo@result$ID) > 0 ) { ggo@result$Description <- sapply(as.vector(ggo@result$Description), function(x) {ifelse(nchar(x)>50, substr(x,1,50),x)},USE.NAMES = FALSE) #nb_max_char = max_str_length_10_first(ggo$Description) #width = width_by_max_char(nb_max_char) name <- paste("GGO_", ontology, "_bar-plot", sep = "") png(name,height = 720, width = 600) p <- barplot(ggo, showCategory=10) print(p) dev.off() ggo <- as.data.frame(ggo) return(ggo) } } # GO over-representation test enrich_GO <- function(geneid, universe, orgdb, ontology, pval_cutoff, qval_cutoff,plot) { ego<-enrichGO(gene=geneid, universe=universe, OrgDb=orgdb, ont=ontology, pAdjustMethod="BH", pvalueCutoff=pval_cutoff, qvalueCutoff=qval_cutoff, readable=TRUE) # Plot bar & dot plots #if there are enriched GopTerms if (length(ego$ID)>0){ ego@result$Description <- sapply(ego@result$Description, function(x) {ifelse(nchar(x)>50, substr(x,1,50),x)},USE.NAMES = FALSE) #nb_max_char = max_str_length_10_first(ego$Description) #width = width_by_max_char(nb_max_char) if ("dotplot" %in% plot ){ dot_name <- paste("EGO_", ontology, "_dot-plot", sep = "") png(dot_name,height = 720, width = 600) p <- dotplot(ego, showCategory=10) print(p) dev.off() } if ("barplot" %in% plot ){ bar_name <- paste("EGO_", ontology, "_bar-plot", sep = "") png(bar_name,height = 720, width = 600) p <- barplot(ego, showCategory=10) print(p) dev.off() } ego <- as.data.frame(ego) return(ego) } else { warning(paste("No Go terms enriched (EGO) found for ",ontology,"ontology"),immediate. = TRUE,noBreaks. = TRUE,call. = FALSE) } } clean_ids <- function(ids){ ids = gsub(" ","",ids) ids = ids[which(ids!="")] ids = ids[which(ids!="NA")] ids = ids[!is.na(ids)] return(ids) } check_ids <- function(vector,type) { uniprot_pattern = "^([OPQ][0-9][A-Z0-9]{3}[0-9]|[A-NR-Z][0-9]([A-Z][A-Z0-9]{2}[0-9]){1,2})$" entrez_id = "^([0-9]+|[A-Z]{1,2}_[0-9]+|[A-Z]{1,2}_[A-Z]{1,4}[0-9]+)$" if (type == "entrez") return(grepl(entrez_id,vector)) else if (type == "uniprot") { return(grepl(uniprot_pattern,vector)) } } get_args <- function(){ args <- commandArgs(TRUE) if(length(args)<1) { args <- c("--help") } # Help section if("--help" %in% args) { cat("clusterProfiler Enrichment Analysis Arguments: --input_type: type of input (list of id or filename) --input: input --ncol: the column number which contains list of input IDs --header: true/false if your file contains a header --id_type: the type of input IDs (UniProt/EntrezID) --universe_type: list or filename --universe: background IDs list --uncol: the column number which contains background IDs list --uheader: true/false if the background IDs file contains header --universe_id_type: the type of universe IDs (UniProt/EntrezID) --species --onto_opt: ontology options --go_function: groupGO/enrichGO --level: 1-3 --pval_cutoff --qval_cutoff --text_output: text output filename --plot : type of visualization, dotplot or/and barplot \n") q(save="no") } # Parse arguments parseArgs <- function(x) strsplit(sub("^--", "", x), "=") argsDF <- as.data.frame(do.call("rbind", parseArgs(args))) args <- as.list(as.character(argsDF$V2)) names(args) <- argsDF$V1 return(args) } main <- function() { #get args from command args <- get_args() #save(args,file="/home/dchristiany/proteore_project/ProteoRE/tools/cluster_profiler/args.Rda") #load("/home/dchristiany/proteore_project/ProteoRE/tools/cluster_profiler/args.Rda") go_represent=str2bool(args$go_represent) go_enrich=str2bool(args$go_enrich) if (go_enrich){ plot = unlist(strsplit(args$plot,",")) } suppressMessages(library(args$species, character.only = TRUE, quietly = TRUE)) # Extract OrgDb if (args$species=="org.Hs.eg.db") { orgdb<-org.Hs.eg.db } else if (args$species=="org.Mm.eg.db") { orgdb<-org.Mm.eg.db } else if (args$species=="org.Rn.eg.db") { orgdb<-org.Rn.eg.db } # Extract input IDs input_type = args$input_type id_type = args$id_type if (input_type == "text") { input = unlist(strsplit(strsplit(args$input, "[ \t\n]+")[[1]],";")) } else if (input_type == "file") { filename = args$input ncol = args$ncol # Check ncol if (! as.numeric(gsub("c", "", ncol)) %% 1 == 0) { stop("Please enter the right format for column number: c[number]") } else { ncol = as.numeric(gsub("c", "", ncol)) } header = str2bool(args$header) # Get file content file = read_file(filename, header) # Extract Protein IDs list input = unlist(sapply(as.character(file[,ncol]),function(x) rapply(strsplit(x,";"),c),USE.NAMES = FALSE)) } input = clean_ids(input) ## Get input gene list from input IDs #ID format Conversion #This case : from UNIPROT (protein id) to ENTREZ (gene id) #bitr = conversion function from clusterProfiler if (id_type=="Uniprot" & any(check_ids(input,"uniprot"))) { any(check_ids(input,"uniprot")) idFrom<-"UNIPROT" idTo<-"ENTREZID" suppressMessages(gene<-bitr(input, fromType=idFrom, toType=idTo, OrgDb=orgdb)) gene<-unique(gene$ENTREZID) } else if (id_type=="Entrez" & any(check_ids(input,"entrez"))) { gene<-unique(input) } else { stop(paste(id_type,"not found in your ids list, please check your IDs in input or the selected column of your input file")) } ontology <- strsplit(args$onto_opt, ",")[[1]] ## Extract GGO/EGO arguments if (go_represent) {level <- as.numeric(args$level)} if (go_enrich) { pval_cutoff <- as.numeric(args$pval_cutoff) qval_cutoff <- as.numeric(args$qval_cutoff) # Extract universe background genes (same as input file) if (!is.null(args$universe_type)) { universe_type = args$universe_type if (universe_type == "text") { universe = unlist(strsplit(strsplit(args$input, "[ \t\n]+")[[1]],";")) } else if (universe_type == "file") { universe_filename = args$universe universe_ncol = args$uncol # Check ncol if (! as.numeric(gsub("c", "", universe_ncol)) %% 1 == 0) { stop("Please enter the right format for column number: c[number]") } else { universe_ncol = as.numeric(gsub("c", "", universe_ncol)) } universe_header = str2bool(args$uheader) # Get file content universe_file = read_file(universe_filename, universe_header) # Extract Protein IDs list universe <- unlist(sapply(universe_file[,universe_ncol], function(x) rapply(strsplit(x,";"),c),USE.NAMES = FALSE)) } universe = clean_ids(input) universe_id_type = args$universe_id_type ##to initialize if (universe_id_type=="Uniprot" & any(check_ids(universe,"uniprot"))) { idFrom<-"UNIPROT" idTo<-"ENTREZID" suppressMessages(universe_gene<-bitr(universe, fromType=idFrom, toType=idTo, OrgDb=orgdb)) universe_gene<-unique(universe_gene$ENTREZID) } else if (universe_id_type=="Entrez" & any(check_ids(universe,"entrez"))) { universe_gene<-unique(unlist(universe)) } else { if (universe_type=="text"){ print(paste(universe_id_type,"not found in your background IDs list",sep=" ")) } else { print(paste(universe_id_type,"not found in the column",universe_ncol,"of your background IDs file",sep=" ")) } universe_gene = NULL } } else { universe_gene = NULL } } else { universe_gene = NULL } ##enrichGO : GO over-representation test for (onto in ontology) { if (go_represent) { ggo<-repartition_GO(gene, orgdb, onto, level, readable=TRUE) if (is.list(ggo)){ggo <- as.data.frame(apply(ggo, c(1,2), function(x) gsub("^$|^ $", NA, x)))} #convert "" and " " to NA output_path = paste("cluster_profiler_GGO_",onto,".tsv",sep="") write.table(ggo, output_path, sep="\t", row.names = FALSE, quote = FALSE ) } if (go_enrich) { ego<-enrich_GO(gene, universe_gene, orgdb, onto, pval_cutoff, qval_cutoff,plot) if (is.list(ego)){ego <- as.data.frame(apply(ego, c(1,2), function(x) gsub("^$|^ $", NA, x)))} #convert "" and " " to NA output_path = paste("cluster_profiler_EGO_",onto,".tsv",sep="") write.table(ego, output_path, sep="\t", row.names = FALSE, quote = FALSE ) } } } if(!interactive()) { main() }