macでインフォマティクス

macでインフォマティクス

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

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

2018 12/10 タイトル訂正

2020 9/6

2020 9/7追記2020 9/11 わかりにくい説明を修正

2020 9/11 簡単な並列処理例追記

2020 9/12.9/15  taxonkit コマンド修正, わかりにくい部分を修正

2020 10/3 taxonkitのコマンドを修正

2020 10/10 コマンド微修正

 

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

 

ncbi-blast-dbsに関するツイート

 

 

インストール

mac os10.14でテストした。

依存(一部)

本体 GIthub

gem install ncbi-blast-dbs

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

 > 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検索データベースを構築するため、全データベースをダウンロードする。

ncbi-blast-dbs nt nr

 

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"と指定する。 

#1 NTなら全てのnt分割ディレクトリを解凍。00から27まであるなら
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; do tar -xvf nt.${i}.tar.gz; done

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

#3a blasting (クエリの配列が結果が多いならタブ区切り出力)
blastn -db nt-database/nt -query query.fasta -outfmt 6 -out out.txt -num_threads 20 -evalue 1e-10

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

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

#4 それから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

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

 

2、blast検索が終わったら、得られたGenbank accession IDテーブルから、lineage情報を得る(リンク先の補足)。ということで、まずはデータベースを作る。

#prepare database ((parepare database step1)

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

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

ダウンロードしたファイルを統合して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に変換。

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

#7 複数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

#8 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 alternative 
#7-1 splitする。blast best hit出力が10万行あり、1000ずつに分割したいなら
split -l 1000 blast_best-hit_genbankID_list split_

#7-2 listを作る
ls split_* > file_list

#7-3 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
'

#7-4 統合する
cat taxid_* > taxid_all
=> step8へ

 

 引用

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