PythonでBLASTする！（ローカルBLASTの実行）【Python】

配列相同性検索のBLASTをPythonを用いて実行する方法を解説します。ここでは、特にローカルBLASTをPythonを用いて実行する方法の説明です。

バイオインフォマティクス環境

以下の方法で構築した環境であることを前提として説明します。

Windows10 + WSL (WindowsユーザーのためのPythonを用いたゲノム解析環境)

※ 検証はしていませんが、プログラム自体は基本的にmacOS、Linuxでも動作すると思われます。

ここで使用したアプリケーションのバージョンは以下の通りです。

Python 3.7.6
Biopython 1.76　　→ Biopythonについて
BLAST 2.9.0

Contents

1 PythonからローカルBLASTを実行する方法
2 DBの作成
- 2.1 コマンド例１
- 2.2 コマンド例２
3 BLAST検索の実行
4 BLAST結果の解析
5 スポンサーリンク

PythonからローカルBLASTを実行する方法

BiopythonのBio.Blast.Applicationsモジュールのクラスを用いることで、PythonからローカルBLASTを実行することができます。なお、あらかじめBLASTをPCにインストールしてパスを通しておく必要があります。もちろん、Anaconda (miniconda)を用いてインストールすればパスも通してくれるので、インストールするだけでOKです。

参考：インストール (PythonでBLASTする！（基礎編：BLASTとは？）)

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'])

参考：FTPでダウンロードする (Windows環境でのゲノムデータの取得)

次に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結果の解析は以下の記事をご覧ください。

PythonでBLASTする！（BLAST結果のXML解析）