Skip to main content eteppo

Import Tidy KEGG Annotations in R

Published: 2023-08-04
Updated: 2023-08-04

Tidy data is so convenient and easy to understand that one hopes to be able to work with it always. Unfortunately sometimes it’s just wasteful or others don’t agree with your mild obsession.

KEGG is an important biological knowledge base and there is already an R package that takes care of fetching data from KEGG into R. Unfortunately this data is not tidy at all. So here’s how you might want to battle with the results to get tidy KEGG metadata for your analyses.

library(tidyverse)
library(assertthat)

# So much can go wrong so you should keep asserting that the data makes sense.
# This is for checking that all rows in a dataframe are distinct.
all_distinct <- function(data, ...) {
  assert_that(is.data.frame(data))
  assert_that(nrow(data) > 0, ncol(data) > 0)
  if (...length() > 0) {
    data <- select(data, ...)
  }
  output <- nrow(data) == nrow(distinct(data))
  return(output)
}

# For checking that there are no missing values in the supplied columns.
all_observed <- function(data, ...) {
	assert_that(is.data.frame(data) | is.matrix(data))
	assert_that(nrow(data) > 0, ncol(data) > 0)
	if (...length() > 0) {
		data <- select(data, ...)
	}
	return(sum(is.na(data)) == 0)
}

# Are all supplied columns of the data unique?
all_unique <- function(data, ...) {
  assert_that(is.data.frame(data))
  assert_that(nrow(data) > 0, ncol(data) > 0)
  if (...length() > 0) {
    data <- select(data, ...)
  }
  output <- all(map_lgl(data, is_unique))
  return(output)
}

# Are all values in the vector unique?
is_unique <- function(x) {
  assert_that(length(x) > 0)
  output <- length(x) == length(unique(remove_na(x)))
  return(output)
}

import_tidy_kegg_rest <- function(databases = "pathway") {

  if (FALSE) {
    # Valid database names for KEGGREST.
    databases <- c(
      "pathway",
      "disease",
      "brite",
      "module",
      "orthology", 
      "genome",
      "genes",
      "compound",
      "glycan",
      "reaction",
      "rclass",
      "enzyme",        
      "network", 
      "variant",
      "drug", 
      "dgroup"
    )
  }

  # Mapping from sets to units.
  kegg_hsa_set_unit_map <- kegg_link_hsa(databases)
  # Mapping between unit IDs. NCBI IDs are commonly used in datasets.
  kegg_hsa_geneid_map <- kegg_conv_hsa(c("NCBI_GeneID" = "ncbi-geneid"))

  # Note: Here I only include human pathways even if organism-specific databases are available.
  kegg_set_names <- kegg_list(databases = databases, organisms = "hsa") %>%
    select(database, set_id = entry_id, set_name = entry_name) %>%
    mutate(set_id = str_c("path:", set_id))

  tidy_kegg <- kegg_hsa_set_unit_map %>%
    left_join(kegg_hsa_geneid_map, by = "unit_id") %>%
    left_join(kegg_set_names, by = c("database", "set_id")) %>%
    # Satisfyingly tidy annotations here.
    select(unit_id, NCBI_GeneID, database, set_id, set_name)
  
  assert_that(all_distinct(tidy_kegg))

  return(tidy_kegg)

}

kegg_link_hsa <- function(databases) {

  assert_that(is.character(databases))
  assert_that(length(databases) > 0)

  get_link <- function(database) {
    # KEGGREST-package has functions for different "mappings" we need.
    mapping_chr <- KEGGREST::keggLink(database, "hsa")
    mapping_df <- tibble::tibble(
        set_id = unname(mapping_chr),
        unit_id = names(mapping_chr)
      ) %>%
      mutate(database = database)
    return(mapping_df)
  }

  mappings <- databases %>%
    map_dfr(get_link)

  assert_that(all_distinct(mappings))
  assert_that(all_observed(mappings))

  return(mappings)

}

kegg_conv_hsa <- function(target_id) {
  assert_that(is_string(target_id))
  assert_that(is_string(names(target_id)))

  response <- KEGGREST::keggConv("hsa", target_id) 
  target_ids_chr <- response %>%
    names() %>%
    str_remove(target_id) %>%
    str_remove("^:")
  source_ids_chr <- unname(response)
  id_mapping <- tibble::tibble(
      target_id = target_ids_chr,
      source_id = source_ids_chr
    ) %>%
    magrittr::set_colnames(c(names(target_id), "unit_id"))

  assert_that(all_distinct(id_mapping))
  assert_that(all_unique(id_mapping))

  return(id_mapping)
}

kegg_list <- function(databases, organisms = NULL) {

  assert_that(is.character(databases))
  assert_that(length(databases) > 0)
  assert_that(is.null(organisms) | is.character(organisms))

  get_list <- function(database, organisms) {
    entries_chr <- KEGGREST::keggList(database)
    if (!is.null(organisms) & database == "pathway") {
      entries_org_chr <- c()
      for (organism in organisms) {
        entries_chr <- KEGGREST::keggList(database, organism = organism)
        entries_org_chr <- append(entries_org_chr, entries_chr)
      }
      entries_chr <- append(entries_chr, entries_org_chr)
    }
    entries_df <- tibble::tibble(
        entry_id = names(entries_chr),
        entry_name = unname(entries_chr)
      ) %>%
      mutate(database = database)
    return(entries_df)
  }

  entry_name_data <- databases %>%
    map_dfr(get_list, organisms = organisms) %>%
    distinct(entry_id, entry_name, .keep_all = TRUE)

  # All mappings should be unique but they might not be.
  if (not(all_unique(entry_name_data, entry_id:entry_name))) {
    n_onetomany <- entry_name_data %>%
        select(entry_id, entry_name) %>%
        group_by(entry_id) %>%
        summarize(n = n()) %>%
        filter(n > 1) %>%
        nrow()
    n_manytoone <- entry_name_data %>%
        select(entry_id, entry_name) %>%
        group_by(entry_name) %>%
        summarize(n = n()) %>%
        filter(n > 1) %>%
        nrow()
    n_missing <- nrow(non_complete(entry_name_data))
    warning_msg <- str_c(
        glue("Metadata mapping KEGG IDs to names had {n_onetomany} one-to-many,"),
        glue("{n_manytoone} many-to-one, and {n_missing} missing entries."),
        "Only the first name of each ID was selected.",
        sep = " "
    )
    warning(warning_msg)

    entry_name_data <- entry_name_data %>%
        distinct(entry_id, .keep_all = TRUE)
  }

  assert_that(all_unique(entry_name_data, entry_id))
  return(entry_name_data)

}

# You can manage your files in a common named list.
# I name the path to the resulting file here "kegg_annotations_file".
import_tidy_kegg <- function(output_files) {

  kegg_annotations_file <- chuck(output_files, "kegg_annotations_file")

  # If the file already exists, it will just read it. It works like a cache.
  if (file.exists(kegg_annotations_file)) {
    tidy_kegg <- kegg_annotations_file %>%
      read_tsv(col_types = cols(.default = "c"))

    return(tidy_kegg)
  }

  tidy_kegg <- import_tidy_kegg_rest()

  assert_that(all_distinct(tidy_kegg))
  write_tsv(tidy_kegg, file = kegg_annotations_file)

  return(tidy_kegg)

}