Biopython 简明教程

Biopython - Overview of BLAST

BLAST 代表 Basic Local Alignment Search Tool 。它在生物序列之间找到相似区域。Biopython 提供了 Bio.Blast 模块来处理 NCBI BLAST 操作。你可以在本地连接或互联网连接上运行 BLAST。

让我们在以下部分简要了解这两个连接 -

Running over Internet

Biopython 提供 Bio.Blast.NCBIWWW 模块来调用 BLAST 的在线版本。要做到这一点，我们需要导入以下模块 -

>>> from Bio.Blast import NCBIWWW

NCBIWW 模块提供 qblast 函数来查询 BLAST 在线版本， https://blast.ncbi.nlm.nih.gov/Blast.cgi 。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/api.html).
   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/api.html

通常，qblast 函数的参数基本上类似于你可以在 BLAST 网页上设置的不同参数。这使得 qblast 函数易于理解，并且减少了使用它的学习曲线。

Connecting and Searching

为了理解连接和搜索 BLAST 在线版本的过程，让我们通过 Biopython 对本地序列文件中的简单序列搜索进行在线 BLAST 服务器搜索。

Step 1 − 在 Biopython 目录中创建一个名为 blast_example.fasta 的文件，并输入以下序列信息作为输入

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

>sequence B ggtaagtcctctagtacaaacacccccaatattgtgatataattaaaattatattca
tattctgttgccagaaaaaacacttttaggctatattagagccatcttctttgaagcgttgtc

Step 2 − 导入 NCBIWWW 模块。

>>> from Bio.Blast import NCBIWWW

Step 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'

Step 4 − 现在，调用 qblast 函数，将序列数据作为主要参数传递。其他参数表示数据库 (nt) 和内部程序 (blastn)。

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

blast_results 保存了我们的搜索结果。可以将其存储到文件中以便稍后使用，还可以进行解析以获取详细信息。我们将在下一部分学习如何执行此操作。

Step 5 − 同样，也可以使用 Seq 对象执行相同的操作，而不是使用整个 fasta 文件，如下所示 −

>>> 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())

现在，调用 qblast 函数，将 Seq 对象 record.seq 作为主参数传递。

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

BLAST 将自动为您的序列分配一个标识符。

Step 6 − result_handle 对象将拥有整个结果，可以将其保存到文件中以便稍后使用。

>>> with open('results.xml', 'w') as save_file:
>>>   blast_results = result_handle.read()
>>>   save_file.write(blast_results)

我们将在后面的部分中了解如何解析结果文件。

Running Standalone BLAST

此部分说明如何在本地系统上运行 BLAST。如果你在本地系统上运行 BLAST，可能会更快，还可以让你创建自己的数据库来搜索序列。

Connecting BLAST

通常，不建议在本地运行 BLAST，因为它体积过大，需要额外的精力来运行软件，并且涉及成本。在线 BLAST 足以满足基本和高级目的。当然，有时候你可能需要在本地安装它。

考虑到你正在进行频繁的在线搜索，这可能需要大量时间和高网络容量，如果你有专有序列数据或 IP 相关问题，那么建议在本地安装它。

要做到这一点，我们需要执行以下步骤 −

Step 1 − 使用给定的链接下载并安装最新的 blast 二进制文件 − ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/

Step 2 − 使用以下链接下载并解压最新且必要的数据库 − ftp://ftp.ncbi.nlm.nih.gov/blast/db/

BLAST 软件在他们的网站上提供了大量数据库。让我们从 blast 数据库网站下载 alu.n.gz 文件并将其解压到 alu 文件夹中。此文件采用 FASTA 格式。要在我们的 blast 应用程序中使用此文件，我们需要首先将文件从 FASTA 格式转换为 blast 数据库格式。BLAST 提供 makeblastdb 应用程序来执行此转换。

使用以下代码片段 −

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

运行上述代码将解析输入文件 alu.n 并创建 BLAST 数据库，如多个文件 alun.nsq、alun.nsi 等。现在，我们可以查询该数据库来查找序列。

我们已经在本地服务器上安装了 BLAST，并且还有示例 BLAST 数据库 alun 来对照其进行查询。

Step 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 文件收集的；因此，它与我们的数据库匹配。

Step 4 − BLAST 软件提供了许多应用程序来搜索数据库，我们使用 blastn。 blastn application requires minimum of three arguments, db, query and out. 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()

其中，第一个是 BLAST 输出的处理程序，而第二个是 BLAST 命令生成的可能的错误输出。

由于我们已提供输出文件作为命令行参数（out = “results.xml”）并将输出格式设置为 XML（outfmt = 5），因此输出文件将保存于当前的工作目录。

Parsing BLAST Result

通常，BLAST 输出使用 NCBIXML 模块解析为 XML 格式。为此，我们需要导入以下模块：

>>> from Bio.Blast import NCBIXML

现在， open the file directly using 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)