井上潤：LocalBlast

Local BLAST

2015 年 5 月 12 日　改訂
井上潤

Local Blast は，NCBI が提供している web version の Blast search を，ローカル (個人の PC) で行うソフトウェアです．私は local Blast は他のソフトウェアと違って特別なものかと思っていましたが，普通のソフトウェアとほとんど同じように使うことができます．こちらのページからダウンロードできます．

ここで紹介している従来の local Blast (legacy Blast と呼ばれるようになったみたいです) をよりスピードアップし，使いやすくした BLAST+ が出たようです．トップページを見ると，今後legacy blast に新しい機能を付け足して行くことはないようです．まだわかりませんが，将来的には BLAST+ にとって変わられるかも知れません．少なくともインストールは BLAST+ の方がずっと簡単です [2019 年 6 月]．

結果の判断

Blast search の結果は，主にScore, e-value, アライメントで判断します (下図)．

Score
クエリ配列 (検索に用いる配列) と見つかった配列で作成したアライメントの質を表します．高い程よいです．スコアの計算は難しいらしいです．こちらを参照しました．

Evalue
E-value は，「現在のデータベースにおいて，まったく偶然に同じスコアになる配列の数の期待値」，だそうです．小さいほどクエリと見つかった配列の類似性が高いと言えます．E-value が小さいほど偶然には起こりえないことを示します．つまり，E-value が小さい場合には，互いの配列の類似性が高いことになります．
　ちなみに，ヒトのタンパク質転写遺伝子の塩基配列をクエリにして哺乳類全体から遺伝子を集める場合，e^-15 を目安にとりあえず Blast 検索している論文が多いです．一般的な現象かわかりませんが，e value を緩くしたら (e^-5)，解析速度が極端に遅くなりました．

アライメント
クエリ配列と見つかった配列のアライメントです．

準備: .ncbirc ファイル

ソフトウェアのダウンロード (統合 TV)
Mac であれば，"blast-2.2.2X-universal-macosx.tar.gz" をこちらからダウンロードして解凍してください．executables/release フォルダ内です．
* 私にはなぜかわからないのですが，この ftp サイトはフォルダに入ると接続解除になってしまうことがあります．このような場合は，ブラウザを FireFox にかえたり，office や家など場所を変えてダウンロードを試みてください．

設定 (統合 TV)
[もしかしたら，最近の local blast は .ncbirc を設定する必要はないかもしれません．新しい Mac に 2.2.25 をダウンロードしたら，PATH の設定だけで動きました (2011 年 6 月)]

ホームディレクトリに ".ncbirc" を作成し，以下の内容を保存する．

[NCBI]
Data=/Volumes/Seagate/bio_Blast/blast/data
[BLAST]
BLASTDB=/Volumes/Seagate/bio_Blast/blast/db

"Data=" の後には，"BLOSUM62" などのデータ (ダウンロードして得られたもの) が入っているアドレスを書きます．アドレスは pwd で調べてください．
"BLASTDB=" の後には，"nr_fol" などのデータ (ダウンロードして得られたもの) が入っているアドレスを書きます．この設定をしておくことで，例えば，

blastall -p blastn -d /Volumes/Seagate/bio_Blast/blast/db/nt_fol/nt_all/nt -i query.fasta -l sequences.gi -o 1_human.out

の，-d オプションを，"-d nt_fol/nt_all/nt" と省略することができます．

次に，アプリケーションの入っているフォルダにパスを通します．シェルが bash の場合は，".bash_profile" に以下の行を付け加えてください．もちろん，bin フォルダのあるアドレスを pwd で調べて下さい．

export PATH=$PATH:/Users/inouejun/bio/blast/blast-2.2.20/bin

データベースのフォーマット

小さなデータベースで練習
formatdb というソフトウェアを使います．こちらに入って，db/FASTA から mito.aa.gz をダウンロード・解凍してください．作成された mito.aa を開くと，fasta 形式のファイルがまとめて入っていることがわかります．これを

formatdb -p T -i mito.aa

としてデータベース化します．mito.aa.psq, mito.aa.pin, mito.aa.phr というファイルが作成されます．これをデータベースとして指定して blast search に用います．

blastall -p blastp -d mito.aa -i queryfile -v 1 -e 1e-50 -m 9 -o outfile

注意1: フォーマットした後で，作成されたデータベースを他のディレクトリに移すと解析が動かないことがあります．この場合は，移動先のディレクトリでフォーマットをやり直す必要があります．

注意2: Mac の場合，formatdb はダウンロードして得られた bin フォルダに入っていて，コンパイルせずにそのまま使えます．2.2.23+ には formatdb が入っていなかったので，2.2.20 のものを使っています．

NCBI・塩基配列データベース
DDBJ などのデータベースを一気にダウンロードするには，少なくとも HD に 50 GB 程度のスペースが必要だと思います．2010 年 11 月の nt データベースは 12 GB でした．

nt.00.tar.gz ～ nt.07.tar.gz，合計 8 種類 (Version によって数は異なります) のファイルをこちらからダウンロードします．ダウンロードの後はすべてのファイルをダブルクリックによって解凍するだけで，フォーマットが自動的に行われます．ファイルが大きいためでしょうか，解凍には割と時間がかかります．

8 種類のデータベースを一気に検索するには，工夫が必要なようです．どうも，解凍して得られたすべてのファイルを一つのフォルダにまとめる必要があるみたいです．ここでは nt_all というフォルダに解凍して得られたファイルをすべて集めています．そして，

blastall -p blastn -d nt_fol/nt_all/nt -i 3human.fast -o All_human.out

を実行したら，すべてのデータベースを自動的に検索してくれたようです．Helpdesk によると，-d nt によって，nt.nal エイリアスを介してすべてのファイルが読まれるようですが，解凍後に得られたファイルを一つのフォルダにまとめるような指示は得られていませんでした (つまり私独自の工夫なので，注意してください)．

NCBI・アミノ酸配列のデータベース
nr.00.tar.gz ～ nr.03.tar.gz，合計 4 種類のファイルをこちらからダウンロードします．上述にある塩基配列のデータベース同様，ダウンロードの後はすべてのファイルをダブルクリックによって解凍するだけで，フォーマットが自動的に行われます．nr データベースは 2010 年 11 月の時点で 11 GB でした．

4 種類のデータベースを一気に検索するには，やはり工夫が必要なようです．解凍して得られたすべてのファイルを一つのフォルダにまとめます．ここでは nt_all というフォルダに解凍して得られたファイルをすべて集めています．

blastall -p blastp -d nr_fol/nr_all/nr -i query_3human.fast -l sequences.gi -b 0 -o b0NonGenomeProjects.out

オリジナルのデータベース: NCBI からダウンロード
ファスタ形式のファイルを一つにまとめておけば，オリジナルのデータベースを作成して blast search を行うことができます．実際，NCBI のデータベースには FASTA というフォルダがあり，ファイルはすべてファスタ形式で配布されています．

q というファスタ形式のファイルをデータベース化すると，q.fa.nhr, q.fa.nin, q.fa.nsq という 3 つのファイルが新たに作成されます．

./formatdb -p F -i q.fa

./
この場合は formatdb がクエリ配列 q.fa と同じディレクトリに入っています．

-p
F で塩基配列，T でアミノ酸配列のデータベースをフォーマットします．ここをよく間違えます．エラーメッセージが出るので，注意してください．

-i
q.fa は配列の入ったファスタ形式のデータベースです．データベースは NCBI のサイトからダウンロードして使いますが，ファイルが fasta になっていれば，Ensembl やあるいは独自のファイルをデータベースとして使うことが可能です．

　例えば，NCBI のサイトで種名をクエリに検索することで，その種だけのデータベースをダウンロードすることができます (ただし，gi number を使った方法の方がずっと効率が良いです．以下「特定の種/分類群から配列を集める」を参照してください)．
　Search のポップアップタグを Protein にして「Canis_lupus_familiaris[Organism]」とし，イヌのアミノ酸配列配列をファスタ形式で一気にダウンロードします (Nucleotide は時間がかかるのであまりやりませんが，EST は382618件をダウンロードするのに 5 分位でした．Protein もそれほど時間がかかりません)．それを上記の手順によってデータベースのフォーマットを行います．

注意: タンパク質のデータを Fasta 形式でダウンロードしたところ，> 種名の行はあったものの，配列が空であった場合がありました．その場合は空の行を削除することで，フォーマットが行きました．エラーメッセージは以下です．空の行はエディターで「>.*?$\r\r」をクエリに検索すると見つかりました．

[inouejun:Canis_lupus_familiaris]$ formatdb -p T -i Dog41612NCBIprot.fasta
[NULL_Caption] WARNING: Cannot add sequence number 2067 (lcl|2067_Dog41612NCBIprot.fasta) because it has zero-length.
[NULL_Caption] FATAL ERROR: Fatal error when adding sequence to BLAST database.

DDBJ のデータベース
ある程度分類群ごとに別れているので，NCBI のデータベースよりも使いやすい場合があります．DDBJ の表紙のページにある FTP ( ftp.ddbj.nig.ac.jp ) から，ftp を使ってファイルをダウンロードしてきます．例えばアミノ酸配列であれば，ddbj_database > dad フォルダに，分類群ごとに分けられたアミノ酸配列のデータベースがあります．fasta と GenBank format の両方が利用できるようです．

実際にプログラム blastall を走らせる

基本的な使い方

blastall -i q.fa -d db.fa -p blastn > blast_output.txt

blastall
プログラムです．

-p
クエリ配列やデータベースが，塩基配列かアミノ酸配列かによって異なります．こちらをご覧下さい．

-i
q.fa がクエリ配列を入れたファイル名です．

blast_output.txt
アウトプットファイルです．

より細かい設定を行う

blastall -p blastp -d nr_fol/nr_all/nr -i query_testhuman.fast -l sequences.gi -b 0 -e 1e-12 -m 9 -o output.txt

-b 数字
出力に含まれるアライメントされた配列数を制限します．

-e 実数
E-value の閾値を指定します．「-e 0.1」であれば，E-value が 0.1 以下でヒットした配列のみが出力されます．「-e 2e-20」は 2 x 10^-20 を表します．「-e 1e-12」は 1 x 10^-12 を表します．

-v 数字
検索結果リストとして出力する数を指定します．

-m 9
結果をタブ区切りテキストとして出力します．
*-b 0 だと結果が出てきません．

-T T
html 形式の output が出力されます．ブラウザで開くと，accession number がサイトにリンクされています．

以上，詳しくは，統合テレビを見てください．NCBI のオフィシャルな説明はこちらです．

特定の種/分類群から配列を集める

ある種の配列だけ集める
NCBI の website を用いて gi list を作成しておけば，複数のクエリ配列それぞれに最も類似した配列だけを，特定の種についてだけ blast 検索することができます．検索結果が得られたら，そのに書かれている accession number をもとに，fastacmd によって配列を得ることができます．

cDNA データを哺乳類だけから集める
1.
Local Blast に取り込ませるバッチファイルを web version からとってきます．哺乳類，GI number を nucleotide か protei database からダウンロードします:

塩基配列:　　　http://www.ncbi.nlm.nih.gov/sites/entrez?db=nuccore
アミノ酸配列:　http://www.ncbi.nlm.nih.gov/sites/entrez?db=protein

以下のクエリを Box に書いて GO を押してください:

biomol_mrna[prop] AND mammalia[orgn] NOT human[orgn]

ここでは，mammalia として哺乳類だけに絞り，NOT homo としてヒトの番号を含めないようにしています．詳しくは Entrez Help をご覧下さい．とくに，Writing Advanced Search Statements にある Table 3 をクリックすると，必要な情報がまとまっています．

*注意:
biomol_mrna[prop] の効力は種によって異なります．この検索方法は NCBI の help desk に聞きました．確かにヒトやネズミなどアノテーションが進んでいる種では，得られる gi number から余分なレコードを除外してくれるようです．しかしマイナーな種では，実際に Ensembl ですでに全ゲノムがデータベース化されているタンパク質遺伝子の数をこの検索方法で調べると，NCBI のデータベースからはほんの数十しか遺伝子が引っかかってこないことがよくあります．2009 年の 9月の時点では，ヒトやネズミ，チンパンジー以外には適用できないでしょう．おそらく，塩基配列は登録されていても，NCBI で Ensembl の情報がデータベースに入っておらず，タンパク質遺伝子としての情報が登録されていない場合がかなりあるのだと，私は解釈しています (2009 年 9 月)．

Display の横にあるメニューを "Summary" から "GI list" に変更します．同様に "Send to" を "file" に変更して，GI list をダウンロードします．sequences.gi というファイルがダウンロードされているはずです (2009 年 3 月の時点で，626048 ヒット，5.6 MB でした)．

2.
ターミルを開いて，Local Blast にかける配列が入ったファイル (ここでは query.fasta) があるフォルダに移ります．以下のように "-l sequences.gi" をコマンドに組み込みます．query.fast は検索したい配列です．もちろん，哺乳類全体から cDNA を集めたい場合は，ヒトなどの cDNA 配列などをクエリとするのが妥当です．

blastall -p blastn -d nt_fol/nt_all/nt -i query.fasta -l sequences.gi -o 1_human.out

参考になるサイト:

http://210.86.230.110/bioinfo/material/20070821_blastalone/blast_practical_nlm.pdf
http://www.ncbi.nlm.nih.gov/entrez/query/static/help/Summary_Matrices.html
http://fasta.bioch.virginia.edu/cshl/pdf/ncbi_part2_110708.pdf

fastacmd: Accession number を使って sequence を取り出す

fastacmd は Accession number と領域がわかっていれば，目的とする塩基配列を出力します．

fastacmd -d nt_fol/nt_all/nt -s NM_002107.3 -L 116,526

-d
データベース

-s NM_002107.3
Accession number．"-s NP_000240,NP_024931"として複数の検索できます．

-L 116,526
該当する配列の切り出し．

fastacmd -d nr_fol/nr_all/nr -i BlastResList > out

-i
gi number をリストアップしたファイルを指定することもできます.

注意: ">" で示される name line に同じ配列が登録されている gi 番号が列挙されます．

リンク
詳しくは，こちらの「4.2.2 Specific sequence and subsequence retrieval」をご覧下さい．
NCBI の公式マニュアル
 RCC

fastacmd: sequence の分類情報を得る

(以下のように書きましたが，あまり良い方法とは言えない気がします．NCBI の Taxonomy - Site Guide に taxonomy database をダウンロードする方法が書かれています (FTP: NCBI Taxonomy)．これを試すのも手ですが，いずれにしろこちらで手動で得られるような細かい分類情報は得られないようです．[2015 年 1 月])

fastacmd で得られる配列には，登録が古いと生物情報が ">" で示される name line に出てこない場合があります．この場合，-T オプションで分類情報だけを得ることになります．

taxdb.tar.gz をダウンロードして解凍すると，taxdb フォルダの内部に taxdb.bti と taxdb.btd という 2 種類のファイルができます．これを Blast のダウンロードで得られた bin, data, db, doc のうちの db フォルダに入れておきます．すると，.ncbirc に記載した [BLAST] のアドレスが自動的に読まれます．

fastacmd -d nr_fol/nr_all/nr -T -s 197102708

-T: 分類だけを得るオプション．
-s: gi (accession) 番号．

以下のようにクエリとした gi number の分類情報だけでなく，同じ配列が得られている分類群が網羅されます．ここでは gi number 「197102708」で検索したところ，他に 2 レコードが出てきました．

[Jun:human20_8e70]$ fastacmd -d nr_fol/nr_all/nr -T -s 197102708
NCBI sequence id: gi|197102708|ref|NP_001125664.1|
NCBI taxonomy id: 9601
Common name: Sumatran orangutan
Scientific name: Pongo abelii

NCBI sequence id: gi|75055033|sp|Q5RAS5.1|TSN6_PONAB
NCBI taxonomy id: 9601
Common name: Sumatran orangutan
Scientific name: Pongo abelii

NCBI sequence id: gi|55728792|emb|CAH91135.1|
NCBI taxonomy id: 9601
Common name: Sumatran orangutan
Scientific name: Pongo abelii

注意: NCBI にある Local Blast 用のデータベースは web version に比べて更新が遅いです．このためデータが充実していないことがよくあります．例えば，Web version では出てくるのに，gi number をもとに種名や配列を local blast で調べても情報が得られないことが頻繁にあります．このような場合は Bioperl を使って，すべての情報が書いてある GenBank format をダウンロードした方が良いです．

アミノ酸の Accession number/配列から塩基配列を得る．

blastn と gi list を用いる
一つのやり方として，アミノ酸配列をクエリ配列として，塩基配列データベースを local blast によって検索する方法が考えられます．

blastall -i queryAA.fas -d NCBI/nt/nt -o output.txt -p tblastn -l gi -e 1 -b 2 -m 9

-p tblastn
アミノ酸配列をクエリとして塩基配列データベースを検索．

-l gi
対象種の塩基配列 gi number リストを，NCBI のウェブサイトを用いてあらかじめ作成して，gi というファイルに保存してある．

-e 1
E-value が 1 以下の配列だけリストアップ．

-b 2
結果を 2 つだけ表示させる．

-m 9
結果をタブ区切りにする．

上記の方法で，アミノ酸配列をクエリ配列として local Blast を使って cDNA 配列を自動的に得る perl script を書きました．こちらです．

Local Blast では GenBank format を得ることはできない
Local blast にアミノ酸の Accession number を与えても，直接塩基配列の情報を得ることはできません．また，Local blast を使って Accession number から GenBank format を得ることはできません (help desk に問い合わせました)．help desk によると efetch を使う必要があるそうですが，私はどのように使ったら良いのかわかりませんでした．

This is NOT possible - blast database contains only FASTA sequences.

You will need to use the GI and efetch to get the sequences. For
information on how to use efetch, see:
http://www.ncbi.nlm.nih.gov/entrez/eutils/

BioPerl を使って GenBank フォーマットを得る
アミノ酸配列のGenBank フォーマットを得ます．

#!/usr/bin/perl -w

use Bio::DB::GenBank;
use Bio::SeqIO;

$db_obj = Bio::DB::GenBank->new;
$seq_obj = $db_obj->get_Seq_by_acc("68270954");
$seqio_obj = Bio::SeqIO->new(-file => '>gi68270954.gbk', -format => 'genbank' );
$seqio_obj->write_seq($seq_obj);

get_Seq_by_acc　の部分は 68270954 のような gi number でも，AAY88974.1 のような GenBank の number でも得られる結果は同じです．得られた GenBank フォーマットを見ると，CDS の欄に「/coded_by="NM_002107.3:116..526"」のように，コードしている塩基配列の Accession number と領域が明示されています．

#!/usr/bin/perl -w

use Bio::DB::GenBank;
use Bio::Species;

$db_obj = Bio::DB::GenBank->new;
$seq_obj = $db_obj->get_Seq_by_acc("68270954");

print "display_id\t", $seq_obj->display_id, "\n";
print "desc\t", $seq_obj->desc, "\n";
print "display_name\t", $seq_obj->display_name, "\n";
print "accession\t", $seq_obj->accession, "\n";
print "primary_id\t", $seq_obj->primary_id, "\n";
print "seq_version\t", $seq_obj->seq_version, "\n";
print "keywords\t", $seq_obj->keywords, "\n";
print "length\t", $seq_obj->length, "\n";
my $species_obj = $seq_obj->species;
print "species\t", $species_obj->binomial(), "\n";
print "seq\t", $seq_obj->seq, "\n";

こちらのウェブサイトを参照させていただきました．BioPerl のチュートリアルが勉強になるみたいです．

Local Blast 解析用のスクリプト

いくつかスクリプトを書きました．ほとんどのスクリプトは Local Blast や mafft，pal2nal などが自身の PC にインストールされていて PATH が切られていないと動かないです．

ある種の gi number を集める
こちらをご覧下さい．BioPerl を使っています．あまり難しくありません．

ある種の塩基，アミノ酸，EST データベースを検索して配列を得る
種名を与えると，3 つのデータベースを検索してヒットした配列を新たにアライメントに加えます．あらかじめ cDNA アライメントが必要です．逆鎖にヒットした配列は相補鎖に変換されてアライメントされます．mafft や pal2nal などが自身の PC で PATH が切られた状態で使える必要があります．000_commands.pl の最初にあるコマンドを打つだけで動くはずです．未だ難解なプログラムなので，これから改善してゆきます．
BLseqRetrieverNCBIdb_fol.tar.gz [2010 年 12 月]．

Perl script: ID の配列を DB からとってくる

Blast の結果で得られた ID の配列を DB からとってくるスクリプトを書きました．
seqRet_fol.tar.gz

Perl script: 二つの protDB を比較して対応する ID リストを作る

例えば，古い DB の ID が変更されていた場合，最新の DB で対応する配列の ID が必要になることがよくあります．このような場合に，対応関係を示したリストと配列が得られます (result.txt)．このプログラムは local Blast を使うので，解凍して得られたフォルダに blastall (ここでは Version 2.2.20 を使用) と formatdb を入れてください．その後，

perl corresFinder.pl

と入力してください．データベースも自動的にフォーマットされます．
corresFinder.tar.gz
[2011 年 11 月]

リンク集

バイオインフォマティクス入門

ダウンロードとコンパイルは，これをみながらやりました．

Local blast の使い方

まずはこれをみた方が良いです．第 1 回～第 5 回まで，すべて役に立ちます．これに限らず，統合 TV シリーズはとても優れています．

NCBI, Learning Center (NCBI 公式 HP)

細かいパラメータ設定について．

blastall (NCBI 公式 HP)

上のサイトの一部．blastall に関するパラメーター設定について．

ダウンロードとコンパイル (NCBI 公式 HP)

ダウンロードできるデータベースの表もあります．

MacOSX 設定の細かい説明 (NCBI 公式 HP)

少しわかりにくいかもしれません．

BLAST output の fasta file を取得する

統合 TV シリーズ

ウィキペディア

Blast について，なかなか良い説明が書いてあります．

Protein vs nucleotide seqch with tblastn

Entrez Help

gi number を得るのに必要な情報が書いてあります．例えば，哺乳類だけの遺伝子を検索したいときに，gi number として，まずこちらから哺乳類だけの gi number を集めます．このときに，「哺乳類だけ」あるいは「遺伝子の長さ」などといった条件をどのように設定するか書いてあるページです．

GenBank flat file

GenBank の flat file をダウンロードする ftp サイトです．README.GenBank に初歩的な説明があり，gbrel.txt にそれぞれのファイルの説明があります．もとのページはこちらです．

Kerfeld and Scott (2011)

このページでは，NCBI の Helpdesk や UCL の友人に教えていただいた情報も紹介しています．教えてくださった皆様に感謝しています．