興味のある遺伝子を他のどんな種類の細菌が保有しているかBLASTNと{rentrez}を使って調べる
0. はじめに
ここでは興味を持った細菌の遺伝子、もしくは遺伝子配列を、どんな種類の他の細菌が持っているのかを知りたい時にやってみたことをまとめておきます。
という感じになります。
尚、この記事を書くに当たってはTwitterで@sai_kai_csさんにめちゃくちゃ助けてもらいました。ありがとうございます。
これは、keggのDBを間違えているか(右の箱をnr-ntにしていればOK)、keggのnrのversionの違いかと思います。
— nishi (@sai_kai_cs) 2018年1月9日
nrであるならばどっちもほぼおなじだと思うので、どちらでもよいかと…
例としてフラクトオリゴ糖という糖を分解するのに必要な酵素である beta-fructofuranosidase ( K01193 ) を使っていきます。これはもともとBifidobacteriumなどが多く有することがわかっております。ここでは、論文などを読んでBifidobacterium breve UCC2003というstrainの持つbeta-fructofuranosidaseに興味を持ったとします。そこでBifidobacterium breve UCC2003以外のstrain、もしくは他のspeciesで、beta-fructofuranosidaseの遺伝子を有する細菌がいないか、調べる方法を記載してみます。
この記事では既知の遺伝子の配列を使っていますが、配列さえ同定できていれば、自分で同定した未知の遺伝子でも同じワークフローで検索できると思います。
1. KEGGとNCBIでシーケンスを取得
KEGGで検索するとBifidobacterium breve UCC2003のbeta-fructofuranosidaseは下記であることがわかります。 KEGG T02534: Bbr_0020
Bbr_0020 Glycosyl hydrolases family 32, Beta-fructosidase or sucrose-6-phosphate hydrolase [KO:K01193] [EC:3.2.1.26]
このページの下の方に NT seq というボタンがあり、クリックすると遺伝子の配列が記載されたページに飛びます。そこに記載されているのがこの遺伝子の配列です。
>bbru:Bbr_0020 K01193 beta-fructofuranosidase [EC:3.2.1.26] | (GenBank) Glycosyl hydrolases family 32, Beta-fructosidase or sucrose-6-phosphate hydrolase (N) atgccaaagctgtactaccaattccctggcacgtggttcggcgattgcatgccgttcggc aagggcaacgagttctttctgtatcaccagcgagacaatcgcaatccggaaccattcggc gaaccgttcggctgggatctggccacaacctctgatttcgtgcattacaaggattgcggt gtagccgtgccgcgcggcagggacgacgagcaggatcagttcatattcgcgggatcggtg tttgaaggggaagggcaataccatatcttctacaccgggtacaaccgcgattacccccag cgggacaagccggcgcaggtgctgatgcacgcagtcagcgatgatttgtaccactggacg aaaacccagaacgccatcaccttcactccgcaagaaggctacgatccagacgattggcgc gacccgtgggtcattcgcgacgacgaacatgaccagtatctgctaatccttggagcccgt aagatcggcccgaaaacgcagcagaccggccgcacggtgaagttcacgtctaaagacctc aagaactggacgttcgaaggtgacttctgggctccgaacctatacacgatgcatgagatg cctgatctgttcaagatgggcgagtggtggtatcacatcgttacggaatacagcgatcgg cacaagatggtgtatcgcatggccaagagtctagaaggcccgtggatcgcacccgatgat gacgcgtttgacggcagcgcctattatgccggacggaccttcgaactaggtggacggcga atcctcttcggctgggttgccaccaaggaagactgcgatgacgtaaagaacttcgaatgg gctggcaccttcgtccctcacgaggtctaccagcgtaccgacggcacacttggcgtgcgc attcccgacacggtctgggatgcgttcaacgaacgcgagcaaatcgtgacgaatgcggta atcgacggtactgacgctcgccgtaacgccgtgctggccgagggctgcggagacctgtat tcgtttgaagcggacgtcatgttcagcgagggaacgcgatccttcggtctccgcgtctac cgcgatgaggcaagcgagcgtggctatcagttcatcttcccgattggtgagcatcgctac gtattcgaaggcagcccgaacttcccgtggttcaactgtatgaatatcgggctagagcgg ccaatcgaactggttgccgggcgtgaatatcacattcagctcatcgtggacgacaccatc gccactctgtacgtcaatggcgtcgcactcaacgcgcgtatgtaccagcacccaggcgac gcattaagcatgttcgtcaccgacggcacgctcaccgtcacgaatgcgaccatcagtcgt ggactcaaggaggtgtga
2. BLASTNで検索
NCBI BLASTチュートリアル を参考にしました。
ウェブサイト
これらの配列をBLASTで検索します。BLASTについては、たくさんの良い記事があるので適当にググってください。一応、自分のメモをコチラにまとめてあります。BLASTにはいろいろと種類がありますが、この遺伝子を有する細菌を調べたいので、塩基→塩基を検索する blastn を使用します。
BLAST Searchのウェブページから配列をコピペして検索してもいいですし、上記の NT seq のリンク先のページの下の方の BLAST というところをクリックすると、Sequence dataが入力された状態でBLAST Searchが開くので便利です。blastn にチェックが入っていることを確認して Computeをクリックします。
BLASTはNCBIのサイトからも実行できますが、どちらでもパラメーターを調整すれば結果は同じだそうです。後述しますが、Hit tableをcsvで取得できるという点から、NCBIの方がいいかな？と個人的には思っており、この記事ではNCBIでの操作を記載します。
データベースの選択
https://blast.ncbi.nlm.nih.gov/
データベースはnr/ntを選択します。nr/ntとはNucleotide Collectionで、GenBank、EMBL、DDBJ、PDBに含まれている塩基配列の全エントリ（ HTGS0,1,2、EST、GSS、STS、PAT、WGSをのぞく ）およびRefseq配列から重複を除いたものとあります。
検索エンジンの選択
データベースが巨大化するつれ、blastでも検索時間がかなり長くかかるようになってきたということから、より高速に検索することの可能な検索エンジンがblastnから派生して開発されています。NCBI BLASTチュートリアルによると
# megablast 問い合わせ配列と非常に近い塩基配列のみを高速に検索する 約95%程度の塩基一致率以上のものを検出する # discontiguous megablast megablastよりも多少のミスマッチを許したシードから検索を開始する 生物種間を越えた 程度の検索が可能 # blastn 7塩基以上一致するシードからの検索 網羅性は高い
となっております。ここでは網羅性を高めるため、blastnを選択しておきます。
Algorithm parameters
デフォルトでは上位100件しか表示されないようになっているので、Algorithm parametersをクリックして、Max target sequencesを必要に応じて増やしておくと良いと思います。
Descriptionsの見方
BLASTを開始してしばらくすると結果のページに移ります。こんな感じで表示されていると思います。
それぞれの数字の見方ですが大まかに言うと下記のように説明できるようです。
- Max score
- そのエントリと問い合わせ配列とが複数個所でヒットした場合最大のスコア
- Total score
- そのエントリと問い合わせ配列とが複数個所でヒットした場合、そのスコアの合計
- Query cover
- 問い合わせ配列に対して、ヒットした領域がどれだけカバーするか
- E-value
- エントリと問い合わせ配列とのヒットが偶然起こる期待値
- 小さいほどそのヒットが偶然出ないことを表す
- Ident
- エントリと問い合わせ配列がヒットしている領域での（塩基）一致率
Hit tableのdownload
DescriptionのDownloadをクリックしてHit table(csv)を選んでダウンロードして、Rでの操作に移ります。 ( もちろんウェブページのテキストをコピーして処理してもいいと思います ) ここでは 5AMDVY63015-Alignment-HitTable.csvというファイルをダウンロードできたとします。このcsvなのですがかなり不親切なもので、最初は何が書いてあるのかわかりませんでした。Max scoreやE-value、またAccessionなどはよく見ると記載されていますが、細菌の名前の一覧を取得できていないことに気づきます。ですので、これからRでAccession ( GenBankのID ) から細菌の名前を取得する方法を記載していきます。
3. Rでのワークフロー
データのダウンロード
僕のGitHubのNOPSというリポジトリに 5AMDVY63015-Alignment-HitTable.csvを置いてありますので、これをダウンロードして使います。データ処理をするために{tidyverse}をインストールしておいてください。またcsvファイルには特に列名が記載されていませんので、自分でマニュアルで追加してあげます。
library(tidyverse) RCurl::getURL("https://raw.githubusercontent.com/Razumall/NOPS/master/Data/5AMDVY63015-Alignment-HitTable.csv") %>% read.csv(text = ., header = FALSE) %>% dplyr::select(V2, V4, V5, V13, V14) %>% dplyr::select(`Max score`=V14, `E-value`=V13, Ident=V5, Accession=V4) %>%# GenBankのID tbl_df() -> df
| Max score | E-value | Ident | Accession |
|---|---|---|---|
| 2522.0 | 0.00e+00 | 100.000 | CP023197.1 |
| 2522.0 | 0.00e+00 | 100.000 | CP023196.1 |
| 2522.0 | 0.00e+00 | 100.000 | CP023195.1 |
| 2522.0 | 0.00e+00 | 100.000 | CP023194.1 |
| 2522.0 | 0.00e+00 | 100.000 | CP023193.1 |
| 2522.0 | 0.00e+00 | 100.000 | CP023192.1 |
| 2522.0 | 0.00e+00 | 100.000 | CP021558.1 |
| 2522.0 | 0.00e+00 | 100.000 | CP021557.1 |
| 2522.0 | 0.00e+00 | 100.000 | CP021556.1 |
| 2522.0 | 0.00e+00 | 100.000 | CP021555.1 |
こんな感じになります。ここでE-valueが高いものを除外しておきましょう。閾値の設定はいろいろあるようですが、ここではE-value = 0のもののみを選んでおきます。
df %>% filter(`E-value` == 0 ) -> df_E0
すると26 strainに絞り込めました。
{rentrez}のインストールと使い方
Rの{rentrez}を使用します。ここからはからだにいいものの記事を参考にさせていただきました。そもそもEnterzって初めて聞いたのですが
EntrezはPubMed,GenBank, GEO等のNCBIのデータベースに対する、ユーザー向けに作られたータ取得システムです。
ということみたいです。EnterzをRから動かすのが{rentrez}で、RからからPubMedにアクセスしたりGenBankのレコードにアクセスしたりすることを可能にします。
# インストール install.packages("rentrez") library(rentrez)
データベースの一覧を取得する関数はentrez_dbs()なのですが、これを実行すると
entrez_dbs()
[1] "pubmed" "protein" "nuccore" "ipg" "nucleotide" "nucgss" "nucest" "structure" "sparcle" [10] "genome" "annotinfo" "assembly" "bioproject" "biosample" "blastdbinfo" "books" "cdd" "clinvar" [19] "clone" "gap" "gapplus" "grasp" "dbvar" "gene" "gds" "geoprofiles" "homologene" [28] "medgen" "mesh" "ncbisearch" "nlmcatalog" "omim" "orgtrack" "pmc" "popset" "probe" [37] "proteinclusters" "pcassay" "biosystems" "pccompound" "pcsubstance" "pubmedhealth" "seqannot" "snp" "sra" [46] "taxonomy" "biocollections" "unigene" "gencoll" "gtr"
となります。今回は"genome"を使用しますが、ゲノム情報の取得以外にもいろいろなことに使えそうですね。
GenBankのIDから細菌の名前を取得する
GenBankのIDはAccessionの列に格納されています。一番上にある CP023197.1 について調べてみましょう。GenBank IDから菌名を取得するにはentrez_search()を使います。
entrez_search(db = "genome", term = "CP023197.1")
と入力すると
Entrez search result with 1 hits (object contains 1 IDs and no web_history object) Search term (as translated): CP023197.1[All Fields]
と返って来ます。実際にどの菌名が取得できたのか確認するには
entrez_search(db = "genome", term = "CP023197.1") %>% .$ids %>% # マッチしたidを表示 entrez_summary(db = "genome", id = .) %>% .$organism_name
[1] "Bifidobacterium breve"
とすることで得られます。本当は株名までわかればよかったのですが、自分の今の解析ではspecies levelで十分なため、ここで終わりにします。 ( 中途半端ですみません ) あとはこれを各列に実行すればいいわけです。ただdplyr::mutateで各行ごとに上記を事項してSpeciesの結果を渡そうと思ったのですが、複数行を同時に実行できなかったので、lapplyで一行ずつ実行してリストにして最後にdo.call(rbind.data.frame)で1つのデータに戻すという超面倒くさいことをやっています。解決策がわかる方は教えていただければと思います。
lapply( 1:nrow(df_E0), function(i){ df_E0 %>% .[i, ] %>% mutate( Species = entrez_search(db = "genome", term = Accession ) %>% .$ids %>% # マッチしたidを表示 entrez_summary(db = "genome", id = .) %>% .$organism_name) } ) %>% do.call(rbind.data.frame,. ) -> df_E0_Species
これを実行すると
| Max score | E-value | Ident | Accession | Species |
|---|---|---|---|---|
| 2522 | 0 | 100.000 | CP023197.1 | Bifidobacterium breve |
| 2522 | 0 | 100.000 | CP023196.1 | Bifidobacterium breve |
| 2522 | 0 | 100.000 | CP023195.1 | Bifidobacterium breve |
| 2522 | 0 | 100.000 | CP023194.1 | Bifidobacterium breve |
| 2522 | 0 | 100.000 | CP023193.1 | Bifidobacterium breve |
| 2522 | 0 | 100.000 | CP023192.1 | Bifidobacterium breve |
| 2522 | 0 | 100.000 | CP021558.1 | Bifidobacterium breve |
| 2522 | 0 | 100.000 | CP021557.1 | Bifidobacterium breve |
| 2522 | 0 | 100.000 | CP021556.1 | Bifidobacterium breve |
| 2522 | 0 | 100.000 | CP021555.1 | Bifidobacterium breve |
| 2522 | 0 | 100.000 | CP021554.1 | Bifidobacterium breve |
| 2522 | 0 | 100.000 | CP021553.1 | Bifidobacterium breve |
| 2522 | 0 | 100.000 | CP021389.1 | Bifidobacterium breve |
| 2522 | 0 | 100.000 | CP021394.1 | Bifidobacterium breve |
| 2522 | 0 | 100.000 | CP021391.1 | Bifidobacterium breve |
| 2522 | 0 | 100.000 | CP021388.1 | Bifidobacterium breve |
| 2522 | 0 | 100.000 | CP021387.1 | Bifidobacterium breve |
| 2522 | 0 | 100.000 | CP021385.1 | Bifidobacterium breve |
| 2522 | 0 | 100.000 | CP006716.1 | Bifidobacterium breve |
| 2522 | 0 | 100.000 | CP006712.1 | Bifidobacterium breve |
| 2522 | 0 | 100.000 | CP000303.1 | Bifidobacterium breve |
| 2516 | 0 | 99.928 | CP019596.1 | Bifidobacterium breve |
| 2516 | 0 | 99.928 | CP010413.1 | Bifidobacterium breve |
| 1279 | 0 | 80.640 | AP012331.1 | Bifidobacterium scardovii |
| 1072 | 0 | 77.419 | AP012325.1 | Bifidobacterium catenulatum |
| 717 | 0 | 72.304 | CP011786.1 | Bifidobacterium actinocoloniiforme |
となり、多くの株はBifidobacterium breveがbeta-fructofuranosidase ( Glycosyl hydrolases family 32, Beta-fructosidase or sucrose-6-phosphate hydrolase [KO:K01193] [EC:3.2.1.26] ) を持っているようですが、B.scardoviiやB.catenulatumを持っているものもいるようです。長くなりましたが、こんな感じです。少しでもお役に立てれば。
