Biopython-overview-of-blast

提供:Dev Guides
移動先:案内検索

Biopython-BLASTの概要

BLASTは、 Basic Local Alignment Search Tool の略です。 生物学的配列間の類似領域を見つけます。 Biopythonは、NCBI BLAST操作を処理するBio.Blastモジュールを提供します。 BLASTは、ローカル接続またはインターネット接続で実行できます。

次のセクションでこれらの2つの接続を簡単に理解しましょう-

インターネット経由で実行する

Biopythonは、BLASTのオンラインバージョンを呼び出すBio.Blast.NCBIWWWモジュールを提供します。 これを行うには、次のモジュールをインポートする必要があります-

>>> from Bio.Blast import NCBIWWW

NCBIWWモ​​ジュールは、BLASTオンラインバージョンhttps://blast.ncbi.nlm.nih.gov/Blast.cgiを照会するqblast関数を提供します。 qblastは、オンラインバージョンでサポートされているすべてのパラメーターをサポートしています。

このモジュールに関するヘルプを取得するには、以下のコマンドを使用し、機能を理解します-

>>> help(NCBIWWW.qblast)
Help on function qblast in module Bio.Blast.NCBIWWW:
qblast(
   program, database, sequence,
   url_base = 'https://blast.ncbi.nlm.nih.gov/Blast.cgi',
   auto_format = None,
   composition_based_statistics = None,
   db_genetic_code =  None,
   endpoints = None,
   entrez_query = '(none)',
   expect = 10.0,
   filter = None,
   gapcosts = None,
   genetic_code = None,
   hitlist_size = 50,
   i_thresh = None,
   layout = None,
   lcase_mask = None,
   matrix_name = None,
   nucl_penalty = None,
   nucl_reward = None,
   other_advanced = None,
   perc_ident = None,
   phi_pattern = None,
   query_file = None,
   query_believe_defline = None,
   query_from = None,
   query_to = None,
   searchsp_eff = None,
   service = None,
   threshold = None,
   ungapped_alignment = None,
   word_size = None,
   alignments = 500,
   alignment_view = None,
   descriptions = 500,
   entrez_links_new_window = None,
   expect_low = None,
   expect_high = None,
   format_entrez_query = None,
   format_object = None,
   format_type = 'XML',
   ncbi_gi = None,
   results_file = None,
   show_overview = None,
   megablast = None,
   template_type = None,
   template_length = None
)

   BLAST search using NCBI's QBLAST server or a cloud service provider.

   Supports all parameters of the qblast API for Put and Get.

   Please note that BLAST on the cloud supports the NCBI-BLAST Common
   URL API (http://ncbi.github.io/blast-cloud/dev/apil).
   To use this feature, please set url_base to 'http://host.my.cloud.service.provider.com/cgi-bin/blast.cgi' and
   format_object = 'Alignment'. For more details, please see 8. Biopython – Overview of BLAST

https://blast.ncbi.nlm.nih.gov/Blast.cgi?PAGE_TYPE = BlastDocs&DOC_TYPE = CloudBlast

   Some useful parameters:

   - program blastn, blastp, blastx, tblastn, or tblastx (lower case)
   - database Which database to search against (e.g. "nr").
   - sequence The sequence to search.
   - ncbi_gi TRUE/FALSE whether to give 'gi' identifier.
   - descriptions Number of descriptions to show. Def 500.
   - alignments Number of alignments to show. Def 500.
   - expect An expect value cutoff. Def 10.0.
   - matrix_name Specify an alt. matrix (PAM30, PAM70, BLOSUM80, BLOSUM45).
   - filter "none" turns off filtering. Default no filtering
   - format_type "HTML", "Text", "ASN.1", or "XML". Def. "XML".
   - entrez_query Entrez query to limit Blast search
   - hitlist_size Number of hits to return. Default 50
   - megablast TRUE/FALSE whether to use MEga BLAST algorithm (blastn only)
   - service plain, psi, phi, rpsblast, megablast (lower case)

   This function does no checking of the validity of the parameters
   and passes the values to the server as is. More help is available at:
   https://ncbi.github.io/blast-cloud/dev/apil

通常、qblast関数の引数は、BLAST Webページで設定できるさまざまなパラメーターに基本的に類似しています。 これにより、qblast関数の理解が容易になり、使用するための学習曲線が縮小されます。

接続と検索

BLASTオンラインバージョンの接続と検索のプロセスを理解するために、Biopythonを介してオンラインBLASTサーバーに対して単純なシーケンス検索(ローカルシーケンスファイルで利用可能)を行いましょう。

ステップ1 *-Biopythonディレクトリに *blast_example.fasta という名前のファイルを作成し、入力として以下のシーケンス情報を提供します

Example of a single sequence in FASTA/Pearson format:
>sequence A ggtaagtcctctagtacaaacacccccaatattgtgatataattaaaattatattcatat
tctgttgccagaaaaaacacttttaggctatattagagccatcttctttgaagcgttgtc

>sequence B ggtaagtcctctagtacaaacacccccaatattgtgatataattaaaattatattca
tattctgttgccagaaaaaacacttttaggctatattagagccatcttctttgaagcgttgtc
  • ステップ2 *-NCBIWWWモジュールをインポートします。
>>> from Bio.Blast import NCBIWWW

ステップ3 *-python IOモジュールを使用して、シーケンスファイル *blast_example.fasta を開きます。

>>> sequence_data = open("blast_example.fasta").read()
>>> sequence_data
'Example of a single sequence in FASTA/Pearson format:\n\n\n> sequence
A\nggtaagtcctctagtacaaacacccccaatattgtgatataattaaaatt
atattcatat\ntctgttgccagaaaaaacacttttaggctatattagagccatcttctttg aagcgttgtc\n\n'
  • ステップ4 *-ここで、メインデータとしてシーケンスデータを渡すqblast関数を呼び出します。 もう1つのパラメーターは、データベース(nt)と内部プログラム(blastn)を表します。
>>> result_handle = NCBIWWW.qblast("blastn", "nt", sequence_data)
>>> result_handle
<_io.StringIO object at 0x000001EC9FAA4558>
*blast_results* は検索結果を保持します。 後で使用するためにファイルに保存したり、解析して詳細を取得したりできます。 次のセクションでその方法を学習します。
  • ステップ5 *-以下に示すように、fastaファイル全体を使用するのではなく、Seqオブジェクトを使用しても同じ機能を実行できます-
>>> from Bio import SeqIO
>>> seq_record = next(SeqIO.parse(open('blast_example.fasta'),'fasta'))
>>> seq_record.id
'sequence'
>>> seq_record.seq
Seq('ggtaagtcctctagtacaaacacccccaatattgtgatataattaaaattatat...gtc',
SingleLetterAlphabet())

ここで、メインパラメータとしてSeqオブジェクトであるrecord.seqを渡すqblast関数を呼び出します。

>>> result_handle = NCBIWWW.qblast("blastn", "nt", seq_record.seq)
>>> print(result_handle)
<_io.StringIO object at 0x000001EC9FAA4558>

BLASTは、シーケンスに識別子を自動的に割り当てます。

  • ステップ6 *-result_handleオブジェクトには結果全体が含まれ、後で使用するためにファイルに保存できます。
>>> with open('results.xml', 'w') as save_file:
>>>   blast_results = result_handle.read()
>>>   save_file.write(blast_results)

結果ファイルの解析方法については、後のセクションで説明します。

スタンドアロンBLASTの実行

このセクションでは、ローカルシステムでBLASTを実行する方法について説明します。 ローカルシステムでBLASTを実行すると、高速になり、独自のデータベースを作成してシーケンスを検索することもできます。

BLASTの接続

一般に、BLASTをローカルで実行することは、サイズが大きく、ソフトウェアを実行するために余分な労力が必要で、コストがかかるため、お勧めしません。 オンラインBLASTは、基本的および高度な目的には十分です。 もちろん、ローカルにインストールする必要がある場合があります。

多くの時間と大量のネットワークを必要とするオンラインで頻繁に検索を行っていることを考慮してください。独自のシーケンスデータまたはIP関連の問題がある場合は、ローカルにインストールすることをお勧めします。

これを行うには、以下の手順に従う必要があります-

BLASTソフトウェアは、サイトに多くのデータベースを提供します。 ブラストデータベースサイトからftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/alu.n.gz[alu.n.gz]ファイルをダウンロードし、aluフォルダーに解凍します。 このファイルはFASTA形式です。 ブラストアプリケーションでこのファイルを使用するには、まずファイルをFASTA形式からブラストデータベース形式に変換する必要があります。 BLASTは、この変換を行うmakeblastdbアプリケーションを提供します。

以下のコードスニペットを使用してください-

cd/path/to/alu
makeblastdb -in alu.n -parse_seqids -dbtype nucl -out alun

上記のコードを実行すると、入力ファイルalu.nが解析され、複数のファイルalun.nsq、alun.nsiなどとしてBLASTデータベースが作成されます。 これで、このデータベースにクエリを実行してシーケンスを見つけることができます。

ローカルサーバーにBLASTをインストールし、サンプルのBLASTデータベース、 alun を照会します。

  • ステップ3 *-データベースを照会するサンプルシーケンスファイルを作成しましょう。 ファイルsearch.fsaを作成し、以下のデータをそこに入れます。
>gnl|alu|Z15030_HSAL001056 (Alu-J)
AGGCTGGCACTGTGGCTCATGCTGAAATCCCAGCACGGCGGAGGACGGCGGAAGATTGCT
TGAGCCTAGGAGTTTGCGACCAGCCTGGGTGACATAGGGAGATGCCTGTCTCTACGCAAA
AGAAAAAAAAAATAGCTCTGCTGGTGGTGCATGCCTATAGTCTCAGCTATCAGGAGGCTG
GGACAGGAGGATCACTTGGGCCCGGGAGTTGAGGCTGTGGTGAGCCACGATCACACCACT
GCACTCCAGCCTGGGTGACAGAGCAAGACCCTGTCTCAAAACAAACAAATAA
>gnl|alu|D00596_HSAL003180 (Alu-Sx)
AGCCAGGTGTGGTGGCTCACGCCTGTAATCCCACCGCTTTGGGAGGCTGAGTCAGATCAC
CTGAGGTTAGGAATTTGGGACCAGCCTGGCCAACATGGCGACACCCCAGTCTCTACTAAT
AACACAAAAAATTAGCCAGGTGTGCTGGTGCATGTCTGTAATCCCAGCTACTCAGGAGGC
TGAGGCATGAGAATTGCTCACGAGGCGGAGGTTGTAGTGAGCTGAGATCGTGGCACTGTA
CTCCAGCCTGGCGACAGAGGGAGAACCCATGTCAAAAACAAAAAAAGACACCACCAAAGG
TCAAAGCATA
>gnl|alu|X55502_HSAL000745 (Alu-J)
TGCCTTCCCCATCTGTAATTCTGGCACTTGGGGAGTCCAAGGCAGGATGATCACTTATGC
CCAAGGAATTTGAGTACCAAGCCTGGGCAATATAACAAGGCCCTGTTTCTACAAAAACTT
TAAACAATTAGCCAGGTGTGGTGGTGCGTGCCTGTGTCCAGCTACTCAGGAAGCTGAGGC
AAGAGCTTGAGGCTACAGTGAGCTGTGTTCCACCATGGTGCTCCAGCCTGGGTGACAGGG
CAAGACCCTGTCAAAAGAAAGGAAGAAAGAACGGAAGGAAAGAAGGAAAGAAACAAGGAG
AG

シーケンスデータはalu.nファイルから収集されます。したがって、データベースと一致します。

ステップ4 *-BLASTソフトウェアは、データベースを検索するための多くのアプリケーションを提供し、blastnを使用します。 *blastnアプリケーションには、db、query、outの少なくとも3つの引数が必要です。 db は、検索対象のデータベースを指します。 query は一致するシーケンスで、 out は結果を保存するファイルです。 今、この単純なクエリを実行するために以下のコマンドを実行します-

blastn -db alun -query search.fsa -out results.xml -outfmt 5

上記のコマンドを実行すると、以下に示すように results.xml ファイルで検索され、出力が得られます(部分的にデータ)

<?xml version = "1.0"?>
<!DOCTYPE BlastOutput PUBLIC "-//NCBI//NCBI BlastOutput/EN"
   "http://www.ncbi.nlm.nih.gov/dtd/NCBI_BlastOutput.dtd">
<BlastOutput>
   <BlastOutput_program>blastn</BlastOutput_program>
   <BlastOutput_version>BLASTN 2.7.1+</BlastOutput_version>
   <BlastOutput_reference>Zheng Zhang, Scott Schwartz, Lukas Wagner, and Webb
      Miller (2000), "A greedy algorithm for aligning DNA sequences", J
      Comput Biol 2000; 7(1-2):203-14.
   </BlastOutput_reference>

   <BlastOutput_db>alun</BlastOutput_db>
   <BlastOutput_query-ID>Query_1</BlastOutput_query-ID>
   <BlastOutput_query-def>gnl|alu|Z15030_HSAL001056 (Alu-J)</BlastOutput_query-def>
   <BlastOutput_query-len>292</BlastOutput_query-len>
   <BlastOutput_param>
      <Parameters>
         <Parameters_expect>10</Parameters_expect>
         <Parameters_sc-match>1</Parameters_sc-match>
         <Parameters_sc-mismatch>-2</Parameters_sc-mismatch>
         <Parameters_gap-open>0</Parameters_gap-open>
         <Parameters_gap-extend>0</Parameters_gap-extend>
         <Parameters_filter>L;m;</Parameters_filter>
      </Parameters>
   </BlastOutput_param>
   <BlastOutput_iterations>
      <Iteration>
         <Iteration_iter-num>1</Iteration_iter-num><Iteration_query-ID>Query_1</Iteration_query-ID>
         <Iteration_query-def>gnl|alu|Z15030_HSAL001056 (Alu-J)</Iteration_query-def>
         <Iteration_query-len>292</Iteration_query-len>
         <Iteration_hits>
         <Hit>
            <Hit_num>1</Hit_num>
            <Hit_id>gnl|alu|Z15030_HSAL001056</Hit_id>
            <Hit_def>(Alu-J)</Hit_def>
            <Hit_accession>Z15030_HSAL001056</Hit_accession>
            <Hit_len>292</Hit_len>
            <Hit_hsps>
               <Hsp>
                 <Hsp_num>1</Hsp_num>
                  <Hsp_bit-score>540.342</Hsp_bit-score>
                  <Hsp_score>292</Hsp_score>
                  <Hsp_evalue>4.55414e-156</Hsp_evalue>
                  <Hsp_query-from>1</Hsp_query-from>
                  <Hsp_query-to>292</Hsp_query-to>
                  <Hsp_hit-from>1</Hsp_hit-from>
                  <Hsp_hit-to>292</Hsp_hit-to>
                  <Hsp_query-frame>1</Hsp_query-frame>
                  <Hsp_hit-frame>1</Hsp_hit-frame>
                  <Hsp_identity>292</Hsp_identity>
                  <Hsp_positive>292</Hsp_positive>
                  <Hsp_gaps>0</Hsp_gaps>
                  <Hsp_align-len>292</Hsp_align-len>

                  <Hsp_qseq>
                     AGGCTGGCACTGTGGCTCATGCTGAAATCCCAGCACGGCGGAGGACGGCGGAAGATTGCTTGAGCCTAGGAGTTTG
                     CGACCAGCCTGGGTGACATAGGGAGATGCCTGTCTCTACGCAAAAGAAAAAAAAAATAGCTCTGCTGGTGGTGCATG
                     CCTATAGTCTCAGCTATCAGGAGGCTGGGACAGGAGGATCACTTGGGCCCGGGAGTTGAGGCTGTGGTGAGCC
                     ACGATCACACCACTGCACTCCAGCCTGGGTGACAGAGCAAGACCCTGTCTCAAAACAAACAAATAA
                  </Hsp_qseq>

                  <Hsp_hseq>
                     AGGCTGGCACTGTGGCTCATGCTGAAATCCCAGCACGGCGGAGGACGGCGGAAGATTGCTTGAGCCTAGGA
                     GTTTGCGACCAGCCTGGGTGACATAGGGAGATGCCTGTCTCTACGCAAAAGAAAAAAAAAATAGCTCTGCT
                     GGTGGTGCATGCCTATAGTCTCAGCTATCAGGAGGCTGGGACAGGAGGATCACTTGGGCCCGGGAGTTGAGG
                     CTGTGGTGAGCCACGATCACACCACTGCACTCCAGCCTGGGTGACAGAGCAAGACCCTGTCTCAAAACAAAC
                     AAATAA
                  </Hsp_hseq>

                  <Hsp_midline>
                     |||||||||||||||||||||||||||||||||||||||||||||||||||||
                     |||||||||||||||||||||||||||||||||||||||||||||||||||||
                     |||||||||||||||||||||||||||||||||||||||||||||||||||||
                     |||||||||||||||||||||||||||||||||||||||||||||||||||||
                     |||||||||||||||||||||||||||||||||||||||||||||||||||||
                     |||||||||||||||||||||||||||
                  </Hsp_midline>
               </Hsp>
            </Hit_hsps>
         </Hit>
         .........................
         .........................
         .........................
         </Iteration_hits>
         <Iteration_stat>
            <Statistics>
               <Statistics_db-num>327</Statistics_db-num>
               <Statistics_db-len>80506</Statistics_db-len>
               <Statistics_hsp-lenv16</Statistics_hsp-len>
               <Statistics_eff-space>21528364</Statistics_eff-space>
               <Statistics_kappa>0.46</Statistics_kappa>
               <Statistics_lambda>1.28</Statistics_lambda>
               <Statistics_entropy>0.85</Statistics_entropy>
            </Statistics>
         </Iteration_stat>
      </Iteration>
   </BlastOutput_iterations>
</BlastOutput>

上記のコマンドは、以下のコードを使用してPython内で実行できます-

>>> from Bio.Blast.Applications import NcbiblastnCommandline
>>> blastn_cline = NcbiblastnCommandline(query = "search.fasta", db = "alun",
outfmt = 5, out = "results.xml")
>>> stdout, stderr = blastn_cline()

ここで、最初のものは爆風出力へのハンドルであり、2番目は爆風コマンドによって生成される可能性のあるエラー出力です。

出力ファイルをコマンドライン引数(out =“ results.xml”)として提供し、出力形式をXML(outfmt = 5)として設定しているため、出力ファイルは現在の作業ディレクトリに保存されます。

BLAST結果の解析

一般に、BLAST出力はNCBIXMLモジュールを使用してXML形式として解析されます。 これを行うには、次のモジュールをインポートする必要があります-

>>> from Bio.Blast import NCBIXML

ここで、 Python open method を使用してファイルを直接開き、以下に示すように use NCBIXML parse method を使用します-

>>> E_VALUE_THRESH = 1e-20
>>> for record in NCBIXML.parse(open("results.xml")):
>>>     if record.alignments:
>>>        print("\n")
>>>        print("query: %s" % record.query[:100])
>>>        for align in record.alignments:
>>>           for hsp in align.hsps:
>>>              if hsp.expect < E_VALUE_THRESH:
>>>                 print("match: %s " % align.title[:100])

これにより、次のような出力が生成されます-

query: gnl|alu|Z15030_HSAL001056 (Alu-J)
match: gnl|alu|Z15030_HSAL001056 (Alu-J)
match: gnl|alu|L12964_HSAL003860 (Alu-J)
match: gnl|alu|L13042_HSAL003863 (Alu-FLA?)
match: gnl|alu|M86249_HSAL001462 (Alu-FLA?)
match: gnl|alu|M29484_HSAL002265 (Alu-J)

query: gnl|alu|D00596_HSAL003180 (Alu-Sx)
match: gnl|alu|D00596_HSAL003180 (Alu-Sx)
match: gnl|alu|J03071_HSAL001860 (Alu-J)
match: gnl|alu|X72409_HSAL005025 (Alu-Sx)

query: gnl|alu|X55502_HSAL000745 (Alu-J)
match: gnl|alu|X55502_HSAL000745 (Alu-J)