NCBIで全データを一度にblast解析し、得られたリストをEntrez Directでアノテーションに変換する。

2020 10/9 リンク追加

複数の配列のblast解析を行う場合、ローカルでデータベースなどを構築して進めるのが一つの手である。しかしローカルだとデータベースの更新や、データサイズが問題になる（例えばnrのデータも2015年にダウンロードすると200GBを超えていた）。

ネットワーク越しのランで十分な速度が得られる限りは、専門のアノテーターが管理している最新のデータベースを使うのがベストである。そのためにコマンドラインからネットワーク越しにblast解析を行う方法も用意されているが、特定のポートを開ける必要があり、環境によってはやや敷居が高い。今回はNCBIにプラウザ越しにアクセスして、数百のcontigのblast解析を行う方法を書いていく。

後半では、全blast結果をテキストでダウンロードし、それからbest hitのIDを取り出し生物名を得る流れをまとめている。ダウンロードされたファイルにはヒットした生物のEntrez_IDしか書かれていないので、生物名を得るためにUNIXのコマンドとNCBIのツールを組み合わせている。

生物系研究者の方はSafariなどのブラウザを使ってNCBIにアクセスされている人が多いと思います。その時、全配列を同時にblastにかけて結果をダウンロードする方法を知っていると便利です。一例ですが参考にしてみてください。

8/20 追記

NCBIのNucleotide blastに移動する。

Nucleotide blastを選択。

f:id:kazumaxneo:20170810175809j:plain

ウィンドウ下のファイルを選択、をクリック。

生物を限定するならorganismの項目に学名を書く。

de novoアセンブルして得たcontigファイルを選択。

f:id:kazumaxneo:20170810180112j:plain

BLAST解析を実行する。

途中でページ更新に失敗した場合、１ページ戻り右上のRecent resultsを選択（recent resultsリンク）。

f:id:kazumaxneo:20170820135533j:plain

IDをクリックすれば再開できる。

f:id:kazumaxneo:20170820135722j:plain

ジョブが終わっている場合、右端の downloadでダウンロード可能（xmlファイルのみ）。

数分で結果が得られた。

f:id:kazumaxneo:20170810180231j:plain

上の方のダウンロードをクリック。

f:id:kazumaxneo:20170810180301j:plain

Hit Table(text) を選択してダウンロード。

ダウンロードされたファイルにはhitした全ての部位の情報がタブ形式で保存されている。

usr $ head alignment.txt

# blastn

# Iteration: 0

# Query: NODE_1_length_40202_cov_39.6477

# RID: SS3HY4K7014

# Database: nr

# Fields: query id, subject ids, query acc.ver, subject acc.ver, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score

# 264 hits found

NODE_1 gi|545742605|emb|FR873486.1| NODE_1 FR873486.1 83.535 2976 407 41 9041 11995 12939 15852 0.0 2704

NODE_1 gi|545742605|emb|FR873486.1| NODE_1 FR873486.1 84.876 2215 295 22 34599 36796 838 3029 0.0 2198

NODE_1 gi|545742605|emb|FR873486.1| NODE_1 FR873486.1 85.036 1517 197 19 12592 14100 17670 19164 0.0 1517

hitしたIDのところはEntrezのIDとembのIDだが、このままでは何がヒットしたのか分かりにくい。IDをコンバートしてわかりやすく修正する。

query_IDとsubject_IDだけ取り出す。

cut -f 1-2 Alignment.txt > subject_and_query_ID

head -15 subject_and_query_ID

# blastn

# Iteration: 0

# Query: NODE_1_length_40202_cov_39.6477

# RID: SS3HY4K7014

# Database: nr

# Fields: query id, subject ids, query acc.ver, subject acc.ver, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score

# 264 hits found

NODE_1 gi|545742605|emb|FR873486.1|

１つのcontigに複数当たっている。ここではbest hitだけで十分なので、uniqコマンドを使い１列目のNODE名が重複する行を消す。awkにパイプしてコメント行の削除もワンライナーで行う。

uniq -f 2 subject_and_query_ID | awk '! /^#/' > subject_and_query_ID_unique

user$ head -5 subject_and_query_ID_unique

NODE_1 gi|545742605|emb|FR873486.1|

NODE_2 gi|332139268|gb|HQ009547.1|

NODE_3 gi|531033897|emb|HF586473.1|

NODE_4 gi|190702157|gb|EF710641.1|

NODE_5 gi|54109776|emb|AJ632317.1|

重複が除けたので、ここからsubject_IDだけ保存する。

cut -f 2 subject_and_query_ID_unique > subject_ID_unique

user$ head subject_ID_unique

gi|545742605|emb|FR873486.1|

gi|332139268|gb|HQ009547.1|

gi|531033897|emb|HF586473.1|

gi|190702157|gb|EF710641.1|

gi|54109776|emb|AJ632317.1|

gi|332139193|gb|HQ009535.1|

gi|531033897|emb|HF586473.1|

gi|531034048|emb|HF586474.1|

gi|531033897|emb|HF586473.1|

Entrez IDだけ保存する。

cut -f 2 -d "|" subject_ID_unique > Entrez_ID

user$ head Entrez_ID

545742605

332139268

531033897

190702157

54109776

332139193

531033897

531034048

531033897

IDリストが入手できた。

ここから先はDAVIDの変換ツールを使っても良いが、今回はEntrez Directのコマンドを使う。ターミナルで以下のコマンドをペーストしてツールをダウンロード。

cd ~
  perl -MNet::FTP -e 
    '$ftp = new Net::FTP("ftp.ncbi.nlm.nih.gov", Passive => 1);
     $ftp->login; $ftp->binary;
     $ftp->get("/entrez/entrezdirect/edirect.tar.gz");'
  gunzip -c edirect.tar.gz | tar xf -
  rm edirect.tar.gz
  export PATH=$PATH:$HOME/edirect
  ./edirect/setup.sh

ホームにedirect/ができる。パスを通しておく。

echo 'export PATH=$PATH:/Users/user/edirect/' >> ~/.bash_profile
source ~/.bash_profile

efetchコマンドでxmlファイルを出力し、xmlからxtractコマンドで<tile>を抽出する。IDを１つだけ探すなら以下のようなコマンドになる。

efetch -id ID -db nucleotide -format docsum | xtract -pattern DocumentSummary -element 'Title' #IDのところはEntrez_IDの番号を入力。

実行例

user$ efetch -id 545742605 -db nucleotide -format docsum | xtract -pattern DocumentSummary -element 'Title'

Cotesia congregata bracovirus segment 28, complete sequence

前半の出力は下のようなxmlファイルである。

f:id:kazumaxneo:20170820132227j:plain

Titleの項目が生物名に使えそうである。

f:id:kazumaxneo:20170820132237j:plain

この行のタグの間をEntrez Directのxtractコマンドで抽出している。

今回はEntrez_IDのリストに対して実行するので、ループで順番に処理させる。

while read line; do title=$(efetch -id $line -db nucleotide -format docsum |  xtract -pattern DocumentSummary -element 'Title'); echo "$line      $title"; done<Entrez_ID > title_list.txt

user$ head -5 output

545742605 Cotesia congregata bracovirus segment 28, complete sequence

332139268 Cotesia vestalis bracovirus segment c24, complete sequence

531033897 Cotesia congregata sequence containing Cotesia congregata bracovirus proviral locus 2 (PL2)

190702157 Cotesia sesamiae Mombasa bracovirus clone BAC 6L7, complete sequence

54109776 Cotesia congregata virus complete genome, segment Circle14

IDから<Title>部分の情報に変換できた。

このような流れで、全contigをblastにかけてベストヒットのアノテーションリストを入手できる。

集計する。1列目は数値順、２列めは辞書順でソート。

sort -k 1,1n -k 2 title_list.txt |uniq -c |sort -r > ranking.txt

591 1228866639 Papilio machaon Ycf2-like protein

7 1228866639 PREDICTED: Helianthus annuus

7 1228808943Medicago truncatula cytochrome b559

2 1044533238 PREDICTED: Vigna angularis LINE-1

省略

上のような集計結果を出力できた。大半のcontigはPapilio machaonにbest hitしていた。

今回はリストを入手しbest hitだけ残してblast解析を行ったが、NCBIのblast結果からxmlファイルもダウンロードできる。ケースによってはそちらを使った方が早い。また行程がやや煩雑なので、頻繁に行うならもう少し工夫してスクリプト化した方が良いと思います。blastp、blastx等も流れは同じです。

追記

https://ncbiinsights.ncbi.nlm.nih.gov/2018/06/14/5-new-youtube-videos-blast-medgen-pubchem-api-keys/

NCBIの中の人のプレゼンテーションのです。投稿日時は2018/05/17。ほぼような同じことをやっていますが、ずっと包括的な説明になっていています。23分あります。詳しく知りたい方はどうぞ。