Accessing KEGG database from R/Bioconductor

KEGG database is a great resource for biological pathway information, which is an essential part of genome/transcriptome analysis where biological interpretation are formed.

For high throughput studies, it is preferred to access KEGG database programmatically. KEGG recently released the REST service to accomodate such needs. More information about REST service in KEGG can be found from here.

There is also a KEGGREST Bioconductor package available. Since KEGGREST is not working well under my computing environment, I wrote my own code. Below are some examples.

First, a simple example of getting an organism’s basic information.

getOrganismInfo <- function(organism){
  KEGG_INFO_BASE <- "http://rest.kegg.jp/info/"
  info_REST_url <- paste(KEGG_INFO_BASE, organism, sep="")
  info <- readLines(info_REST_url)
  info
}

Calling the above function for human (“hsa” as the organism code), you will get the following:

> getOrganismInfo(“hsa”)
[1] “T01001 Homo sapiens (human) KEGG Genes Database”
[2] “hsa Release 65.0+/02-21, Feb 13″
[3] ” Kanehisa Laboratories”
[4] ” 26,241 entries”

The second example demonstrate how to list KEGG pathway ids and names for a specified organism:

mapPathwayToName <- function(organism) {
  KEGG_PATHWAY_LIST_BASE <- "http://rest.kegg.jp/list/pathway/"
  pathway_list_REST_url <- paste(KEGG_PATHWAY_LIST_BASE, organism, sep="")

  pathway_id_name <- data.frame()

  for (line in readLines(pathway_list_REST_url)) {
    tmp <- strsplit(line, "\t")[[1]]
    pathway_id <- strsplit(tmp[1], organism)[[1]][2]
    pathway_name <- tmp[2]
    pathway_name <- strsplit(pathway_name, "\\s+-\\s+")[[1]][1]
    pathway_id_name[pathway_id, 1] = pathway_name

  }

  names(pathway_id_name) <- "pathway_name"
  pathway_id_name
}

Running this code with “hsa” lists all the human pathways stored in KEGG.

> human_pathways head(human_pathways)
pathway_name
00010 Glycolysis / Gluconeogenesis
00020 Citrate cycle (TCA cycle)
00030 Pentose phosphate pathway
00040 Pentose and glucuronate interconversions
00051 Fructose and mannose metabolism
00052 Galactose metabolism

The third example shows how to map a list of genes to pathways. Note since not all the genes have functions assigned, we’ll only get pathways for some of the genes.

mapGeneToPathway <- function(organism) {
  KEGG_PATHWAY_LINK_BASE <- "http://rest.kegg.jp/link/pathway/"
  pathway_link_REST_url <- paste(KEGG_PATHWAY_LINK_BASE, organism, sep="")
  
  gene_pathway <- data.frame()

  for (line in readLines(pathway_link_REST_url)) {
    tmp <- strsplit(line, "\t")[[1]]
    gene <- tmp[1]
    gene <- strsplit(gene, ":")[[1]][2]  
    pathway_id<- strsplit(tmp[2], organism)[[1]][2]

    if (is.null(gene_pathway[gene, 1])) {
      gene_pathway[gene,1] = pathway_id
    } else {
      if (is.na(gene_pathway[gene,1])) {
        gene_pathway[gene,1] = pathway_id
      } else {
        gene_pathway[gene,1] = paste(gene_pathway[gene, 1], pathway_id, sep=";")
      }
    }
  }
  names(gene_pathway) <- "pathway_id"
  gene_pathway
}

Try it with “hsa”:

> gene_to_pathway <- mapGeneToPathway("hsa")
> head(gene_to_pathway, n=3)
pathway_id
10 00232;00983;01100;05204
100 00230;01100;05340
1000 04514;05412

Note that some gene (10) are mapped to 4 pathway id. Combine this function with the second example, you will be able to add pathway names to the table.