Mercurial > repos > proteore > proteore_kegg_pathways_coverage

comparison compute_kegg_pathways.R @ 0:8883a7173cba draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit 63302cb49d4f0f4dbc9ae141d20704822588f54e-dirty
author proteore
date Mon, 12 Nov 2018 10:59:49 -0500
parents
children 7004924a3686
comparison
equal deleted inserted replaced
-1:000000000000 0:8883a7173cba
1 options(warn=-1) #TURN OFF WARNINGS !!!!!!
2
3 suppressMessages(library(KEGGREST))
4
5 get_args <- function(){
6
7 ## Collect arguments
8 args <- commandArgs(TRUE)
9
10 ## Default setting when no arguments passed
11 if(length(args) < 1) {
12 args <- c("--help")
13 }
14
15 ## Help section
16 if("--help" %in% args) {
17 cat("Pathview R script
18 Arguments:
19 --help Print this test
20 --input tab file
21 --id_list id list ',' separated
22 --id_type type of input ids (kegg-id, uniprot_AC,geneID)
23 --id_column number og column containg ids of interest
24 --nb_pathways number of pathways to return
25 --header boolean
26 --output output path
27 --species species used to get specific pathways (hsa,mmu,rno)
28
29 Example:
30 Rscript keggrest.R --input='P31946,P62258' --id_type='uniprot' --id_column 'c1' --header TRUE \n\n")
31
32 q(save="no")
33 }
34
35 parseArgs <- function(x) strsplit(sub("^--", "", x), "=")
36 argsDF <- as.data.frame(do.call("rbind", parseArgs(args)))
37 args <- as.list(as.character(argsDF$V2))
38 names(args) <- argsDF$V1
39
40 return(args)
41 }
42
43 str2bool <- function(x){
44 if (any(is.element(c("t","true"),tolower(x)))){
45 return (TRUE)
46 }else if (any(is.element(c("f","false"),tolower(x)))){
47 return (FALSE)
48 }else{
49 return(NULL)
50 }
51 }
52
53 read_file <- function(path,header){
54 file <- try(read.csv(path,header=header, sep="\t",stringsAsFactors = FALSE, quote="\"", check.names = F),silent=TRUE)
55 if (inherits(file,"try-error")){
56 stop("File not found !")
57 }else{
58 return(file)
59 }
60 }
61
62 get_pathways_list <- function(species){
63 ##all available pathways for the species
64 pathways <-keggLink("pathway", species)
65 tot_path<-unique(pathways)
66
67 ##formating the dat into a list object
68 ##key= pathway ID, value = genes of the pathway in the kegg format
69 pathways_list <- sapply(tot_path, function(pathway) names(which(pathways==pathway)))
70 return (pathways_list)
71 }
72
73 get_list_from_cp <-function(list){
74 list = strsplit(list, "[ \t\n]+")[[1]]
75 list = list[list != ""] #remove empty entry
76 list = gsub("-.+", "", list) #Remove isoform accession number (e.g. "-2")
77 return(list)
78 }
79
80 geneID_to_kegg <- function(vector,species){
81 vector <- sapply(vector, function(x) paste(species,x,sep=":"),USE.NAMES = F)
82 return (vector)
83 }
84
85 kegg_mapping<- function(kegg_id_list,id_type,ref_ids) {
86
87 #convert to KEGG ID
88 #if (id_type!="kegg-id"){
89 # id_list <- unique(sapply(id_list, function(x) paste(id_type,":",x,sep=""),USE.NAMES = F))
90 # if (length(id_list)>250){
91 # id_list <- split(id_list, ceiling(seq_along(id_list)/250))
92 # id_list <- sapply(id_list, function(x) keggConv("genes",x))
93 # kegg_id_list <- unique(unlist(id_list))
94 # } else {
95 # kegg_id_list <- unique(keggConv("genes", id_list))
96 # }
97 #} else {
98 # kegg_id_list <- unique(id_list)
99 #}
100
101 #mapping
102 map<-lapply(ref_ids, is.element, unique(kegg_id_list))
103 names(map) <- sapply(names(map), function(x) gsub("path:","",x),USE.NAMES = FALSE) #remove the prefix "path:"
104
105 in.path<-sapply(map, function(x) length(which(x==TRUE)))
106 tot.path<-sapply(map, length)
107
108 ratio <- (as.numeric(in.path[which(in.path!=0)])) / (as.numeric(tot.path[which(in.path!=0)]))
109 ratio <- as.numeric(format(round(ratio*100, 2), nsmall = 2))
110
111 ##useful but LONG
112 ## to do before : in step 1
113 path.names<-names(in.path[which(in.path!=0)])
114 name <- sapply(path.names, function(x) keggGet(x)[[1]]$NAME,USE.NAMES = FALSE)
115
116 res<-data.frame(I(names(in.path[which(in.path!=0)])), I(name), ratio, as.numeric(in.path[which(in.path!=0)]), as.numeric(tot.path[which(in.path!=0)]))
117 res <- res[order(as.numeric(res[,3]),decreasing = TRUE),]
118 colnames(res)<-c("pathway_ID", "Description" , "Ratio IDs mapped/total IDs (%)" ,"nb KEGG genes IDs mapped in the pathway", "nb total of KEGG genes IDs present in the pathway")
119
120 return(res)
121
122 }
123
124 #get args from command line
125 args <- get_args()
126
127 #save(args,file="/home/dchristiany/proteore_project/ProteoRE/tools/kegg_pathways_identification/args.Rda")
128 #load("/home/dchristiany/proteore_project/ProteoRE/tools/kegg_pathways_identification/args.Rda")
129
130 ###setting variables
131 header = str2bool(args$header)
132 if (!is.null(args$id_list)) {id_list <- get_list_from_cp(args$id_list)}
133 if (!is.null(args$input)) {
134 csv <- read_file(args$input,header)
135 ncol <- as.numeric(gsub("c", "" ,args$id_column))
136 id_list <- as.vector(csv[,ncol])
137 id_list <- id_list[which(!is.na(id_list))]
138 }
139 if (args$id_type == "ncbi-geneid") {
140 id_list <- geneID_to_kegg(id_list,args$species)
141 }
142
143
144 #get pathways of species with associated KEGG ID genes
145 pathways_list <- get_pathways_list(args$species)
146
147 #mapping on pathways
148 res <- kegg_mapping(id_list,args$id_type,pathways_list)
149 if (nrow(res) > as.numeric(args$nb_pathways)) { res <- res[1:args$nb_pathways,] }
150
151 write.table(res, file=args$output, quote=FALSE, sep='\t',row.names = FALSE, col.names = TRUE)
152