macでインフォマティクス

macでインフォマティクス

HTS (NGS) 関連のインフォマティクス情報についてまとめています。

nrなどのNCBIデータベースをダウンロードする ncbi-blast-dbs

2018 12/10 タイトル訂正

2020 9/7追記2020 9/11 わかりにくい説明を修正、9/11 簡単な並列処理例追記、9/12.9/15  taxonkit コマンド修正, わかりにくい部分を修正、10/3 taxonkitのコマンドを修正、10/10 コマンド微修正、10/28 誤字修正

2022 1/6 例のパラメータ(e-value)修正、07/10  refseq proteinの例を追加

2023/10/17 追記

 

 ncbi-blast-dbsはデータベースファイルを並行してダウンロードすることで、NCBIのデータベースをローカルに用意するのにかかる時間を短縮する。使用するスレッド数は自動的に決定される。 MD5チェックサムが検証され、ダウンロード時にデータベースボリュームが抽出される。 データベースボリュームは特定の順序でダウンロードされない。 新しいバージョンがサーバー上で使用可能な場合はボリュームが更新され、破損している場合は再ダウンロードされる。 中止されたダウンロードは安全に再開される。

 

 

インストール

mac os10.14でテストした。

依存(一部)

本体 GIthub

gem install ncbi-blast-dbs

#md5sumコマンドがないならhomebrewで導入
brew install md5sha1sum

#あるいはdockerのrubyのイメージを使って導入(rubyオフィシャルイメージだと失敗した)
docker run -itv $PWD:/data --rm -w /data circleci/ruby bash
$ gem install ncbi-blast-dbs

 > ncbi-blast-dbs -h

$ ncbi-blast-dbs 

 

refseq_protein, 16SMicrobial, cdd_delta, env_nr, env_nt, est, est_human, est_mouse, est_others, gss, gss_annot, htgs, human_genomic, nr, landmark, nt, other_genomic, pataa, patnt, pdbaa, pdbnt, ref_prok_rep_genomes, ref_viroids_rep_genomes, ref_viruses_rep_genomes, refseq_genomic, refseq_rna, refseqgene, sts, swissprot, taxdb, tsa_nr, tsa_nt, vector, est_human_blob, taxdump

そのまま叩くとダウンロードできるデータベース一覧が表示される。 

 

実行方法

例えばローカル環境にBLAST検索データベースを構築するため、全データベースをダウンロードする。

#NTは核酸(nucleotide collection)、NRはタンパク質(non-redundant translated sequecne)
ncbi-blast-dbs nt nr

#refseq protein DB(解凍後150GBほど。ダウンロード時には圧縮ファイルの解凍が行われるため、500GB近く空きが必要)
mdkir Refseq_protein_DB && cd Refseq_protein_DB/
ncbi-blast-dbs refseq_protein
#ダウンロードが終われば、BLASTPサーチのDBとしてそのまま使用できる。(.tar.gzは不要)
cd ../
blastp -db Refseq_protein_DB/refseq_protein -query query.fasta -outfmt 7 -out out.txt -num_threads 20 -evalue 1e-1

 

 

NCBIは、FTPサーバからダウンロードする際、メールアドレスを送ることになっている(らしい)。メールアドレスを送りダウンロードするには、コマンドを以下のように変える。

email="my email address here" ncbi-blast-dbs nr

 

どのようなデータベースがダウンロードできるかと、blastの使い方についてはこちらのチュートリアルをみて下さい。

 

追記

実際にローカルでblastサーチを行う例

得られたデータベースを使ってblastサーチを行い、生物種の情報を得るまでの流れをまとめた。

nt.00.tar.gzだけを使うなら、解凍したnt.00/データベースを以下のように指定。

blastn -db nt.00/nt.00 -query query.fasta -outfmt 6 -out out.txt -num_threads 20 -evalue 1e-100

 

1、分割されたデータベースを全て指定したい場合、リンク先の記事のようにループさせるか、分割ファイルを1箇所に集めて"-db nt"と指定する。 

#step1 NTなら全てのnt分割ディレクトリを解凍。00から52まであるなら
for i in 00 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52; do tar -xvf nt.${i}.tar.gz; done

#step2 分割ファイルを集める
mkdir nt-database
mv nt* nt-database/
#不要なら分割tar.gzを消す
rm nt*tar.gz

#step3a blasting (クエリの配列数が多くて結果が肥大しそうなら、タブ区切りで1行1出力)
blastn -db nt-database/nt -query query.fasta -outfmt 6 -out out.txt -num_threads 20 -evalue 1e-10

#step3b blasting (クエリの配列が数個程度なら"-html"をつけてHTML出力。アラインメントまで見る)
blastn -db nt-database/nt -query query.fasta -out out.txt -num_threads 20 -evalue 1e-10 -html

#step3c クエリの配列が大量にあり、かつ出力後にベストヒットだけを確実に取り出すなら"-outfmt 7"にしてベストヒット順に出力(-outfmt 6はベストヒット順になってない)
blastn -db nt-database/nt -query query.fasta -outfmt 7 -out out.txt -num_threads 20 -evalue 1e-10

#3c => step4 それからawkなどで最初の行を取り出す(参考サイト)。sortを使うと文字とE-value(指数)の同時ソートでつまづくので、 -outfmt 7出力して下の方がおすすめ。
#create header
echo qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore |sed s/' '/$'\t'/g > best_hits.txt
#then, extract best hit
cat out.txt | awk '/hits found/{getline;print}' | grep -v "#" >> best_hits.txt

#step5 必要に応じてcutで必要な列だけ取り出す
cut -f 1,2,11 best_hits.txt > best_hits_truncated.txt

 

2、blast検索が終わったら、得られたGenbank accession IDテーブルから、lineage情報を得る(リンク先の補足)。そのためのデータベースを作る(2は初回のみの作業)。

#prepare database (parepare database step1)

nucl_*accession2taxid.gzをダウンロードする。

ftp://ftp.ncbi.nih.gov/pub/taxonomy/accession2taxid/

> wget ftp://ftp.ncbi.nih.gov/pub/taxonomy/accession2taxid/nucl_*accession2taxid.gz && gzip -dv *gz

次のコマンドを走らせ、ダウンロードしたファイルを統合して1つのデータベースにする。

#prepare database (parepare database step2)
#accession2taxidを結合、”sort -u”で重複を除きつつソート
cat nucl_gb.accession2taxid nucl_wgs.accession2taxid dead_nucl.accession2taxid dead_wgs.accession2taxid |sort -u > concatenated_accession2taxid 

#prepare database (parepare database step3)
#1列目とGI numberは不要(2列目はバージョンつき。".1"ならバージョン1)
cut -f 2,3 concatenated_accession2taxid > genbank_AC_db

 

3、ルックアップする表が得られたので、この表に対してripgrep検索してGenbak IDからtaxIDを取り出す。最後にtaxonkitでlineageに変換。

#step6 blast結果をIDだけに整形(step5と重複しているが2列目だけ取り出している)
cut -f 2 best_hits.txt > blast_best-hit_genbankID_list

#step7 複数IDをループ処理.echoで1列目にクエリのID、2列目に見つかったIDを出力
cat blast_best-hit_genbankID_list | while read line
do
echo -n -e "$line\t" >> taxid_list ;rg -m 1 $line genbank_AC_db |cut -f 2 >> taxid_list
done

#step8 taxid => lineage (taxonkit使用)。2列目がtaxid(-i 2
taxonkit lineage -i 2 taxid_list >> lineage_list

#step8 alternative
#step8のtaxonkitは入力リストが多すぎるとほとんどの行を無視してしまうことがある。その場合は分割して処理(* コメントをいただいてます。確認してください)
split -l 1000 taxid_all taxidsplit_
ls taxidsplit_* > file_list

cat file_list | while read line; do
taxonkit lineage -i 2 $line >> lineage
done

 

step7 alternativeblast出力が大きい場合、grep(ripgrep) 検索のstep7が律速になる。その場合、分割して少し工夫するだけで並列処理できる。

#step7 alternative1 
#7-1 splitする。blast best hit出力が10万行あり、1000ずつに分割したいなら
split -l 1000 blast_best-hit_genbankID_list split_

#step7 alternative2 listを作る
ls split_* > file_list

#step7 alternative3 20ずつ並列処理。ここではparallelを使う。(*1)
cat file_list | parallel -j 20 '
cat {} | while read line; do
echo -n -e "$line\t" >> taxid_{.}
rg -m 1 $line genbank_AC_db | cut -f 2 >> taxid_{.}
done
'

#step7 alternative4 統合する
cat taxid_* > taxid_all
=> step8へ

 

コメント

ダウンロードしたデータベースが認識されない時は、ダウンロードに部分的に失敗しているか(md5を確認)、インストールしたblastのバージョンがデータベースのバージョンと対応していない可能性があります。blastを最新版に更新してみて下さい。

 引用

GitHub - yeban/ncbi-blast-dbs: Fast download BLAST databases from NCBI.

 

関連ツール


 

 

 

*1

TR3990xを用いてテストしたところ、64並列処理でも問題なく動作した。

(メモリ128GB、ubuntu18.04環境)

 

 

メモ

ダウンロード直後にも関わらずデータベースとのバージョンが異なるというエラーが出る場合、blastコマンドの方が古い可能性が高い。最新の安定板blastをインストールすれば解決する。

> conda remove blast

> conda install -c bioconda -y blast

 

 

追記

このブログ記事とは関係ありませんが、blastでmax_target_seqs 1を指定する場合の挙動について調べた結果が論文の補足データとしてまとめられています。

”アプリケーションによっては、「ベストヒット」とは、E値が最も低い配列、ビットスコアが最も高い配列、同一性の割合が最も高い配列、データベースやクエリシーケンスのカバレッジが最も高い配列のことを指します。BLASTがどのように定義されているかに関わらず、「ベストヒット」を正しく識別することを期待するのは不合理です。”