日本東部ソフトウェア　バイオインフォマティクス情報　関連情報

ここではバイオインフォマティクスの現場で使われるツールや、トピックスを扱っていきます。初回は配列比較ツールのデファクトスタンダードであるBlastの導入方法をご紹介します。

◆Blastとは

BlastとはNCBI(National Center for Biotechnology Information)が公開しているツールです。 Basic Local Alignment Search Toolの略で、任意の配列同士を比較するためのモジュール群を含みます。例えば、blastnというプログラムは自身が保有する塩基配列情報(=質問配列)が、ヒトゲノム (=データベース配列)のどこに位置するのかを割り出すことなどに用いられます。いまヒトゲノムと書きましたが、比較対象にするデータはユーザーが任意に設定可能です。

◆Blastのモジュール

Blastは比較する配列の種類によってプログラムを使い分ける必要があります。
以下のようなプログラムがあります。
blastn ・・・塩基配列対塩基配列の比較を行うためのプログラム。
blastp ・・・アミノ酸配列対アミノ酸配列の比較を行うためのプログラム。
blastx ・・・塩基配列(質問配列)対アミノ酸配列(データベース配列)の比較を行うためのプログラム。
tblastn ・・・アミノ酸配列(質問配列)対塩基配列(データベース配列)の比較を行うためのプログラム。
tblastx ・・・塩基配列対塩基配列の比較を、アミノ酸配列に翻訳して行うためのプログラム。
これらのプログラムはblastallというモジュールを実行するときに指定します。

また、Blastは通常1回の実行で質問配列は1種類しか実施することは出来ませんが、megablast というプログラムを使うことで複数種の質問配列を一度に処理することができます。

他にもゲノム全体のような大きな配列の中で質問配列がどこに位置するのかではなく、自身が保有する質問配列同士を比較するためのモジュールとしてbl2seqがあります。

今回は研究室用サーバーや自宅サーバーなどのLinux環境でBlastを使えるようにしたいときの手順をご紹介します。
以下、Linux環境で(bシェルで)作業する前提で書きます。

1.Blastの入手

NCBIの公式サイトから入手できます。 Blastには通常のblastとインターネット経由で実行できるblastがありますが、通常版をダウンロードします。通常版のblastでも、OSごとに別々に管理されているので、自分の環境にあったものを選択します。ここではLinuxの32bit版("linux-ia32"の"blast")をダウンロードします。約22.6MBあります。

DL元URL:
ftp://ftp.ncbi.nlm.nih.gov/blast/executables/LATEST/blast-2.2.12-ia32-linux.tar.gz
(2005.9.30現在)

2.作業用フォルダの作成と展開

Blast環境をつくりたいLinuxマシンにログインします。つぎにFTPなどを使ってダウンロードしたblastを、とりあえず適当なディレクトリに格納します。ここでは自分用の一般ユーザのホームディレクトリに格納することにします。
つぎに、blastのファイルをツールプログラムを格納するディレクトリに移動させます。ツールプログラムを格納するディレクトリは、決まった場所があるならそこに移動させてください。
特に無い場合、自分のホームディレクトリにいるものとして、そこに"tools"というディレクトリを作成し、toolsで管理することにします。

$ mkdir tools

と入力し、ディレクトリを作ります($は注記が無ければ一般ユーザプロンプトを意味します。以後同様とします。…要は、Linuxのコマンド自体は"mkdir"以後です。ユーザプロンプトというコトバがわからなかったら、調べてみてください)。

つぎに

$ mv blast-2.2.12-ia32-linux.tar.gz ./tools/

と入力し、blastを移動させます。そして

$ cd tools/

と入力してblastの移動先に移動します。つぎに

$ tar zxvf blast-2.2.12-ia32-linux.tar.gz

と入力してください。圧縮されていたファイルが展開されます。

$ ls

と入力したら、元ファイルのほかにblast-2.2.12というディレクトリが出来ているはずです。ここのなかに移動します。

$ blast-2.2.12/bin/

と入力してください。移動したらそのディレクトリ内に"blastall"があるはずです。確認のために実行してみましょう。

$ ./blastall

と入力してみてください。ダウンロード時などにファイルが壊れていなかったら、 blastのコマンドオプション情報がずらっと表示されるはずです。

ここまで完了したら、blastの本体の準備はできました。次に、Linuxのコマンドと同様にblastallを使えるようにパスを通してしまいましょう。

3.blastallのパス

ここまでの作業でblastallの本体は準備が出来ました。ですが、このままではほかのLinuxコマンドのように使えません (使う度に毎回blastallがあるディレクトリまで指定して入力しなければなりません。例えば「blastall」と入力するだけで使いたいのに、「/home/suzuki/tools/blast-2.2.12/bin/blastall」などと毎回入力しなければならない。これはあまりに面倒です。)。

パスの通し方は少しLinuxに慣れた方ならご存知と思うのですが、例えば急に研究室で専用のLinuxサーバを用意し、他から干渉されない形でblast環境を構築し始めたばかりの方も居られるかもしれません。念のため書いておきます。

仮に、自分が「suzuki」というユーザで、blastallの本体を「/home/suzuki/tools/blast-2.2.12/bin/blastall」に格納してあるとします。

(1)自分だけ設定したい場合

自分のホームディレクトリ直下の「.bashrc」というファイルにパスを追記します。追記する内容は2行

#blast 20050930
PATH="$PATH":/home/suzuki/tools/blast-2.2.12/bin/

です。1行目はコメント、2行目はパスです。追記して保存したら

$ cd
$ source .bashrc

とコマンド入力してください。コマンドの1行目が自分のホームディレクトリに戻るコマンドで、 2行目が追記した内容を有効化するコマンドです。

(2)みんなで使えるように設定したい場合

「/etc/bashrc」というファイルにパスを追記します。追記する内容は2行

#blast 20050930
PATH="$PATH":/home/suzuki/tools/blast-2.2.12/bin/

です。1行目はコメント、2行目はパスです。

またこの場合は念のため、blastallの実行の権限を変更しておきましょう。

$ chmod 755 /home/suzuki/tools/blast-2.2.12/bin/blastall

とコマンドを入力してください。ここまでの作業で、新規にログインするすべてのユーザが「blastall」の入力でblastallを実行できるようになっています。

4.blast用データベースの準備

blastは本体のみでは使えません。例えばDNA対DNAの比較を行う場合は、DNAのデータベースを準備する必要があります。ここでいうデータベースとはOracleやPostgreSQLなどのDBMSではなく、 DNAやアミノ酸配列の文字情報のファイルをblastパッケージ内のツールで加工したファイル群をさします。

(1)元配列の入手

まず、NCBIのFTPサイトなどから目的の配列情報を取得します。例えば、ヒトのBuild35.1の第22染色体のDNA情報を例に話を進めます。

配列情報はNCBIで公開されています。
ftp://ftp.ncbi.nih.gov/genomes/H_sapiens/CHR_22/
(2005.9.30現在)

ここからファイル名に"mfa"という文字列を含むファイルを取得します。つまり今回の目的のファイルは
ftp://ftp.ncbi.nih.gov/genomes/H_sapiens/CHR_22/hs_ref_chr22.mfa.gz
(2005.9.30現在)
です。

他の動物で構築したい場合は、ftp://ftp.ncbi.nih.gov/genomes/に移動して探してください。例えば、M. musculus(マウス)の第18番染色体のDNA情報なら
ftp://ftp.ncbi.nih.gov/genomes/M_musculus/CHR_18/mm_ref_chr18.mfa.gz
になるでしょう。

ここで取得するmfaファイルの"mfa"はMalti Fastaの略です。 1つのファイル中にFasta形式のデータが複数件はいっています。

Fasta形式とは、1行目に配列の情報(配列の名称やゲノム中の位置、ほか)があり、 2行目以降に配列情報を記述した書式です。

(2)元配列の格納と展開

取得した元配列を格納するディレクトリを作成します。このディレクトリがデータベースのディレクトリになります。

今回はユーザ「suzuki」のホームディレクトリ直下に「hs_ncbidb/B35.1」というディレクトリを作成します。

ユーザsuzukiでログインし、ホームディレクトリで

$ mkdir -p hs_ncbidb/B35.1

とコマンドを入力します。そして、上で取得したhs_ref_chr22.mfa.gzをB35.1ディレクトリへ移動させます。その後にこのディレクトリに移動してファイルを展開します。

$ mv hs_ref_chr22.mfa.gz hs_ncbidb/B35.1/
$ cd hs_ncbidb/B35.1/
$ gunzip hs_ref_chr22.mfa.gz

とコマンドを入力します。 1行目がファイルの移動。2行目がユーザプロンプトの移動。3行目がファイルの展開です。

(3)formatdbの実行

blastallで使うデータベースの準備には「formatdb」を使います。 formatdbはblastallと同じ場所に格納されているので、blastallのパスを通している場合、コマンドのように使えます。

今回はこのformatdbを

$ formatdb -i hs_ref_chr22.mfa -p F -o T

とコマンド入力して実行します。
実行するときのオプションですが、「-i」は入力ファイルの指定です。「-p」は入力ファイルの配列の種類の指定でアミノ酸配列のときT、塩基配列のときFと指定します。今回は塩基配列なので「-p F」と指定します。「-o」はインデックスファイル作成の指定です。blast用のデータベースを準備するときはTを指定します。
この他のオプションを調べたいときは

$ formatdb --help

と入力してください。

formatdbを実行すると、「.nin」「.nni」「.nnd」「.nhr」「.nsd」「.nsi」「.nsq」の7種類の拡張子のファイルが作成されます。拡張子の前は入力ファイルと同じです。ここまでの作業でDNA対DNAのblast用のデータベースは準備できました。

実際に使うときは、データベースの指定個所で上記の7種類の拡張子より前の共通する部分までを指定します。今回の場合は具体的には「/home/suzuki/hs_ncbidb/B35.1/hs_ref_chr22.mfa」を指定します。

また、人間のような情報量が多い生物のmfaファイルは染色体単位で公開されています。ところがblastを実行したいときは全染色体に対して実行したい場合もあります。そのような時は例えば以下のようなスクリプトを書いて、染色体単位のmfaファイルをひとつのファイルにまとめ、 formatdbを実行してデータベースを準備してください。

スクリプト例:
#! /bin/sh

cat hs_ref_chr1.mfa > hs_ref_chrALL.mfa
cat hs_ref_chr2.mfa >> hs_ref_chrALL.mfa
cat hs_ref_chr3.mfa >> hs_ref_chrALL.mfa
(中略)
cat hs_ref_chr22.mfa >> hs_ref_chrALL.mfa
cat hs_ref_chrX.mfa >> hs_ref_chrALL.mfa
cat hs_ref_chrY.mfa >> hs_ref_chrALL.mfa

5.blastの実行

ここまでで準備したblastallを用い、DNA対DNAのblastである"blastn"を実行してみます。

調べたい配列をFasta形式で準備します。通常は実験などで配列情報が取得できているのでしょうが、今回は実験はしておりません。従い、NCBIから取得した"ABCD1P4"という遺伝子の配列を使って実行してみます。

fasta形式のデータ例:
>gi|29789831|ref|NG_001265.1| Homo sapiens ATP-binding cassette(ABCD1P4)
CCCTGCCTAGGTTGGGAGGCTATGTGTGACTGGAAGGATGTCCTGCCGGGTGGCAAGAAG
CAGAGAATCGGCATGGCCTGCATGTTCTACCACAGGTGAGCACTCAGGGCTGGCAGGCTC
CCTGGGGTCCCCTGGAAGGAGAAGTAGCAGCTGTGGGGAGGCCTGGGCTCAGTGGAGCCT
GAGCCGGGTTGGGGTGTTGGGCCCTGGAGGGTGCACAGACTCTCCTCTAGGCCCGGACCC
CCAGGCCCAAGTACACCCTCCTGGATGAATGCACCAGTGCCGTGAGCATCGACGTGGAAG
GCAAGATCTTCCAGGCGGCCAAGGACGCAGGCATTGCCCTGCTCTCCATCACCCACCGGC
CCTCCCTGTGGTAGGTGCCCTGTCTCCCTNCCTGGGGTGAGTGGGAGTGGCTGCCTGAGG
GGAGGAGGTGGCCTGTTGGGCCAGGCGGCAGCAGCAGGCGGCTGTCGTCAGCAGCCCTCG
TGCCGTGCCCCTGACCCTGTCCCTCTCCTGGCCAGGGAGTACCACACACACTTGCTACAG
TTCGATGGGGAGGGCGGCTGGAAGTTCGAGAAGCTGGACTCAGCTGCCCGCCTGAGTCTG
ACGGAGGAGAAGCAGCGGCTGGAGCAGCAGCTGGCAGGCATTCCCAAGATGCAGCGGCAC
CTCCAGGAGCTCTGCCAAATCCTGGGCGAGGCCGTGGCCCCAGCGCATGTGCCGGCACCT
AGCCCGCAAGGCCCTGGTGGCCTCCAGGGTGCCTCCACCTGACGCCACCCTCCCCAGCCC
CTGCCCCGCC

FASTA形式ですが、blastに供する場合は1行目の書式は「行頭に">"をつける」以外は決まり事はありません。また、2行目以降の改行は必須ではありません。また、2行目以降の配列データは大文字と小文字で結果に違いはありません。

テストデータとして以下の内容で実行してみます。

fasta形式のデータ例その2:
> ABCD1P4 1to50
CCCTGCCTAGGTTGGGAGGCTATGTGTGACTGGAAGGATGTCCTGCCGGG

この2行のみのデータで「test_abcd1p4.txt」という名前でファイルに保存します。そして、以下のコマンドを入力します。

$ blastall -p blastn -d ./hs_ncbidb/B35.1/hs_ref_chr22.mfa -i test_abcd1p4.txt > test_abcd1p4.txt.out

「-p」は実行プログラムのオプションです。DNA対DNAの比較の場合は「blastn」を指定します。
「-d」はデータベースの指定のオプションです。ここまでの作業で作成したデータベースを指定します。
「-i」は入力ファイルのオプションです。入力ファイルのパスを指定します。

出力は何も指定しない場合は標準出力されます。そのため、上の記述ではリダイレクトで「test_abcd1p4.txt.out」に出力しています。例えば出力が大きくなるのが予想される場合は「>」以降を「| gzip > test_abcd1p4.txt.out.gz」とすれば圧縮形式で保存できます。因みに「-o」オプションでも出力ファイルを指定できます。

この他にもいくつかオプションが指定可能です。詳細はblastallを何のオプションも指定しないで実行したときにヘルプ情報が表示されますので、そちらを参照してください。

6.appendix

blast実行結果例:

BLASTN 2.2.12 [Aug-07-2005]

Reference: Altschul, Stephen F., Thomas L. Madden, Alejandro A. Schaffer, 
Jinghui Zhang, Zheng Zhang, Webb Miller, and David J. Lipman (1997), 
"Gapped BLAST and PSI-BLAST: a new generation of protein database search
programs",  Nucleic Acids Res. 25:3389-3402.

Query= ABCD1P4 1to100
         (100 letters)

Database: hs_ref_chr22.mfa 
           12 sequences; 34,972,128 total letters

Searching............done

                                                                 Score    E
Sequences producing significant alignments:                      (bits) Value

ref|NT_011519.10|Hs22_11676 Homo sapiens chromosome 22 genomic c...   190   1e-48
ref|NT_011520.10|Hs22_11677 Homo sapiens chromosome 22 genomic c...    34   0.15 
ref|NT_011521.3|Hs22_11678 Homo sapiens chromosome 22 genomic co...    30   2.3  
ref|NT_019197.4|Hs22_19353 Homo sapiens chromosome 22 genomic co...    28   9.2  
ref|NT_011525.6|Hs22_11682 Homo sapiens chromosome 22 genomic co...    28   9.2  
ref|NT_011523.10|Hs22_11680 Homo sapiens chromosome 22 genomic c...    28   9.2  

>ref|NT_011519.10|Hs22_11676 Homo sapiens chromosome 22 genomic contig
          Length = 3661581

 Score =  190 bits (96), Expect = 1e-48
 Identities = 99/100 (99%)
 Strand = Plus / Plus

                                                                         
Query: 1     ccctgcctaggttgggaggctatgtgtgactggaaggatgtcctgccgggtggcaagaag 60
             |||||||||||||||||||||||||||||||||||||| |||||||||||||||||||||
Sbjct: 23002 ccctgcctaggttgggaggctatgtgtgactggaaggacgtcctgccgggtggcaagaag 23061

                                                     
Query: 61    cagagaatcggcatggcctgcatgttctaccacaggtgag 100
             ||||||||||||||||||||||||||||||||||||||||
Sbjct: 23062 cagagaatcggcatggcctgcatgttctaccacaggtgag 23101



 Score = 34.2 bits (17), Expect = 0.15
 Identities = 17/17 (100%)
 Strand = Plus / Plus

                                
Query: 24      gtgtgactggaaggatg 40
               |||||||||||||||||
Sbjct: 3573407 gtgtgactggaaggatg 3573423


(中略)


>ref|NT_011523.10|Hs22_11680 Homo sapiens chromosome 22 genomic contig
          Length = 2589913

 Score = 28.2 bits (14), Expect = 9.2
 Identities = 14/14 (100%)
 Strand = Plus / Plus

                           
Query: 24    gtgtgactggaagg 37
             ||||||||||||||
Sbjct: 83698 gtgtgactggaagg 83711



 Score = 28.2 bits (14), Expect = 9.2
 Identities = 14/14 (100%)
 Strand = Plus / Minus

                            
Query: 23     tgtgtgactggaag 36
              ||||||||||||||
Sbjct: 981706 tgtgtgactggaag 981693


  Database: hs_ref_chr22.mfa
    Posted date:  Sep 29, 2005  2:57 PM
  Number of letters in database: 34,972,128
  Number of sequences in database:  12
  
Lambda     K      H
    1.37    0.711     1.31 

Gapped
Lambda     K      H
    1.37    0.711     1.31 


Matrix: blastn matrix:1 -3
Gap Penalties: Existence: 5, Extension: 2
Number of Hits to DB: 2275
Number of Sequences: 12
Number of extensions: 2275
Number of successful extensions: 37
Number of sequences better than 10.0: 6
Number of HSP's better than 10.0 without gapping: 10
Number of HSP's successfully gapped in prelim test: 0
Number of HSP's that attempted gapping in prelim test: 0
Number of HSP's gapped (non-prelim): 37
length of query: 100
length of database: 34,972,128
effective HSP length: 16
effective length of query: 84
effective length of database: 34,971,936
effective search space: 2937642624
effective search space used: 2937642624
T: 0
A: 0
X1: 11 (21.8 bits)
X2: 15 (29.7 bits)
S1: 12 (24.3 bits)
S2: 14 (28.2 bits)