はじめに

だいぶ時間が経過してしまいましたが、前回、こんな記事を書きました。

pediatricsurgery.hatenadiary.jp

PiphillinなどでpathwayやKEGG orthologyの予測メタゲノムデータを得た後の解析は多々あると思いますが、変動しているpathwayの中で自分に興味深いものがあったとします。その中で、どこの遺伝子が実際に動いているのか？ということを明らかにする方法を記事にしてみたいと思います。

結構ググったんですけど、まとまった情報が出てこなかったのでまとめてみます。間違いや他にもっと良い情報あるだろうということがあれば教えていただきたいです。ここでは例として、胆汁酸代謝に興味があるとして、一次胆汁酸と二次胆汁酸の合成のpathwayである"Primary bile acid biosynthesis"と"Secondary bile acid biosynthesis"について調べていこうと思います。

1. KEGGRESTのインストール

いろいろ調べた結果、どうやらRパッケージの{KEGGREST}を使用するのが良さそうです。KEGGRESTはBioconductorの中に入っているので、biostatisticsなどを参考にインストールしてください。

source("https://bioconductor.org/biocLite.R")
biocLite("KEGGREST")
library(KEGGREST)

2. KEGGRESTに含まれるデータベースを確認する

listDatabases()

[1] "pathway"  "brite"    "module"   "disease"  "drug"     "environ"  "ko"       "genome"   "compound" "glycan"   "reaction" "rclass"  
[13] "enzyme"   "organism"

つまり、KEGGに保管されるデータベースがだいたいこのパッケージからアクセスできるということになると思います。これで全部かは僕は知らないです。今回必要なのは、"pathway"と"ko"です。

3. KEGG pathwayの一覧を取得する

KEGG pathwayの一覧を取得するにはkeggList(“データベース名”)を入れます。データベース名はlistDatabases()で返ってきた上記の値を入れて下さい。pathwayのリストを取得します。とっても多いので10個だけ

pathway <- keggList("pathway")[1:10]
pathway

                             path:map00010                              path:map00020                              path:map00030 
            "Glycolysis / Gluconeogenesis"                "Citrate cycle (TCA cycle)"                "Pentose phosphate pathway" 
                             path:map00040                              path:map00051                              path:map00052 
"Pentose and glucuronate interconversions"          "Fructose and mannose metabolism"                     "Galactose metabolism" 
                             path:map00053                              path:map00061                              path:map00062 
       "Ascorbate and aldarate metabolism"                  "Fatty acid biosynthesis"                    "Fatty acid elongation" 
                             path:map00071 
                  "Fatty acid degradation"

4. 自分に興味のあるpathwayを探す

自分に興味のあるpathwayがみなさんあると思います。また単語からpathwayを検索したいときもあると思います。そんな時はkeggFind()を使えば良いようです。複数の検索ワードも使えるようですが、ここでは"Bile"で検索してみたいと思います。

as.data.frame(keggFind("pathway", c("Bile")))

                keggFind("pathway", c("Bile"))
path:map00120   Primary bile acid biosynthesis
path:map00121 Secondary bile acid biosynthesis
path:map04976                   Bile secretion

という結果が返って来ますので、"Primary bile acid biosynthesis" “Secondary bile acid biosynthesis"のpath:map00120とpath:map00121を使用することにします。

5. 興味あるpathwayの情報をKEGGから抽出する

pathwayの情報をKEGGから抽出する際際はkeggGetを使用します。

# 後々のことを考えて欲しいIDの数字を格納しておきます
PathwayID <- c("00120","00121")
paste("path:map",PathwayID,sep="")
query <- keggGet(paste("path:map",PathwayID,sep=""))
Pathway_Table <-cbind.data.frame(Pathway_ID,)
length(query) # 2つ登録したので2になっていればOK
names(query[[1]])

keggGetの結果はlistで返ってきますので、1で1つ1つの情報にアクセスします。どんな情報が含まれているか調べると

[1] "ENTRY"       "NAME"        "DESCRIPTION" "CLASS"       "PATHWAY_MAP" "KO_PATHWAY"  "REFERENCE"

具体的にこれらの引数にどんな情報が割り当てられているかをチェックします。

query[[1]]$ENTRY # "map0120" 
query[[1]]$NAME # "Primary bile acid biosynthesis"
query[[1]]$DESCRPITION # NULL
query[[1]]$CLASS # Metabolism; Lipid metabolism
query[[1]]$PATHWAY_MAP # "Primary bile acid biosynthesis" 
query[[1]]$KO_PATHWAY # "ko01200"
query[[1]]$REFERENCE

恐らくこの中で有用なのは、query1$CLASSで、Metabolism; Lipid metabolismの中に含まれるpathwayであることがわかるということでしょうか？

6. 対応するKEGG pathwayに含まれるorthologを抽出する

pathwayには数十から数百のorhologueが含まれており、これは各pathwayで共有されているものもあります。例えば今回の"Primary bile acid biosynthesis"と"Secondary bile acid biosynthesis"のpathwayに含まれるKEGG orthologyがどれなのかということを明らかにします。

しかし、この作業に関して、私の調べた範囲内では{KEGGREST}では実行することができませんでしたが、Twitterで質問したところこのような返事をアドバイスをいただきました。

これですかねhttps://t.co/YeMRUoseqL
— nishi (@sai_kai_cs) 2017年5月31日

このページにアクセスして、下記にURLを好きなpathwayのIDに変えれば取り込めることになります。さてこれをRに取り込むわけですわけですが、どういったウェブ上のテーブルに対してどれを使えばいいかよくわかっていませんので、下記を参考に手当たり次第試してみます。

結論としては、{XML}のreadHTMLTable()と{RCurl}のgetURL()はダメで、read.csv()とread.delim()で可能でした。read.csv()とread.delim()ってURL指定できるんですね、個人的にはread.delimが汎用性高くて好きです。

僕はfor構文を使い回さないとコード書けないマンなので、冗長なコードをお許し下さい。そして誰か格好いいコードに直して下さい。

# KO_pathXXXXXに各pathwayに含まれるKOを格納する
for(i in 1:length(PathwayID)){
assign(paste("KO_path",PathwayID[i],sep=""),
       read.delim(paste("http://rest.kegg.jp/link/ko/ko",PathwayID[i],sep=""),
                         header=F)[,2]
)}

> KO_path00120 
 [1] ko:K00037 ko:K00251 ko:K00488 ko:K00489 ko:K00659 ko:K01442 ko:K01796 ko:K07430
 [9] ko:K07431 ko:K07439 ko:K07440 ko:K08748 ko:K08764 ko:K10214 ko:K10223 ko:K11992
[17] ko:K12405 ko:K12408
18 Levels: ko:K00037 ko:K00251 ko:K00488 ko:K00489 ko:K00659 ko:K01442 ... ko:K12408

> KO_path00121
[1] ko:K00076 ko:K01442 ko:K15868 ko:K15869 ko:K15870 ko:K15871 ko:K15872 ko:K15873
[9] ko:K15874
9 Levels: ko:K00076 ko:K01442 ko:K15868 ko:K15869 ko:K15870 ko:K15871 ... ko:K15874

7. 各KOのannotationをチェックする

あとはこれらのKEGG orthologyがどのようなannotationになっているのかをチェックします。再びkeggGet()を使用して抽出していきます。

Table_KO <- NULL

for(j in 1:length(PathwayID)){
# 各PathwayIDに含まれるKOの詳細情報を格納する
KO_path_j <- paste("KO_path",PathwayID[j],sep="")
Detail_KO_path_j <-
  keggGet(eval(parse(text=KO_path_j)))

# KOのそれぞれからDefinition ( annotation ) を抽出していく
# keggGet()の結果はlist
for(i in 1:length(Detail_KO_path_j)){
Table_KO_i <-
  cbind.data.frame(paste("path:",PathwayID[j],sep=""),
                 Detail_KO_path_j[[i]]$ENTRY,
                 Detail_KO_path_j[[i]]$DEFINITION)
Table_KO <- rbind.data.frame(Table_KO,Table_KO_i)
}}
Table_KO <- tbl_df(Table_KO)

colnames(Table_KO) <- c("Pathway","Entry","Definition")

     Pathway  Entry
      <fctr> <fctr>
1 path:00120 K00037
2 path:00120 K00251
3 path:00120 K00488
4 path:00120 K00489
5 path:00120 K00659
6 path:00120 K01442
# ... with 1 more variables: Definition <fctr>

これで各KOのDefinitionが隠れてますが、興味があるpathwayに含まれるKOの一覧を抽出する作業が完了しました。Enjoy!

Note of Pediatric Surgery

腸内細菌、R、ときどき小児外科

{KEGGREST}でKEGG pathwayのどのorthologyに差があるのか調べる