配列相同性検索のBLASTをPythonを用いて実行する方法を解説します。ここでは、特にローカルBLASTをPythonを用いて実行する方法の説明です。
バイオインフォマティクス環境
以下の方法で構築した環境であることを前提として説明します。
- Windows10 + WSL (WindowsユーザーのためのPythonを用いたゲノム解析環境)
※ 検証はしていませんが、プログラム自体は基本的にmacOS、Linuxでも動作すると思われます。
ここで使用したアプリケーションのバージョンは以下の通りです。
- Python 3.7.6
- Biopython 1.76 → Biopythonについて
- BLAST 2.9.0
PythonからローカルBLASTを実行する方法
BiopythonのBio.Blast.Applicationsモジュールのクラスを用いることで、PythonからローカルBLASTを実行することができます。なお、あらかじめBLASTをPCにインストールしてパスを通しておく必要があります。もちろん、Anaconda (miniconda)を用いてインストールすればパスも通してくれるので、インストールするだけでOKです。
DBの作成
ローカルBLASTを実行するためには、まず検索対象とするDBを作成する必要があります。BLAST検索用のDBはmakeblastdbコマンドを用いて次のように作成します。
# 塩基配列を指定する場合
makeblastdb -in ファイル名 -dbtype nucl
#アミノ酸配列を指定する場合
makeblastdb -in ファイル名 -dbtype prot
makeblastdbコマンドはBiopythonではNcbimakeblastdbCommandlineクラスのインスタンスが相当します。
NcbimakeblastdbCommandlineクラスのイニシャライザに、makeblastdbコマンドのオプションに相当する引数を指定してインスタンスを作成し、そのインスタンスの__call__メソッドを呼び出して実行する(=インスタンスを関数のように扱って実行する)ことで、makeblastdbコマンドを実行させることができます。
コマンド例１
例えば次のコマンドは
makeblastdb -in ファイル名 -dbtype nucl
Biopythonを用いると、以下のようになります。
まず、makeblastdbコマンドを表すNcbimakeblastdbCommandlineクラスのインスタンス作成します。なお、Linuxコマンドのオプションが「-オプション 〇〇〇」であれば、Biopythonでは「オプション=〇〇〇」として指定すればよいのですが、「-inオプション」だけはinがPythonの予約語にあたってしまっているために、例外的に「input_file=」という形にする必要があります。
ここで作成したインスタンスをprint関数で出力してみると、実行しようとしていたLinuxコマンドが生成されていることが分かります。
from Bio.Blast import Applications
cline = Applications.NcbimakeblastdbCommandline(input_file='ファイル名', dbtype='nucl')
print(cline)
makeblastdb -dbtype nucl -in ファイル名
このインスタンスclineに()をつけて実行することで、__call__メソッドを呼び出してこのコマンドを実行します。デフォルトでは標準出力と標準エラーを返すので、それを受け取り標準出力を表示してみましょう。
stdout, stderr = cline()
print(stdout)
これでBiopythonを用いてmakeblastdbコマンドを実行することができました。
コマンド例２
ヒトゲノム配列をBLAST検索用のDBとして作成してみます。
まずは、EnsemblのFTPサイトからヒトゲノム配列を取得して解凍しましょう。なお、ここで取得するファイルは1GB程度のサイズがありますので、ダウンロードには時間がかかる場合があります。このファイルを解凍すると62.5GB程度の容量となります。
import subprocess
subprocess.run(['curl', '-O', 'ftp://ftp.ensembl.org/pub/current_fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.toplevel.fa.gz'])
subprocess.run(['gunzip', 'Homo_sapiens.GRCh38.dna.toplevel.fa.gz'])
次にmakeblastdbを実行してみますが、今度は-parse_seqids, -hash_indexのオプションをつけてみようと思います。このオプションをつけることで、DBゲノム配列から任意の領域を切り抜くことができるようになります。
ただし、-parse_seqidsオプションや-hash_indexオプションは引数の指定がありません。このように引数の指定のないオプションの場合は、Biopythonでは「オプション=True」というようにして指定します。
では、Homo_sapiens.GRCh38.dna.toplevel.faというファイルができているはずなので、Pythonでmakeblastdbコマンドを作成してみましょう。
from Bio.Blast import Applications
cline = Applications.NcbimakeblastdbCommandline(input_file='Homo_sapiens.GRCh38.dna.toplevel.fa', dbtype='nucl' ,hash_index=True ,parse_seqids=True)
print(cline)
makeblastdb -dbtype nucl -in Homo_sapiens.GRCh38.dna.toplevel.fa -parse_seqids -hash_index
これで目的のコマンドが作成できましたので、clineを関数のように扱って__call__メソッドを呼び出して、コマンドを実行しましょう。コマンドの標準出力をstdout変数で受け取って出力します。
stdout, stderr = cline()
print(stdout)
Building a new DB, current time: 06/06/2020 02:28:31 New DB name: /mnt/d/Python-NGS/Homo_sapiens.GRCh38.dna.toplevel.fa New DB title: Homo_sapiens.GRCh38.dna.toplevel.fa Sequence type: Nucleotide Keep MBits: T Maximum file size: 1000000000B Adding sequences from FASTA; added 639 sequences in 1630.91 seconds.
これでヒトゲノム配列をDBとして作成することができました。
BLAST検索の実行
先ほど作成したDBをもとにここではblastnを実行してみます。Linuxコマンドでは以下のオプションが基本となります。
blastn -query クエリファイル名 -db DBファイル名 -out 出力ファイル名
blastnコマンドは、BiopythonではNcbiblastnCommandlineクラスのインスタンスに相当します。まずは、NcbiblastnCommandlineクラスのインスタンスとしてコマンドを作成します。
from Bio.Blast import Applications
cline = Applications.NcbiblastnCommandline(query='クエリファイル名', db='DBファイル名', out='出力ファイル名')
このインスタンスの__call__メソッドを呼び出すことで、作成したコマンドを実行します。
stdout, stderr = cline()
NcbiblastnCommandlineクラスの使い方もNcbimakeblastdbCommandlineクラスの使い方と同様で、blastnコマンドの「-オプション 〇〇〇」をNcbiblastnCommandlineでは「オプション=〇〇〇」として指定します。
コマンド例
「データベースにおいて、与えられた条件を満たすデータの検索 (Windows環境でのゲノムデータの取得)」で検索したID「1819861774」の配列を取得し、ヒトゲノム配列をDBとしてblastn検索します。
from Bio import TogoWS, SeqIO
with TogoWS.entry('nucleotide', '1819861774', 'fasta') as handle:
record = SeqIO.read(handle, "fasta")
SeqIO.write(record, '1819861774.fasta', 'fasta')
from Bio.Blast import Applications
cline = Applications.NcbiblastnCommandline(query='1819861774.fasta', db='Homo_sapiens.GRCh38.dna.toplevel.fa', out='blastn-out.txt')
stdout, stderr = cline()
これで以下のファイルが出力され、配列相同部分のアラインメントを可視化することができます。
BLASTN 2.9.0+
......(略)......
>19 dna:chromosome chromosome:GRCh38:19:1:58617616:1 REF
Length=58617616
Score = 7694 bits (4166), Expect = 0.0
Identities = 4168/4169 (99%), Gaps = 0/4169 (0%)
Strand=Plus/Plus
Query 1851 GTCTCAGTCCTGGTTCCTCTGATAAGATCGACAAAAGAAACAACAAAATGAGAAGAAGAG 1910
|| |||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct 13094628 GTTTCAGTCCTGGTTCCTCTGATAAGATCGACAAAAGAAACAACAAAATGAGAAGAAGAG 13094687
Query 1911 GTTCCTCGAAAGGGGGGAGAAGAAATTTTGAGAATGGAAAAATCCCCCAGCCCAGCCCAG 1970
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct 13094688 GTTCCTCGAAAGGGGGGAGAAGAAATTTTGAGAATGGAAAAATCCCCCAGCCCAGCCCAG 13094747
......(略)......
出力フォーマットの指定 (-outfmt)
blastnコマンドの出力フォーマットは-outfmtで指定します。-outfmtの主なオプション値は次のようになります。
-outfmt 0 : アラインメントの表示
-outfmt 0 は -outfmt の指定なしと同じ出力となります。Biopythonでは次のように指定します。
cline = Applications.NcbiblastnCommandline(query='1819861774.fasta', db='Homo_sapiens.GRCh38.dna.toplevel.fa', out='blastn-out0.txt', outfmt=0)
-outfmt 5 : XMLファイルで出力
-outfmt 5 はXML形式での出力となります。Biopythonでは次のように指定します。
cline = Applications.NcbiblastnCommandline(query='1819861774.fasta', db='Homo_sapiens.GRCh38.dna.toplevel.fa', out='blastn-out5.xml', outfmt=5)
Biopythonにはここで出力したXMLファイルを解析する関数も用意されています。
-outfmt 6 : 配列相同部分のみタブ区切りで出力
-outfmt 6 は配列相同部分のみをタブ区切りで出力します。Biopythonでは次のように指定します。
cline = Applications.NcbiblastnCommandline(query='1819861774.fasta', db='Homo_sapiens.GRCh38.dna.toplevel.fa', out='blastn-out6.txt', outfmt=6)
NM_001365902.3 19 99.976 4169 1 0 1851 6019 13094628 13098796 0.0 7694 NM_001365902.3 19 100.000 535 0 0 388 922 13025018 13025552 0.0 989 NM_001365902.3 19 100.000 391 0 0 1 391 12995475 12995865 0.0 723 NM_001365902.3 19 87.954 523 59 3 388 908 3381709 3382229 4.81e-171 614 NM_001365902.3 19 100.000 179 0 0 1439 1617 13081677 13081855 5.43e-86 331 NM_001365902.3 19 100.000 151 0 0 1615 1765 13087986 13088136 1.99e-70 279 NM_001365902.3 19 99.324 148 0 1 1172 1319 13075526 13075672 1.55e-66 267 NM_001365902.3 19 99.237 131 0 1 1316 1445 13078610 13078740 4.38e-57 235 NM_001365902.3 19 98.438 128 1 1 1054 1181 13073900 13074026 9.47e-54 224 NM_001365902.3 19 100.000 92 0 0 1766 1857 13090299 13090390 1.25e-37 171 NM_001365902.3 19 100.000 76 0 0 985 1060 13073421 13073496 9.81e-29 141 NM_001365902.3 19 100.000 65 0 0 922 986 13073046 13073110 1.28e-22 121
-outfmt 7 : 配列相同部分のタブ区切り＋コメント部分を出力
-outfmt 7 もタブ区切りですが、クエリ・DBに関するコメント入りのファイルを出力します。Biopythonでは次のように指定します。
cline = Applications.NcbiblastnCommandline(query='1819861774.fasta', db='Homo_sapiens.GRCh38.dna.toplevel.fa', out='blastn-out7.txt', outfmt=7)
# BLASTN 2.9.0+ # Query: NM_001365902.3 Homo sapiens nuclear factor I X (NFIX), transcript variant 4, mRNA # Database: Homo_sapiens.GRCh38.dna.toplevel.fa # Fields: query acc.ver, subject acc.ver, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score # 12 hits found NM_001365902.3 19 99.976 4169 1 0 1851 6019 13094628 13098796 0.0 7694 NM_001365902.3 19 100.000 535 0 0 388 922 13025018 13025552 0.0 989 NM_001365902.3 19 100.000 391 0 0 1 391 12995475 12995865 0.0 723 NM_001365902.3 19 87.954 523 59 3 388 908 3381709 3382229 4.81e-171 614 NM_001365902.3 19 100.000 179 0 0 1439 1617 13081677 13081855 5.43e-86 331 NM_001365902.3 19 100.000 151 0 0 1615 1765 13087986 13088136 1.99e-70 279 NM_001365902.3 19 99.324 148 0 1 1172 1319 13075526 13075672 1.55e-66 267 NM_001365902.3 19 99.237 131 0 1 1316 1445 13078610 13078740 4.38e-57 235 NM_001365902.3 19 98.438 128 1 1 1054 1181 13073900 13074026 9.47e-54 224 NM_001365902.3 19 100.000 92 0 0 1766 1857 13090299 13090390 1.25e-37 171 NM_001365902.3 19 100.000 76 0 0 985 1060 13073421 13073496 9.81e-29 141 NM_001365902.3 19 100.000 65 0 0 922 986 13073046 13073110 1.28e-22 121 # BLAST processed 1 queries
E-valueのカットオフ値の指定 (-evalue)
-evalueオプションを指定することで、低い類似性のものをカットできます。例えば、-evalue 1e-10とするとE-valueが\(1×10^{-2}\)以下の有意性のあったものだけが出力されます。Biopythonでは次のように指定します。
cline = Applications.NcbiblastnCommandline(query='1819861774.fasta', db='Homo_sapiens.GRCh38.dna.toplevel.fa', out='blastn-out.txt', evalue=1e-10)
並列計算のコア数の指定 (-num_threads)
-num_threadsオプションを指定することで、並列計算する際のコア数を指定します。PythonでCPUのコア数を取得し、-num_threadsに指定するプログラムは次のようになります。
import multiprocessing
core = multiprocessing.cpu_count()
from Bio.Blast import Applications
cline = Applications.NcbiblastnCommandline(query='1819861774.fasta', db='Homo_sapiens.GRCh38.dna.toplevel.fa', out='blastn-out.txt', num_threads=core)
BLAST検索は時間がかかるので、並列計算を有効に使いましょう。
BLASTコマンドとBiopythonの関係
ここまでblastnについてみてきましたが、ほかのBLASTコマンドはBiopythonでは次のように対応します。
| BLASTコマンド | Biopython |
|---|---|
| blastn | NcbiblastnCommandline クラス |
| blastp | NcbiblastpCommandline クラス |
| blastx | NcbiblastxCommandline クラス |
| tblastn | NcbitblastnCommandline クラス |
| tblastx | NcbitblastxCommandline クラス |
基本的なコマンドの使い方はいずれも同様です。
BLAST結果の解析
BLAST結果の解析は以下の記事をご覧ください。
コメント