Mercurial > repos > proteore > proteore_kegg_pathways_coverage

comparison kegg_identification.R @ 6:f4e32dee3b28 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload commit 151e7b469b231bbc43c4c39e8e836b05ab6d2253-dirty"
author proteore
date Mon, 17 May 2021 12:29:42 +0000
parents d600ce7f2484
children
comparison
equal deleted inserted replaced
5:86d2f0377e69 6:f4e32dee3b28
1 options(warn=-1) #TURN OFF WARNINGS !!!!!! 1 options(warn = -1) #TURN OFF WARNINGS !!!!!!
2 2
3 suppressMessages(library(KEGGREST)) 3 suppressMessages(library(KEGGREST))
4 4
5 get_args <- function(){ 5 get_args <- function() {
6 6
7 ## Collect arguments 7 ## Collect arguments
8 args <- commandArgs(TRUE) 8 args <- commandArgs(TRUE)
9 9
10 ## Default setting when no arguments passed 10 ## Default setting when no arguments passed
11 if(length(args) < 1) { 11 if (length(args) < 1) {
12 args <- c("--help") 12 args <- c("--help")
13 } 13 }
14 14
15 ## Help section 15 ## Help section
16 if("--help" %in% args) { 16 if ("--help" %in% args) {
17 cat("Pathview R script 17 cat("Pathview R script
18 Arguments: 18 Arguments:
19 --help Print this test 19 --help Print this test
20 --input tab file 20 --input tab file
21 --id_list id list ',' separated 21 --id_list id list ',' separated
22 --id_type type of input ids (kegg-id, uniprot_AC,geneID) 22 --id_type type of input ids (kegg-id, uniprot_AC,geneID)
23 --id_column number og column containg ids of interest 23 --id_column number og column containg ids of interest
24 --nb_pathways number of pathways to return 24 --nb_pathways number of pathways to return
25 --header boolean 25 --header boolean
26 --output output path 26 --output output path
27 --species species used to get specific pathways (hsa,mmu,rno) 27 --species species used to get specific pathways(hsa,mmu,rno)
28 28
29 Example: 29 Example:
30 Rscript keggrest.R --input='P31946,P62258' --id_type='uniprot' --id_column 'c1' --header TRUE \n\n") 30 Rscript keggrest.R --input='P31946,P62258' --id_type='uniprot'
31 31 --id_column 'c1' --header TRUE \n\n")
32 q(save="no") 32
33 } 33 q(save = "no")
34 34 }
35 parseArgs <- function(x) strsplit(sub("^--", "", x), "=") 35
36 argsDF <- as.data.frame(do.call("rbind", parseArgs(args))) 36 parseargs <- function(x) strsplit(sub("^--", "", x), "=")
37 args <- as.list(as.character(argsDF$V2)) 37 argsdf <- as.data.frame(do.call("rbind", parseargs(args)))
38 names(args) <- argsDF$V1 38 args <- as.list(as.character(argsdf$V2))
39 39 names(args) <- argsdf$V1
40
40 return(args) 41 return(args)
41 } 42 }
42 43
43 str2bool <- function(x){ 44 str2bool <- function(x) {
44 if (any(is.element(c("t","true"),tolower(x)))){ 45 if (any(is.element(c("t", "true"), tolower(x)))) {
45 return (TRUE) 46 return(TRUE)
46 }else if (any(is.element(c("f","false"),tolower(x)))){ 47 }else if (any(is.element(c("f", "false"), tolower(x)))) {
47 return (FALSE) 48 return(FALSE)
48 }else{ 49 }else {
49 return(NULL) 50 return(NULL)
50 } 51 }
51 } 52 }
52 53
53 read_file <- function(path,header){ 54 read_file <- function(path, header) {
54 file <- try(read.csv(path,header=header, sep="\t",stringsAsFactors = FALSE, quote="\"", check.names = F),silent=TRUE) 55 file <- try(read.csv(path, header = header, sep = "\t",
55 if (inherits(file,"try-error")){ 56 stringsAsFactors = FALSE, quote = "\"", check.names = F), silent = TRUE)
57 if (inherits(file, "try-error")) {
56 stop("File not found !") 58 stop("File not found !")
57 }else{ 59 }else {
58 return(file) 60 return(file)
59 } 61 }
60 } 62 }
61 63
62 get_pathways_list <- function(species){ 64 get_pathways_list <- function(species) {
63 ##all available pathways for the species 65 ##all available pathways for the species
64 pathways <-keggLink("pathway", species) 66 pathways <- keggLink("pathway", species)
65 tot_path<-unique(pathways) 67 tot_path <- unique(pathways)
66 68
67 ##formating the dat into a list object 69 ##formating the dat into a list object
68 ##key= pathway ID, value = genes of the pathway in the kegg format 70 ##key= pathway ID, value = genes of the pathway in the kegg format
69 pathways_list <- sapply(tot_path, function(pathway) names(which(pathways==pathway))) 71 pathways_list <- sapply(tot_path, function(pathway)
70 return (pathways_list) 72 names(which(pathways == pathway)))
71 } 73 return(pathways_list)
72 74 }
73 get_list_from_cp <-function(list){ 75
74 list = strsplit(list, "[ \t\n]+")[[1]] 76 get_list_from_cp <- function(list) {
75 list = gsub("[[:blank:]]|\u00A0|NA","",list) 77 list <- strsplit(list, "[ \t\n]+")[[1]]
76 list = list[which(!is.na(list[list != ""]))] #remove empty entry 78 list <- gsub("[[:blank:]]|\u00A0|NA", "", list)
77 list = unique(gsub("-.+", "", list)) #Remove isoform accession number (e.g. "-2") 79 list <- list[which(!is.na(list[list != ""]))] #remove empty entry
80 list <- unique(gsub("-.+", "", list))
81 #Remove isoform accession number (e.g. "-2")
78 return(list) 82 return(list)
79 } 83 }
80 84
81 geneID_to_kegg <- function(vector,species){ 85 geneid_to_kegg <- function(vector, species) {
82 vector <- sapply(vector, function(x) paste(species,x,sep=":"),USE.NAMES = F) 86 vector <- sapply(vector, function(x) paste(species, x, sep = ":"),
83 return (vector) 87 USE.NAMES = F)
84 } 88 return(vector)
85 89 }
86 to_keggID <- function(id_list,id_type){ 90
87 if (id_type == "ncbi-geneid") { 91 to_keggid <- function(id_list, id_type) {
88 id_list <- unique(geneID_to_kegg(id_list,args$species)) 92 if (id_type == "ncbi-geneid") {
89 } else if (id_type=="uniprot"){ 93 id_list <- unique(geneid_to_kegg(id_list, args$species))
90 id_list <- unique(sapply(id_list, function(x) paste(id_type,":",x,sep=""),USE.NAMES = F)) 94 }else if (id_type == "uniprot") {
91 if (length(id_list)>250){ 95 id_list <- unique(sapply(id_list, function(x)
92 id_list <- split(id_list, ceiling(seq_along(id_list)/250)) 96 paste(id_type, ":", x, sep = ""), USE.NAMES = F))
93 id_list <- sapply(id_list, function(x) keggConv("genes",x)) 97 if (length(id_list) > 250) {
98 id_list <- split(id_list, ceiling(seq_along(id_list) / 250))
99 id_list <- sapply(id_list, function(x) keggConv("genes", x))
94 id_list <- unique(unlist(id_list)) 100 id_list <- unique(unlist(id_list))
95 } else { 101 } else {
96 id_list <- unique(keggConv("genes", id_list)) 102 id_list <- unique(keggConv("genes", id_list))
97 } 103 }
98 } else if (id_type=="kegg-id") { 104 } else if (id_type == "kegg-id") {
99 id_list <- unique(id_list) 105 id_list <- unique(id_list)
100 } 106 }
101 return (id_list) 107 return(id_list)
102 } 108 }
103 109
104 #take data frame, return data frame 110 #take data frame, return data frame
105 split_ids_per_line <- function(line,ncol){ 111 split_ids_per_line <- function(line, ncol) {
106 112
107 #print (line) 113 #print (line)
108 header = colnames(line) 114 header <- colnames(line)
109 line[ncol] = gsub("[[:blank:]]|\u00A0","",line[ncol]) 115 line[ncol] <- gsub("[[:blank:]]|\u00A0", "", line[ncol])
110 116
111 if (length(unlist(strsplit(as.character(line[ncol]),";")))>1) { 117 if (length(unlist(strsplit(as.character(line[ncol]), ";"))) > 1) {
112 if (length(line)==1 ) { 118 if (length(line) == 1) {
113 lines = as.data.frame(unlist(strsplit(as.character(line[ncol]),";")),stringsAsFactors = F) 119 lines <- as.data.frame(unlist(strsplit(
120 as.character(line[ncol]), ";")), stringsAsFactors = F)
114 } else { 121 } else {
115 if (ncol==1) { #first column 122 if (ncol == 1) { #first column
116 lines = suppressWarnings(cbind(unlist(strsplit(as.character(line[ncol]),";")), line[2:length(line)])) 123 lines <- suppressWarnings(cbind(unlist(strsplit(
117 } else if (ncol==length(line)) { #last column 124 as.character(line[ncol]), ";")), line[2:length(line)]))
118 lines = suppressWarnings(cbind(line[1:ncol-1],unlist(strsplit(as.character(line[ncol]),";")))) 125 } else if (ncol == length(line)) { #last column
126 lines <- suppressWarnings(cbind(line[1:ncol - 1],
127 unlist(strsplit(as.character(line[ncol]), ";"))))
119 } else { 128 } else {
120 lines = suppressWarnings(cbind(line[1:ncol-1], unlist(strsplit(as.character(line[ncol]),";"),use.names = F), line[(ncol+1):length(line)])) 129 lines <- suppressWarnings(cbind(line[1:ncol - 1],
130 unlist(strsplit(as.character(line[ncol]), ";"), use.names = F),
131 line[(ncol + 1):length(line)]))
121 } 132 }
122 } 133 }
123 colnames(lines)=header 134 colnames(lines) <- header
124 return(lines) 135 return(lines)
125 } else { 136 } else {
126 return(line) 137 return(line)
127 } 138 }
128 } 139 }
129 140
130 #create new lines if there's more than one id per cell in the columns in order to have only one id per line 141 #create new lines if there's more than one id per cell in the columns in order
131 one_id_one_line <-function(tab,ncol){ 142 #to have only one id per line
132 143 one_id_one_line <- function(tab, ncol) {
133 if (ncol(tab)>1){ 144
134 145 if (ncol(tab) > 1) {
135 tab[,ncol] = sapply(tab[,ncol],function(x) gsub("[[:blank:]]","",x)) 146
136 header=colnames(tab) 147 tab[, ncol] <- sapply(tab[, ncol], function(x) gsub("[[:blank:]]", "", x))
137 res=as.data.frame(matrix(ncol=ncol(tab),nrow=0)) 148 header <- colnames(tab)
138 for (i in 1:nrow(tab) ) { 149 res <- as.data.frame(matrix(ncol = ncol(tab), nrow = 0))
139 lines = split_ids_per_line(tab[i,],ncol) 150 for (i in seq_len(nrow(tab))) {
140 res = rbind(res,lines) 151 lines <- split_ids_per_line(tab[i, ], ncol)
152 res <- rbind(res, lines)
141 } 153 }
142 }else { 154 } else {
143 res = unlist(sapply(tab[,1],function(x) strsplit(x,";")),use.names = F) 155 res <- unlist(sapply(tab[, 1], function(x) strsplit(x, ";")), use.names = F)
144 res = data.frame(res[which(!is.na(res[res!=""]))],stringsAsFactors = F) 156 res <- data.frame(res[which(!is.na(res[res != ""]))], stringsAsFactors = F)
145 colnames(res)=colnames(tab) 157 colnames(res) <- colnames(tab)
146 } 158 }
147 return(res) 159 return(res)
148 } 160 }
149 161
150 kegg_mapping<- function(kegg_id_list,id_type,ref_ids) { 162 kegg_mapping <- function(kegg_id_list, id_type, ref_ids) {
151 163
152 #mapping 164 #mapping
153 map<-lapply(ref_ids, is.element, unique(kegg_id_list)) 165 map <- lapply(ref_ids, is.element, unique(kegg_id_list))
154 names(map) <- sapply(names(map), function(x) gsub("path:","",x),USE.NAMES = FALSE) #remove the prefix "path:" 166 names(map) <- sapply(names(map), function(x) gsub("path:", "", x),
155 167 USE.NAMES = FALSE) #remove the prefix "path:"
156 in.path<-sapply(map, function(x) length(which(x==TRUE))) 168
157 tot.path<-sapply(map, length) 169 in_path <- sapply(map, function(x) length(which(x == TRUE)))
158 170 tot_path <- sapply(map, length)
159 ratio <- (as.numeric(in.path[which(in.path!=0)])) / (as.numeric(tot.path[which(in.path!=0)])) 171
160 ratio <- as.numeric(format(round(ratio*100, 2), nsmall = 2)) 172 ratio <- (as.numeric(in_path[which(in_path != 0)])) /
161 173 (as.numeric(tot_path[which(in_path != 0)]))
174 ratio <- as.numeric(format(round(ratio * 100, 2), nsmall = 2))
175
162 ##useful but LONG 176 ##useful but LONG
163 ## to do before : in step 1 177 ## to do before : in step 1
164 path.names<-names(in.path[which(in.path!=0)]) 178 path_names <- names(in_path[which(in_path != 0)])
165 name <- sapply(path.names, function(x) keggGet(x)[[1]]$NAME,USE.NAMES = FALSE) 179 name <- sapply(path_names, function(x) keggGet(x)[[1]]$NAME,
166 180 USE.NAMES = FALSE)
167 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)])) 181
168 res <- res[order(as.numeric(res[,3]),decreasing = TRUE),] 182 res <- data.frame(I(names(in_path[which(in_path != 0)])), I(name), ratio,
169 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") 183 as.numeric(in_path[which(in_path != 0)]),
170 184 as.numeric(tot_path[which(in_path != 0)]))
185 res <- res[order(as.numeric(res[, 3]), decreasing = TRUE), ]
186 colnames(res) <- c("pathway_ID", "Description",
187 "Ratio IDs mapped / total IDs (%)",
188 "nb KEGG genes IDs mapped in the pathway",
189 "nb total of KEGG genes IDs present in the pathway")
190
171 return(res) 191 return(res)
172 192
173 } 193 }
174 194
175 #get args from command line 195 #get args from command line
176 args <- get_args() 196 args <- get_args()
177 197
178 #save(args,file="/home/dchristiany/proteore_project/ProteoRE/tools/kegg_identification/args.Rda")
179 #load("/home/dchristiany/proteore_project/ProteoRE/tools/kegg_identification/args.Rda")
180
181 ###setting variables 198 ###setting variables
182 header = str2bool(args$header) 199 header <- str2bool(args$header)
183 if (!is.null(args$id_list)) {id_list <- get_list_from_cp(args$id_list)} #get ids from copy/paste input 200 if (!is.null(args$id_list)) {
184 if (!is.null(args$input)) { #get ids from input file 201 id_list <- get_list_from_cp(args$id_list)
185 csv <- read_file(args$input,header) 202 } #get ids from copy/paste input
186 ncol <- as.numeric(gsub("c", "" ,args$id_column)) 203 if (!is.null(args$input)) { #get ids from input file
187 csv <- one_id_one_line(csv,ncol) 204 csv <- read_file(args$input, header)
188 id_list <- as.vector(csv[,ncol]) 205 ncol <- as.numeric(gsub("c", "", args$id_column))
189 id_list <- unique(id_list[which(!is.na(id_list[id_list!=""]))]) 206 csv <- one_id_one_line(csv, ncol)
207 id_list <- as.vector(csv[, ncol])
208 id_list <- unique(id_list[which(!is.na(id_list[id_list != ""]))])
190 } 209 }
191 210
192 #convert to keggID if needed 211 #convert to keggID if needed
193 id_list <- to_keggID(id_list,args$id_type) 212 id_list <- to_keggid(id_list, args$id_type)
194 213
195 #get pathways of species with associated KEGG ID genes 214 #get pathways of species with associated KEGG ID genes
196 pathways_list <- get_pathways_list(args$species) 215 pathways_list <- get_pathways_list(args$species)
197 216
198 #mapping on pathways 217 #mapping on pathways
199 res <- kegg_mapping(id_list,args$id_type,pathways_list) 218 res <- kegg_mapping(id_list, args$id_type, pathways_list)
200 if (nrow(res) > as.numeric(args$nb_pathways)) { res <- res[1:args$nb_pathways,] } 219 if (nrow(res) > as.numeric(args$nb_pathways)) {
201 220 res <- res[1:args$nb_pathways, ]
202 write.table(res, file=args$output, quote=FALSE, sep='\t',row.names = FALSE, col.names = TRUE) 221 }
203 222
223 write.table(res, file = args$output, quote = FALSE, sep = "\t",
224 row.names = FALSE, col.names = TRUE)