极客教程Biopython BLAST概述
BLAST是 Basic Local Alignment Search Tool的缩写 它可以发现生物序列之间的相似性区域。Biopython提供了Bio.Blast模块来处理NCBI的BLAST操作。你可以在本地连接或通过互联网连接运行BLAST。
让我们在下面的章节中简单了解一下这两种连接方式。
通过互联网运行
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函数容易理解，同时也减少了使用它的学习曲线。
连接和搜索
为了了解连接和搜索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函数，传递序列数据作为主参数。其他参数代表数据库（nt）和内部程序（blastn）。
>>> result_handle = NCBIWWW.qblast("blastn", "nt", sequence_data)
>>> result_handle
<_io.StringIO object at 0x000001EC9FAA4558>
blast_results 持有我们的搜索结果。它可以被保存到一个文件中供以后使用，也可以被解析以获得细节。我们将在下一节中学习如何做到这一点。
第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将为你的序列自动分配一个标识符。
第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相关的问题，那么建议在本地安装它。
要做到这一点，我们需要遵循以下步骤 –
第1步 - 使用给定的链接下载并安装最新的blast二进制文件 –ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/
第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 来查询它。
第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。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>
The above command can be run inside the python using the below code −
>>> from Bio.Blast.Applications import NcbiblastnCommandline
>>> blastn_cline = NcbiblastnCommandline(query = "search.fasta", db = "alun",
outfmt = 5, out = "results.xml")
>>> stdout, stderr = blastn_cline()
这里，第一个是爆炸输出的句柄，第二个是爆炸命令可能产生的错误输出。
由于我们提供了输出文件作为命令行参数（out = “results.xml”）并将输出格式设置为XML（outfmt = 5），输出文件将被保存在当前工作目录中。
解析BLAST结果
一般来说，BLAST输出是使用NCBIXML模块解析为XML格式。要做到这一点，我们需要导入以下模块 –
>>> from Bio.Blast import NCBIXML
现在， 使用python打开方法直接打开文件 ，并 使用NCBIXML解析方法 ，如下图所示
>>> 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)