BLAST代表基本本地对齐搜索工具.它找到了生物序列之间的相似区域. Biopython提供Bio.Blast模块来处理NCBI BLAST操作.您可以在本地连接或通过Internet连接运行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文件，如下所示 :
>>>来自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足以满足基本和高级目的.当然，有时您可能需要在本地安装它.
考虑您在线进行频繁搜索，这可能需要大量时间和高网络量，如果您有专有序列数据或建议使用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 指的是反对搜索的数据库; 查询是要匹配的序列， 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)，因此输出文件将为保存在当前工作目录中.
通常，BLAST输出使用NCBIXML模块解析为XML格式.为此，我们需要导入以下模块 :
>>> from Bio.Blast import NCBIXML
现在，使用python open方法直接打开文件和使用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)
