# HG changeset patch
# User proteore
# Date 1545145240 18000
# Node ID d600ce7f2484cdec735dfc90fe0794f5aa84e5fe
# Parent 42d0805353b690b8f66da1785370b9fa47f9dae4
planemo upload commit bdd7e8a1f08c11db2a9f1b6db5535c6d32153b2b-dirty
diff -r 42d0805353b6 -r d600ce7f2484 compute_kegg_pathways.R
--- a/compute_kegg_pathways.R Wed Sep 19 05:38:52 2018 -0400
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,112 +0,0 @@
-library(KEGGREST)
-
-get_args <- function(){
-
- ## Collect arguments
- args <- commandArgs(TRUE)
-
- ## Default setting when no arguments passed
- if(length(args) < 1) {
- args <- c("--help")
- }
-
- ## Help section
- if("--help" %in% args) {
- cat("Pathview R script
- Arguments:
- --help Print this test
- --input tab file
- --id_list
-id list ',' separated
- --id_type type of input ids (uniprot_AC or geneID)
- --id_column number og column containg ids of interest
- --nb_pathways number of pathways to return
- --header boolean
- --output output path
- --ref ref file (l.hsa.gene.RData, l.hsa.up.RData, l.mmu.up.Rdata)
-
- Example:
- Rscript keggrest.R --input='P31946,P62258' --id_type='uniprot' --id_column 'c1' --header TRUE \n\n")
-
- q(save="no")
- }
-
- 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)
-}
-
-args <- get_args()
-
-#save(args,file="/home/dchristiany/proteore_project/ProteoRE/tools/compute_KEGG_pathways/args.Rda")
-#load("/home/dchristiany/proteore_project/ProteoRE/tools/compute_KEGG_pathways/args.Rda")
-
-##function arguments :
-## id.ToMap = input from the user to map on the pathways = list of IDs
-## idType : must be "UNIPROT" or "ENTREZ"
-## org : for the moment can be "Hs" only. Has to evoluate to "Mm"
-
-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)
- }
-}
-
-
-read_file <- function(path,header){
- file <- try(read.table(path,header=header, sep="\t",stringsAsFactors = FALSE, quote=""),silent=TRUE)
- if (inherits(file,"try-error")){
- stop("File not found !")
- }else{
- return(file)
- }
-}
-
-ID2KEGG.Mapping<- function(id.ToMap,ref) {
-
- ref_ids = get(load(ref))
- map<-lapply(ref_ids, is.element, unique(id.ToMap))
- names(map) <- sapply(names(map), function(x) gsub("path:","",x),USE.NAMES = FALSE) #remove the prefix "path:"
-
- in.path<-sapply(map, function(x) length(which(x==TRUE)))
- tot.path<-sapply(map, length)
-
- ratio<-(as.numeric(in.path[which(in.path!=0)])) / (as.numeric(tot.path[which(in.path!=0)]))
- ratio <- as.numeric(format(round(ratio*100, 2), nsmall = 2))
-
- ##useful but LONG
- ## to do before : in step 1
- path.names<-names(in.path[which(in.path!=0)])
- name <- sapply(path.names, function(x) keggGet(x)[[1]]$NAME,USE.NAMES = FALSE)
-
- 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)]))
- res <- res[order(as.numeric(res[,3]),decreasing = TRUE),]
- colnames(res)<-c("pathway_ID", "Description" , "Ratio IDs mapped/total IDs (%)" ,"nb genes mapped in the pathway", "nb total genes present in the pathway")
-
- return(res)
-
-}
-
-###setting variables
-header = str2bool(args$header)
-if (!is.null(args$id_list)) {id_list <- strsplit(args$id_list,",")[[1]]}
-if (!is.null(args$input)) {
- csv <- read_file(args$input,header)
- ncol <- as.numeric(gsub("c", "" ,args$id_column))
- id_list <- as.vector(csv[,ncol])
-}
-id_type <- toupper(args$id_type)
-
-#mapping on pathways
-res <- ID2KEGG.Mapping(id_list,args$ref)
-if (nrow(res) > as.numeric(args$nb_pathways)) { res <- res[1:args$nb_pathways,] }
-
-write.table(res, file=args$output, quote=FALSE, sep='\t',row.names = FALSE, col.names = TRUE)
-
diff -r 42d0805353b6 -r d600ce7f2484 compute_kegg_pathways.xml
--- a/compute_kegg_pathways.xml Wed Sep 19 05:38:52 2018 -0400
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,153 +0,0 @@
-
-
- bioconductor-keggrest
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-@misc{githubKEGGREST,
- title = {KEGGREST: Client-side REST access to KEGG},
- author = {Dan Tenenbaum},
- year = {2018},
- note = {R package version 1.18.1},
- publisher = {GitHub},
- journal = {GitHub repository},
- url = {https://github.com/Bioconductor/KEGGREST},
-}
-
-
diff -r 42d0805353b6 -r d600ce7f2484 entrez_kegg_list.loc.sample
--- a/entrez_kegg_list.loc.sample Wed Sep 19 05:38:52 2018 -0400
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,3 +0,0 @@
-#name date organism value(path)
-Human (Homo sapiens) 27-07-18 hsa tool-data/l.hsa.gene.RData
-Mouse (Mus musculus) 27-07-18 mmu tool-data/l.mmu.gene.RData
diff -r 42d0805353b6 -r d600ce7f2484 kegg_identification.R
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/kegg_identification.R Tue Dec 18 10:00:40 2018 -0500
@@ -0,0 +1,203 @@
+options(warn=-1) #TURN OFF WARNINGS !!!!!!
+
+suppressMessages(library(KEGGREST))
+
+get_args <- function(){
+
+ ## Collect arguments
+ args <- commandArgs(TRUE)
+
+ ## Default setting when no arguments passed
+ if(length(args) < 1) {
+ args <- c("--help")
+ }
+
+ ## Help section
+ if("--help" %in% args) {
+ cat("Pathview R script
+ Arguments:
+ --help Print this test
+ --input tab file
+ --id_list id list ',' separated
+ --id_type type of input ids (kegg-id, uniprot_AC,geneID)
+ --id_column number og column containg ids of interest
+ --nb_pathways number of pathways to return
+ --header boolean
+ --output output path
+ --species species used to get specific pathways (hsa,mmu,rno)
+
+ Example:
+ Rscript keggrest.R --input='P31946,P62258' --id_type='uniprot' --id_column 'c1' --header TRUE \n\n")
+
+ q(save="no")
+ }
+
+ 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)
+}
+
+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)
+ }
+}
+
+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{
+ return(file)
+ }
+}
+
+get_pathways_list <- function(species){
+ ##all available pathways for the species
+ pathways <-keggLink("pathway", species)
+ tot_path<-unique(pathways)
+
+ ##formating the dat into a list object
+ ##key= pathway ID, value = genes of the pathway in the kegg format
+ pathways_list <- sapply(tot_path, function(pathway) names(which(pathways==pathway)))
+ return (pathways_list)
+}
+
+get_list_from_cp <-function(list){
+ list = strsplit(list, "[ \t\n]+")[[1]]
+ list = gsub("[[:blank:]]|\u00A0|NA","",list)
+ list = list[which(!is.na(list[list != ""]))] #remove empty entry
+ list = unique(gsub("-.+", "", list)) #Remove isoform accession number (e.g. "-2")
+ return(list)
+}
+
+geneID_to_kegg <- function(vector,species){
+ vector <- sapply(vector, function(x) paste(species,x,sep=":"),USE.NAMES = F)
+ return (vector)
+}
+
+to_keggID <- function(id_list,id_type){
+ if (id_type == "ncbi-geneid") {
+ id_list <- unique(geneID_to_kegg(id_list,args$species))
+ } else if (id_type=="uniprot"){
+ id_list <- unique(sapply(id_list, function(x) paste(id_type,":",x,sep=""),USE.NAMES = F))
+ if (length(id_list)>250){
+ id_list <- split(id_list, ceiling(seq_along(id_list)/250))
+ id_list <- sapply(id_list, function(x) keggConv("genes",x))
+ id_list <- unique(unlist(id_list))
+ } else {
+ id_list <- unique(keggConv("genes", id_list))
+ }
+ } else if (id_type=="kegg-id") {
+ id_list <- unique(id_list)
+ }
+ return (id_list)
+}
+
+#take data frame, return data frame
+split_ids_per_line <- function(line,ncol){
+
+ #print (line)
+ header = colnames(line)
+ line[ncol] = gsub("[[:blank:]]|\u00A0","",line[ncol])
+
+ if (length(unlist(strsplit(as.character(line[ncol]),";")))>1) {
+ if (length(line)==1 ) {
+ lines = as.data.frame(unlist(strsplit(as.character(line[ncol]),";")),stringsAsFactors = F)
+ } else {
+ if (ncol==1) { #first column
+ lines = suppressWarnings(cbind(unlist(strsplit(as.character(line[ncol]),";")), line[2:length(line)]))
+ } else if (ncol==length(line)) { #last column
+ lines = suppressWarnings(cbind(line[1:ncol-1],unlist(strsplit(as.character(line[ncol]),";"))))
+ } else {
+ lines = suppressWarnings(cbind(line[1:ncol-1], unlist(strsplit(as.character(line[ncol]),";"),use.names = F), line[(ncol+1):length(line)]))
+ }
+ }
+ colnames(lines)=header
+ return(lines)
+ } else {
+ return(line)
+ }
+}
+
+#create new lines if there's more than one id per cell in the columns in order to have only one id per line
+one_id_one_line <-function(tab,ncol){
+
+ if (ncol(tab)>1){
+
+ tab[,ncol] = sapply(tab[,ncol],function(x) gsub("[[:blank:]]","",x))
+ header=colnames(tab)
+ res=as.data.frame(matrix(ncol=ncol(tab),nrow=0))
+ for (i in 1:nrow(tab) ) {
+ lines = split_ids_per_line(tab[i,],ncol)
+ res = rbind(res,lines)
+ }
+ }else {
+ res = unlist(sapply(tab[,1],function(x) strsplit(x,";")),use.names = F)
+ res = data.frame(res[which(!is.na(res[res!=""]))],stringsAsFactors = F)
+ colnames(res)=colnames(tab)
+ }
+ return(res)
+}
+
+kegg_mapping<- function(kegg_id_list,id_type,ref_ids) {
+
+ #mapping
+ map<-lapply(ref_ids, is.element, unique(kegg_id_list))
+ names(map) <- sapply(names(map), function(x) gsub("path:","",x),USE.NAMES = FALSE) #remove the prefix "path:"
+
+ in.path<-sapply(map, function(x) length(which(x==TRUE)))
+ tot.path<-sapply(map, length)
+
+ ratio <- (as.numeric(in.path[which(in.path!=0)])) / (as.numeric(tot.path[which(in.path!=0)]))
+ ratio <- as.numeric(format(round(ratio*100, 2), nsmall = 2))
+
+ ##useful but LONG
+ ## to do before : in step 1
+ path.names<-names(in.path[which(in.path!=0)])
+ name <- sapply(path.names, function(x) keggGet(x)[[1]]$NAME,USE.NAMES = FALSE)
+
+ 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)]))
+ res <- res[order(as.numeric(res[,3]),decreasing = TRUE),]
+ 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")
+
+ return(res)
+
+}
+
+#get args from command line
+args <- get_args()
+
+#save(args,file="/home/dchristiany/proteore_project/ProteoRE/tools/kegg_identification/args.Rda")
+#load("/home/dchristiany/proteore_project/ProteoRE/tools/kegg_identification/args.Rda")
+
+###setting variables
+header = str2bool(args$header)
+if (!is.null(args$id_list)) {id_list <- get_list_from_cp(args$id_list)} #get ids from copy/paste input
+if (!is.null(args$input)) { #get ids from input file
+ csv <- read_file(args$input,header)
+ ncol <- as.numeric(gsub("c", "" ,args$id_column))
+ csv <- one_id_one_line(csv,ncol)
+ id_list <- as.vector(csv[,ncol])
+ id_list <- unique(id_list[which(!is.na(id_list[id_list!=""]))])
+}
+
+#convert to keggID if needed
+id_list <- to_keggID(id_list,args$id_type)
+
+#get pathways of species with associated KEGG ID genes
+pathways_list <- get_pathways_list(args$species)
+
+#mapping on pathways
+res <- kegg_mapping(id_list,args$id_type,pathways_list)
+if (nrow(res) > as.numeric(args$nb_pathways)) { res <- res[1:args$nb_pathways,] }
+
+write.table(res, file=args$output, quote=FALSE, sep='\t',row.names = FALSE, col.names = TRUE)
+
diff -r 42d0805353b6 -r d600ce7f2484 kegg_identification.xml
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/kegg_identification.xml Tue Dec 18 10:00:40 2018 -0500
@@ -0,0 +1,155 @@
+
+ and coverage
+
+ bioconductor-keggrest
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+@misc{githubKEGGREST,
+ title = {KEGGREST: Client-side REST access to KEGG},
+ author = {Dan Tenenbaum},
+ year = {2018},
+ note = {R package version 1.18.1},
+ publisher = {GitHub},
+ journal = {GitHub repository},
+ url = {https://github.com/Bioconductor/KEGGREST},
+}
+
+
diff -r 42d0805353b6 -r d600ce7f2484 test-data/SPZ.soluble.txt
--- a/test-data/SPZ.soluble.txt Wed Sep 19 05:38:52 2018 -0400
+++ b/test-data/SPZ.soluble.txt Tue Dec 18 10:00:40 2018 -0500
@@ -118,16 +118,13 @@
Q8WXX0
P13639
Q14697
-P55809
-A0AVT1
+P55809;A0AVT1
O14980
Q9BVA1
Q14697
O95202
O75694
-Q16851
-P26640
-P23368
+Q16851;P26640;P23368
P55084
P17174
P07814
diff -r 42d0805353b6 -r d600ce7f2484 tool-data/l.hsa.gene.RData
Binary file tool-data/l.hsa.gene.RData has changed
diff -r 42d0805353b6 -r d600ce7f2484 tool-data/l.hsa.up.RData
Binary file tool-data/l.hsa.up.RData has changed
diff -r 42d0805353b6 -r d600ce7f2484 tool-data/l.mmu.gene.RData
Binary file tool-data/l.mmu.gene.RData has changed
diff -r 42d0805353b6 -r d600ce7f2484 tool-data/l.mmu.up.RData
Binary file tool-data/l.mmu.up.RData has changed
diff -r 42d0805353b6 -r d600ce7f2484 tool_data_table_conf.xml.sample
--- a/tool_data_table_conf.xml.sample Wed Sep 19 05:38:52 2018 -0400
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,12 +0,0 @@
-
-
-
- name,date,organism,value
-
-
-
-
- name,date,organism,value
-
-
-
diff -r 42d0805353b6 -r d600ce7f2484 uniprot_kegg_list.loc.sample
--- a/uniprot_kegg_list.loc.sample Wed Sep 19 05:38:52 2018 -0400
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,3 +0,0 @@
-#name date organism value(path)
-Human (Homo sapiens) 27-07-18 hsa tool-data/l.hsa.up.RData
-Mouse (Mus musculus) 27-07-18 mmu tool-data/l.mmu.up.RData