2.1 生物種コードを検索.

• wget -q -O - http://rest.kegg.jp/list/organism
  • 戻り値をそのまま使いたいので, -qでダウンロードの詳細情報を入れないようにする.
  • -O -内容を標準出力に出す.
• system関数の戻り値としてキャプチャして.ベクトルに入れる.

# ダウンロード
org <- system("wget -q -O - http://rest.kegg.jp/list/organism", intern = T)

# データフレームに変換
org_dat <- data.frame(T_number  = sapply(strsplit(org, "\t"), "[", 1),
                  organism_code  = sapply(strsplit(org, "\t"), "[", 2),
                  species = sapply(strsplit(org, "\t"), "[", 3),
                  stringsAsFactors = F
                  )

2.2 MAP IDと名前

• wget -q -O - http://rest.kegg.jp/list/pathway/cge # cgeは生物種コード
• wget -q -O - http://rest.kegg.jp/list/pathway # 全生物種の場合
  • 全PATHWAYを取得する場合、生物種を指定しない.

# ダウンロード
pathway <- system("wget -q -O - http://rest.kegg.jp/list/pathway/cge", intern = T)　

# データフレームに変換
map_dat <- data.frame(map = sapply(strsplit(pathway, "\t"), "[", 1),
                      name  = sapply(strsplit(pathway, "\t"), "[", 2),
                      stringsAsFactors = F) 

3.1 IDを取得

• KEGG GENEのidとdescriptionが取得できる.
• wget -q -O - http://rest.kegg.jp/list/cge

# ダウンロード
gene <- system("wget -q -O - http://rest.kegg.jp/list/cge", intern = T)

# データフレームに変換
gi_cge <- 
  data.frame(gene = sapply(strsplit(gene, "\t"), "[", 1),
             description = sapply(strsplit(gene, "\t"), "[", 2),
             stringsAsFactors = F
             )

3.2 FASTA配列取得

• wget -q -O - http://rest.kegg.jp/get/cge:100765387/ntseq # 塩基配列
• wget -q -O - http://rest.kegg.jp/get/cge:100765387/aaseq # アミノ酸配列
• idの部分を + で繋げてcge:100682525+cge:100682526+cge:100682527のようにしてもよい.
  • ただし10個以上の配列をつなげても返ってこない. idを10個づつ+で繋げて実行する.
• 特定の生物群をまとめて行う場合, 生物種コードのリストをhttps://www.genome.jp/kegg/catalog/org_list.html から取得する.

# idを10づつに分割してアミノ酸配列を取得する. 
gi <- gi_cge$gene[1:29]

# id chunk per 10 seq 
chunk <- function(x, n){
  nchnk <- ceiling(length(x)/n)
  split(x, cut(seq_along(x), nchnk, labels = FALSE))
}
chnk <- chunk(x = gi, n = 10)

# 10 seqづつwget 
out_fna <- "~/pub/fasta/cge_29.fna" 
fna <- file(out_fna, "w") # ファイルコネクションを書き込みモードで開く

# fasta配列を取得し, wireLinesで順次書き出し
for(i in seq_along(chnk)){
  gids_chnk <- paste(chnk[[i]], collapse = "+")
  com <- paste0("wget -q -O - http://rest.kegg.jp/get/", gids_chnk, "/ntseq")
  fst <- system(com, wait = T, intern = T)
  writeLines(fst, fna)
}
close(fna) # ファイルコネクション閉じる

3.3 KEGG GENEのIDと KEGG PATHWAYのIDとの対応

• wget -q -O - http://rest.kegg.jp/link/pathway/cge

# ダウンロード
genemap <- system("wget -q -O - http://rest.kegg.jp/link/pathway/cge", intern = T) 

# データフレームに変換
map_cge <- 
  data.frame(gene = sapply(strsplit(genemap, "\t"), "[", 1),
             map = sub("path:","",sapply(strsplit(genemap, "\t"), "[", 2)),
             stringsAsFactors = F
             )

# ユニークなGENE IDにして　MAP IDを";"で結合
maps_cge <- map_cge %>% group_by(gene) %>% summarise(maps = paste(map, collapse = ";"))

3.4 KEGG GENEとenzyme code

• wget -q -O - http://rest.kegg.jp/link/ec/cge

# ダウンロード
ec_cge <- system("wget -q -O - http://rest.kegg.jp/link/ec/cge", intern = T)

# データフレームに変換
ec_cge <- data.frame(gene = sapply(strsplit(ec_cge, "\t"), "[", 1),
                     ec = sapply(strsplit(ec_cge, "\t"), "[", 2), 
                     stringsAsFactors = F
                     )

# 複数のenzyme codeが付いているものは";"で結合する. 
ecs_cge <- ec_cge %>% group_by(gene) %>% summarise(ecs = paste(ec, collapse = "; "))

3.5 KEGG GENEとRefseq_proteinの対応

• wget -q -O - http://rest.kegg.jp/conv/ncbi-proteinid/cge

#　ダウンロード
refp <- system("wget -q -O - http://rest.kegg.jp/conv/ncbi-proteinid/cge", intern = T)

# データフレームに変換
refp_cge <- data.frame(gene = sapply(strsplit(refp, "\t"), "[", 1), 
                       ref.protein = sub("ncbi-proteinid:","", 
                                         sapply(strsplit(refp, "\t"), "[", 2)), 
                       stringsAsFactors = F)

# length(unique(refp_cge$gene)) == nrow(refp_cge)

3.6 KEGG GENEとUniprotの対応

• wget -q -O - http://rest.kegg.jp/conv/uniprot/cge
• 真核生物のKEGG GENEはRefSeqを元にしているデータが多い.
• RefSeq Proteinとは1対1対応だが、UniProtのほうは複数のIDが対応する場合がある.

# ダウンロード
unip <- system("wget -q -O - http://rest.kegg.jp/conv/uniprot/cge", intern = T)

# データフレームに変換
unip_cge <- data.frame(gene = sapply(strsplit(unip, "\t"), "[", 1),
             upid = sapply(strsplit(unip, "\t"), "[", 2),
             stringsAsFactors = F)

# 複数のUniProt IDが付いているものは";"で結合する.
unips_cge <- unip_cge %>% group_by(gene) %>% summarise(upids = paste(upid, collapse = "; "))

3.7 KEGG GENEとncbi-geneとの対応

# ダウンロード
ncbi <- system("wget -q -O - http://rest.kegg.jp/conv/ncbi-geneid/cge", intern = T)

# データフレームに変換
ncbi_cge <- data.frame(gene = sapply(strsplit(ncbi, "\t"), "[", 1),
             ncbi_gi = sub("ncbi-geneid:", "", sapply(strsplit(ncbi, "\t"), "[", 2)),
             stringsAsFactors = F)

3.8 KEGG KO (ORTHOLOGY)とKEGG GENEとの対応

• wget -q -O - http://rest.kegg.jp/link/ko/cge
• 1つのGENEのIDに複数のKOが付いている場合がある.

# ダウンロード
ko <- system("wget -q -O - http://rest.kegg.jp/link/ko/cge", intern = T)

# データフレームに変換
ko_cge <- data.frame(gene = sapply(strsplit(ko, "\t"), "[", 1),
             ko = sub("ko:","", sapply(strsplit(ko, "\t"), "[", 2)), 
             stringsAsFactors = F)

# 複数のKO IDが付いているものは";"で結合する.
kos_cge <- ko_cge %>% group_by(gene) %>% summarise(kos = paste(ko, collapse = ";"))

# kos_cge$kos[duplicated(kos_cge$kos)]
# kos_cge$gene[duplicated(kos_cge$gene)]
# ko_cge[ko_cge$gene == "cge:100770162",] # K00504, K18200

6.1 代謝マップ上に乗っている化合物.

# list/cpd
cpd <- system("wget -q -O - http://rest.kegg.jp/list/cpd", intern = T, wait = T)
cpdat <- data.frame(cid = sub("cpd:", "", sapply(strsplit(cpd, "\t"), "[", 1)),
                    description = sapply(strsplit(cpd, "\t"), "[", 2), 
                    stringsAsFactors = F)

# link/pathway/cpd
cpdmap <- system("wget -q -O - http://rest.kegg.jp/link/pathway/cpd", intern = T, wait = T)
cpd_map <- data.frame(cid = sub("cpd:", "", sapply(strsplit(cpdmap, "\t"), "[", 1)),
                    map = sub("path:", "", sapply(strsplit(cpdmap, "\t"), "[", 2)), 
                    stringsAsFactors = F)


# 代謝マップ上に載っている化合物
cpd_maps <- cpd_map %>%
  group_by(cid) %>%
  summarise(maps = paste(map, collapse = ";")) %>%
  left_join(., cpdat, by = "cid")

6.2 PubChem Substanceとのid対応

• KEGG COMPOUNDのPubChemへのリンクはPubChem Compoundではなく, PubChem Substanceのidであることに注意.
  • wget -q -O - http://rest.kegg.jp/conv/pubchem/cpd
• PubChemのftpサイトからデータDL
  • PubChem Substance: エントリー毎にID(SID)が振られるため重複がある. ftp
    
  • PubChem Compound: PubChemのルールに基づいて, 一意なID(CID) ftp
    
  • PubChem Identifier Exchange Service: IDの変換
    
  • Web APIを使ってデータ取得する方法もある.

# KEGG COMPOUND ID と PubChem SID
cpdpub <- system("wget -q -O - http://rest.kegg.jp/conv/pubchem/cpd", intern = T, wait = T)
cpd_pub <- data.frame(
  cid = sub("cpd:", "", sapply(strsplit(cpdpub, "\t"), "[", 1)),
  pubchem_SID = sub("pubchem:", "", sapply(strsplit(cpdpub, "\t"), "[", 2)), 
  stringsAsFactors = F
  )

# PubChem SIDとPubChem CID (PubChem Identifier Exchange Service からDL)
scid <- readr::read_delim("./pubchem_sid_cid.txt", "\t", col_types = "cc")

# データをマージ 
keg_pubc <- dplyr::left_join(cpd_pub, scid, by = "pubchem_SID")

# KEGG COMPOUND IDの冗長性を確認
kegg_cpdup <- keg_pubc[unlist(rsko::dup(keg_pubc$pubchem_CID)),] %>%
  left_join(., cpdat, by = "cid")

7.1 KEGG REACTION IDを取得

# REACTION idと説明
rn <- system("wget -q -O - http://rest.kegg.jp/list/rn/", intern = T, wait = T)
rndat <- data.frame(RID = sapply(strsplit(rn, "\t"), "[", 1), 
                    Description = sapply(strsplit(rn, "\t"), "[", 2),
                    stringsAsFactors = F)

7.2 代謝マップ上に乗っているKEGG REACTION ID取得

• KEGG REACTIONをMAP idで検索する場合IDを変換する必要がある e.g. map00010 -> rn00010
  • wget -q -O - http://rest.kegg.jp/list/pathway/
  • wget -q -O - http://rest.kegg.jp/link/rn/rn00010+rn00020

# 全てのpathway idを取得する 570 pathways
pi <- system("wget -q -O - http://rest.kegg.jp/list/pathway/", intern = T, wait = T) 


# idの変換 
pw_dat <- do.call(rbind, strsplit(pi, "\t")) %>%
  data.frame(., stringsAsFactors = F) %>%
  setNames(., c("map_id","description")) %>%
  mutate(map_id = sub("path:map", "rn", map_id))


# 10 idずづ結合
chunk <- function(x, n){
  nchnk <- ceiling(seq_along(x)/n)
  split(x, nchnk)
}
rnids <- lapply(chunk(pw_dat$map, 10), paste, collapse = "+")

# map idと reaction_idとの対応(reactionがないmapにはNA)
l_rct_id <- lapply(seq_along(rnids), function(i){
  com <- paste0("wget -q -O - http://rest.kegg.jp/link/rn/", rnids[[i]])
  system(com, intern = T, wait = T)
})


map_rids <- do.call(c, l_rct_id) %>% 
  strsplit(., "\t") %>%
  {data.frame(do.call(rbind, .))} %>%
  setNames(., c("map_id", "Rids")) %>%
  mutate(map_id = sub("path:", "", map_id),
         Rids = sub("rn:", "", Rids)) %>%
  group_by(map_id) %>%
  summarise(Rids = paste(Rids, collapse = "; ")) %>%
  right_join(pw_dat, ., by = "map_id")

7.3 REACTION データベースを取得

• wget -q -O - http://rest.kegg.jp/get/R00001
• タグの種類ENTRY, NAME, DEFINITION, EQUATION, REMARK, COMMENT, RCLASS, ENZYME, PATHWAY, MODULE, ORTHOLOGY, REFERENCE
• タグの種類によっては複数行におよぶものがある. ;で結合する
• エントリによっては記載の無いタグもあるので, NAを代入する.

# ユニークなreaction id
# rn <- system("wget -q -O - http://rest.kegg.jp/list/rn/", intern = T, wait = T)
urid <- sub("rn:", "", sapply(strsplit(rn, "\t"), "[", 1))

# KEGG REACTIONのダウンロード(一部, 10づつ50エントリー)
chunk <- function(x, n){
  nchnk <- ceiling(seq_along(x)/n)
  split(x, nchnk)
}
urid_chnk <- lapply(chunk(urid[1:50], 10), paste, collapse = "+")

l_react <- vector("list", length(urid_chnk))
for (i in seq_along(urid_chnk)) {
  com <- paste0("wget -q -O - http://rest.kegg.jp/get/", urid_chnk[[i]])
  l_react[[i]] <- system(com, intern = T, wait = T)
}

# ベクトルの要素を特定の要素を基に分割
vsplt <- function(x, splt){
      pos <- grep(splt, x)
      splt_x <- split(x, findInterval(seq_along(x), pos))
      res <- lapply(splt_x, function(x) x[!x %in% splt])
      res[!sapply(res, identical, character(0))]
}
react <- vsplt(unlist(l_react), "///")

# check
react[[1]]

## [1] "ENTRY       R00001                      Reaction"          
## [2] "NAME        polyphosphate polyphosphohydrolase"            
## [3] "DEFINITION  Polyphosphate + n H2O <=> (n+1) Oligophosphate"
## [4] "EQUATION    C00404 + n C00001 <=> (n+1) C02174"            
## [5] "ENZYME      3.6.1.10"                                      
## [6] "DBLINKS     RHEA: 22455"

7.4 REACTIONデータベースの構文解析とテーブル形式への変換

# # tagの種類
# taglist <- lapply(react, function(x) {v <- sapply(strsplit(x, " "), "[", 1); v[v!=""]})
# tags <- unique(unlist(taglist))

tags <- c("ENTRY", "NAME", "DEFINITION", "EQUATION", "ENZYME", "DBLINKS", "COMMENT",
          "RCLASS", "PATHWAY", "ORTHOLOGY", "MODULE", "REMARK", "REFERENCE")

# 各タグの情報を入れるベクトル
ENTRY <- character(); NAME <- character(); DEFINITION <- character();
EQUATION <- character(); REMARK <- character(); DBLINKS <- character();
COMMENT <- character(); RCLASS <- character(); ENZYME <- character();
PATHWAY <- character(); MODULE <- character();  ORTHOLOGY <- character();
REFERENCE <- character()

# 各レコードを処理
for (j in 1:length(react)) {
  ent <- react[[j]]
  
  # 指定タグが何行目にあるか (どのタグが複数行あるか)
  ent_tags <- sapply(strsplit(ent, "[ ]"), "[", 1)
  extag <- ent_tags[ent_tags != ""] # そのエントリ持っているタグの種類
  extagst <- which(ent_tags != "") # それぞれのタグの開始行
  extaglen <- diff(c(extagst, length(ent) + 1)) # それぞれのタグの記載行の長さ
  extaged <- extagst + extaglen - 1 # それぞれのタグの終了行
  
  # 記載無しのタグはNAを代入
  n.tag <- match(tags, extag)
  if (any(is.na(n.tag))) {
    eval(parse(text = paste0(tags[is.na(n.tag)] , "[[j]] <- NA;")))
  }

  # 各タグごとに文字列を編集
  for (i in 1:length(extaglen)) {
        lines <- ent[extagst[i]:extaged[i]]
        lines <- sub(paste0(extag[i], "[ ]+"), "", lines, fixed = F)
        lines <- gsub("^[ ]+|[ ]+$","", lines)
        # タグ別の編集
        if (extag[i] == "ENTRY") {
          Val <- unlist(strsplit(lines, split = "[ ]+"))[1]
        } else if (extag[i] == "NAME") {
          Val <- paste(lines, collapse = ";")
        } else if (extag[i] == "DEFINITION") {
          Val <- paste(lines, collapse = ";")
        } else if (extag[i] == "EQUATION") {
          Val <- paste(lines, collapse = ";")
        } else if (extag[i] == "REMARK") {
          Val <- paste(lines, collapse = ";")
        } else if (extag[i] == "DBLINKS") {
          Val <- paste(lines, collapse = ";")
        } else if (extag[i] == "COMMENT") {
          Val <- paste(lines, collapse = ";")
        } else if (extag[i] == "RCLASS") {
          Val <- paste(lines, collapse = ";")
        } else if (extag[i] == "ENZYME") {
          Val <- paste(lines, collapse = ";")
        } else if (extag[i] == "PATHWAY") {
          Val <- paste(lines, collapse = ";")
        } else if (extag[i] == "MODULE") {
          Val <- paste(lines, collapse = ";")
        }else if (extag[i] == "ORTHOLOGY") {
          Val <- paste(lines, collapse = ";")
        }else if (extag[i] == "REFERENCE") {
          Val <- paste(lines, collapse = " ")
        }
        # タグリストに代入
        eval(parse(text = paste0(extag[i], "[[j]] <- Val")))
  }
}


# データフレームに変換
react_dat <- data.frame(
  cbind(ENTRY, NAME, DEFINITION, EQUATION,
        REMARK, COMMENT, RCLASS, ENZYME,
        PATHWAY, MODULE,ORTHOLOGY, DBLINKS, REFERENCE), 
  stringsAsFactors = F
  )

7.5 EQUATION(反応式)の取り出し.

• 反応式の書式が一貫していないので注意する.
  • (n-1) C00001, C00046(n+m), 2 C01330(in), 2 C01330(out)

# 酵素反応式を取り出す
options(stringsAsFactors = F)

# 左辺を分解
Lt <- react_dat %$% 
  {setNames(strsplit(sapply(strsplit(EQUATION, " <=> "), "[", 1), " \\+ "), ENTRY)} %>%
  reshape2::melt(.) %>% 
  mutate_if(is.factor, as.character) %>%
  setNames(., c("Cid", "Rid"))

Lc <- lapply(strsplit(Lt$Cid, " "), function(x){
    if (length(x) != 1) {
      gsub("\\(in\\)", "", x)
    } else if (length(x) == 1) {
        if (grepl("\\(.*\\)$", x) & !grepl("\\(in\\)$|\\(out\\)", x)) {
          ai <- gsub("\\(|\\)", "", stringr::str_extract(x, "\\(.*\\)$"))
          xi <- gsub("\\(.*\\)", "", x)
          c(ai, xi) 
        } else if (grepl("\\(in\\)$", x)) {
          xi <- gsub("\\(in\\)", "", x)
          c(1, xi)
        } else {
          c(1, x)
        }
    }
  }) %>%
  {data.frame(do.call(rbind, .), stringsAsFactors = F)} %>%
  setNames(., c("n", "Cid")) %>%
  data.frame(Rid = Lt$Rid, ., equation = "left", stringsAsFactors = F) 

# 右辺を分解  
Rt <- react_dat %$% 
  {setNames(strsplit(sapply(strsplit(EQUATION, " <=> "), "[", 2), " \\+ "), ENTRY)} %>%
  reshape2::melt(.) %>% 
  mutate_if(is.factor, as.character) %>%
  setNames(., c("Cid", "Rid"))

Rc <- lapply(strsplit(Rt$Cid, " "), function(x){
    if (length(x) != 1) {
      gsub("\\(in\\)|\\(out\\)", "", x)
    } else if (length(x) == 1) {
        if (grepl("\\(.*\\)$", x) & !grepl("\\(in\\)$|\\(out\\)", x)) {
          ai <- gsub("\\(|\\)", "", stringr::str_extract(x, "\\(.*\\)$"))
          xi <- gsub("\\(.*\\)", "", x)
          c(ai, xi) 
        } else if (grepl("\\(in\\)|\\(out\\)$", x)) {
          xi <- gsub("\\(in\\)|\\(out\\)", "", x)
          c(1, xi)
        } else {
          c(1, x)
        }
    }
  }) %>%
  {data.frame(do.call(rbind, .))} %>%
  setNames(., c("n", "Cid")) %>%
  data.frame(Rid = Rt$Rid, ., equation = "right", stringsAsFactors = F)