BLAST (新)
嗨，每个人都喜欢 BLAST，对吧？我的意思是，天哪，还能有什么比将你的序列与已知世界中的所有其他序列进行比较更容易的事了？但是，当然，这一部分不是关于 BLAST 有多酷，因为我们已经知道这一点。而是关于 BLAST 的问题——处理大型运行产生的海量数据，以及自动化 BLAST 运行，可能非常困难。
幸运的是，Biopython 的开发人员非常了解这一点，因此他们开发了许多用于处理 BLAST 并使其变得更轻松的工具。本节详细介绍了如何使用这些工具以及如何用它们做一些有用的事情。
处理 BLAST 可以分为两个步骤，这两个步骤都可以在 Biopython 中完成。首先，为你的查询序列运行 BLAST，并获取一些输出。其次，在 Python 中解析 BLAST 输出以进行进一步分析。
你第一次接触运行 BLAST，可能是在 NCBI BLAST 网页 上。实际上，你可以通过多种方式运行 BLAST，这些方式可以分为几种类别。最重要的区别是在本地（在您自己的机器上）运行 BLAST，以及在远程（在另一台机器上，通常是 NCBI 服务器）运行 BLAST。我们将在本章开始部分介绍如何在 Python 脚本中调用 NCBI 在线 BLAST 服务。
通过互联网运行 BLAST
我们在 Bio.Blast 模块中使用 qblast 函数来调用 BLAST 的在线版本。
NCBI 指南 规定
不要比每 10 秒一次更频繁地联系服务器。
不要比每分钟一次更频繁地轮询任何单个 RID。
使用 URL 参数 email 和 tool，以便 NCBI 在出现问题时可以联系你。
如果要提交 50 次以上搜索，请在周末或工作日晚上 9 点到次日凌晨 5 点之间运行脚本。
Blast.qblast 会自动遵循前两点。为了满足第三点，请设置 Blast.email 变量（默认情况下，Blast.tool 变量已设置为 "biopython"）
>>> from Bio import Blast
>>> Blast.tool
'biopython'
>>> Blast.email = "[email protected]"
BLAST 参数
qblast 函数有三个非可选参数
第一个参数是要用于搜索的 BLAST 程序，以小写字符串形式表示。程序及其选项在 NCBI BLAST 网页 上有说明。目前，
qblast仅适用于 blastn、blastp、blastx、tblast 和 tblastx。第二个参数指定要搜索的数据库。同样，该选项在 NCBI 的 BLAST 帮助页面 上提供。
第三个参数是一个字符串，包含你的查询序列。它可以是序列本身、fasta 格式的序列，或类似 GI 号码的标识符。
qblast 函数还接受许多其他可选参数，这些参数基本上类似于你可以在 BLAST 网页上设置的不同参数。我们这里只重点介绍其中一些。
url_base参数设置用于通过互联网运行 BLAST 的基本 URL。默认情况下，它连接到 NCBI，但你可以使用它来连接到云中运行的 NCBI BLAST 实例。有关更多详细信息，请参阅qblast函数的文档。qblast函数可以以多种格式返回 BLAST 结果，你可以使用可选的format_type关键字进行选择："XML"、"HTML"、"Text"、"XML2"、"JSON2"或"Tabular"。默认值为"XML"，因为这是解析器期望的格式，在下面的第 解析 BLAST 输出 节中进行了说明。expect参数设置期望值或 e 值阈值。
有关可选 BLAST 参数的更多信息，请参阅 NCBI 自身的文档，或者 Biopython 中内置的文档。
>>> from Bio import Blast
>>> help(Blast.qblast)
请注意，NCBI BLAST 网站上的默认设置与 QBLAST 上的默认设置略有不同。如果得到不同的结果，则需要检查参数（例如，期望值阈值和间隙值）。
例如，如果你有一个核苷酸序列，你想要使用 BLASTN 对其进行搜索，并且你了解你的查询序列的 GI 号码，你可以使用
>>> from Bio import Blast
>>> result_stream = Blast.qblast("blastn", "nt", "8332116")
或者，如果我们已经将查询序列保存在 FASTA 格式的文件中，我们只需要打开文件并将此记录读入为字符串，并将其用作查询参数即可
>>> from Bio import Blast
>>> fasta_string = open("m_cold.fasta").read()
>>> result_stream = Blast.qblast("blastn", "nt", fasta_string)
我们也可以将 FASTA 文件读入为 SeqRecord，然后仅提供序列本身
>>> from Bio import Blast
>>> from Bio import SeqIO
>>> record = SeqIO.read("m_cold.fasta", "fasta")
>>> result_stream = Blast.qblast("blastn", "nt", record.seq)
仅提供序列意味着 BLAST 将自动为你的序列分配一个标识符。你可能更喜欢对 SeqRecord 对象调用 format 来创建 FASTA 字符串（其中将包含现有标识符）
>>> from Bio import Blast
>>> from Bio import SeqIO
>>> records = SeqIO.parse("ls_orchid.gbk", "genbank")
>>> record = next(records)
>>> result_stream = Blast.qblast("blastn", "nt", format(record, "fasta"))
如果你将序列保存在非 FASTA 文件格式中，你可以使用 Bio.SeqIO（参见第 序列输入/输出 章）提取这些序列，这种方法更有意义。
保存 BLAST 结果
无论你为 qblast() 函数传递什么参数，你都应该以 bytes 数据流的形式获取结果（默认情况下为 XML 格式）。下一步将是解析 XML 输出到代表搜索结果的 Python 对象（第 解析 BLAST 输出 节），但你可能希望先保存输出文件的本地副本。我发现这在调试提取 BLAST 结果中信息的代码时特别有用（因为重新运行在线搜索速度很慢，并且浪费了 NCBI 的计算机时间）。
我们需要稍微注意一下，因为我们可以使用 result_stream.read() 仅读取一次 BLAST 输出——再次调用 result_stream.read() 会返回一个空 bytes 对象。
>>> with open("my_blast.xml", "wb") as out_stream:
... out_stream.write(result_stream.read())
...
>>> result_stream.close()
执行完此操作后，结果将保存在 my_blast.xml 文件中，并且 result_stream 中的所有数据都被提取了（因此我们将其关闭了）。但是，BLAST 解析器的 parse 函数（在第 解析 BLAST 输出 节中进行了说明）接受文件类对象，因此我们可以将保存的文件打开为输入 bytes
>>> result_stream = open("my_blast.xml", "rb")
现在我们已经将 BLAST 结果重新读入数据流，我们准备对它们做点什么，因此这直接带我们进入了解析部分（参见下面的第 解析 BLAST 输出 节）。你现在可能想要跳到那里……。
以其他格式获取 BLAST 输出
通过在调用 qblast 时使用 format_type 参数，你可以以除 XML 之外的其他格式获取 BLAST 输出。以下是如何读取 JSON 格式的 BLAST 输出的示例。使用 format_type="JSON2"，Blast.qblast 提供的数据将以压缩的 JSON 格式保存
>>> from Bio import Blast
>>> from Bio import SeqIO
>>> record = SeqIO.read("m_cold.fasta", "fasta")
>>> result_stream = Blast.qblast("blastn", "nt", record.seq, format_type="JSON2")
>>> data = result_stream.read()
>>> data[:4]
b'PK\x03\x04'
这是 ZIP 文件的魔数。
>>> with open("myzipfile.zip", "wb") as out_stream:
... out_stream.write(data)
...
13813
请注意，我们将数据读入和写入为 bytes。现在打开我们创建的 ZIP 文件
>>> import zipfile
>>> myzipfile = zipfile.ZipFile("myzipfile.zip")
>>> myzipfile.namelist()
['N5KN7UMJ013.json', 'N5KN7UMJ013_1.json']
>>> stream = myzipfile.open("N5KN7UMJ013.json")
>>> data = stream.read()
这些数据是 bytes，因此我们需要对其进行解码以获取字符串对象
>>> data = data.decode()
>>> print(data)
{
"BlastJSON": [
{"File": "N5KN7UMJ013_1.json" }
]
}
现在打开 ZIP 文件中包含的第二个文件，以 JSON 格式获取 BLAST 结果
>>> stream = myzipfile.open("N5KN7UMJ013_1.json")
>>> data = stream.read()
>>> len(data)
145707
>>> data = data.decode()
>>> print(data)
{
"BlastOutput2": {
"report": {
"program": "blastn",
"version": "BLASTN 2.14.1+",
"reference": "Stephen F. Altschul, Thomas L. Madden, Alejandro A. ...
"search_target": {
"db": "nt"
},
"params": {
"expect": 10,
"sc_match": 2,
"sc_mismatch": -3,
"gap_open": 5,
"gap_extend": 2,
"filter": "L;m;"
},
"results": {
"search": {
"query_id": "Query_69183",
"query_len": 1111,
"query_masking": [
{
"from": 797,
"to": 1110
}
],
"hits": [
{
"num": 1,
"description": [
{
"id": "gi|1219041180|ref|XM_021875076.1|",
...
我们可以使用 Python 标准库中的 JSON 解析器将 JSON 数据转换为常规的 Python 字典
>>> import json
>>> d = json.loads(data)
>>> print(d)
{'BlastOutput2': {'report': {'program': 'blastn', 'version': 'BLASTN 2.14.1+',
'reference': 'Stephen F. Altschul, Thomas L. Madden, Alejandro A. Schäffer,
Jinghui Zhang, Zheng Zhang, Webb Miller, and David J. Lipman (1997),
"Gapped BLAST and PSI-BLAST: a new generation of protein database search programs",
Nucleic Acids Res. 25:3389-3402.',
'search_target': {'db': 'nt'}, 'params': {'expect': 10, 'sc_match': 2,
'sc_mismatch': -3, 'gap_open': 5, 'gap_extend': 2, 'filter': 'L;m;'},
'results': {'search': {'query_id': 'Query_128889', 'query_len': 1111,
'query_masking': [{'from': 797, 'to': 1110}], 'hits': [{'num': 1,
'description': [{'id': 'gi|1219041180|ref|XM_021875076.1|', 'accession':
'XM_021875076', 'title':
'PREDICTED: Chenopodium quinoa cold-regulated 413 plasma membrane protein 2-like (LOC110697660), mRNA',
'taxid': 63459, 'sciname': 'Chenopodium quinoa'}], 'len': 1173, 'hsps':
[{'num': 1, 'bit_score': 435.898, 'score': 482, 'evalue': 9.02832e-117,
'identity': 473, 'query_from'
...
本地运行 BLAST
简介
在本地运行 BLAST（而不是通过互联网，参见第 通过互联网运行 BLAST 节）至少有两个主要优势
本地 BLAST 可能比通过互联网运行 BLAST 更快；
本地 BLAST 允许你创建自己的数据库来搜索序列。
在本地运行 BLAST 的另一个原因可能是处理专有或未公开的序列数据。您可能不被允许重新发布这些序列，因此将它们提交给 NCBI 作为 BLAST 查询将不可行。
不幸的是，也有一些主要的缺点 - 安装所有组件并使其正确设置需要一些努力。
本地 BLAST 需要安装命令行工具。
本地 BLAST 需要设置（并可能保持更新）大型 BLAST 数据库。
独立 NCBI BLAST+
“新的”NCBI BLAST+ 套件于 2009 年发布。它取代了旧的 NCBI “传统” BLAST 包（见 其他版本的 BLAST）。
本节将简要介绍如何从 Python 中使用这些工具。如果您已经阅读或尝试过第 对齐工具 节中的对齐工具示例，那么这一切都应该看起来相当简单。首先，我们构建一个命令行字符串（就像您在运行独立 BLAST 时手动在命令行提示符中输入的那样）。然后，我们可以在 Python 中执行此命令。
例如，如果您有一个包含基因核苷酸序列的 FASTA 文件，您可能想对非冗余 (NR) 蛋白质数据库运行 BLASTX（翻译）搜索。假设您（或您的系统管理员）已下载并安装了 NR 数据库，您可以运行
$ blastx -query opuntia.fasta -db nr -out opuntia.xml -evalue 0.001 -outfmt 5
这应该对 NR 数据库运行 BLASTX，使用 \(0.001\) 的期望截止值，并将 XML 输出生成到指定的输出文件（然后我们可以解析该文件）。在我的电脑上，这大约需要六分钟 - 这说明保存输出到文件是一个好主意，这样您就可以根据需要重复任何分析。
在 Python 中，我们可以使用 subprocess 模块来构建命令行字符串并运行它
>>> import subprocess
>>> cmd = "blastx -query opuntia.fasta -db nr -out opuntia.xml"
>>> cmd += " -evalue 0.001 -outfmt 5"
>>> subprocess.run(cmd, shell=True)
在这个例子中，BLASTX 应该不会向终端输出任何内容。您可能需要检查输出文件 opuntia.xml 是否已创建。
正如您可能从本教程中前面的示例中回忆到的，opuntia.fasta 包含七个序列，因此 BLAST XML 输出应包含多个结果。因此，使用 Bio.Blast.parse() 解析它，如下面的第 解析 BLAST 输出 节所述。
其他版本的 BLAST
NCBI BLAST+（用 C++ 编写）最初于 2009 年发布，以取代原始的 NCBI “传统” BLAST（用 C 编写），后者不再更新。您可能还会遇到 华盛顿大学 BLAST (WU-BLAST) 及其后续版本 高级生物计算 BLAST（AB-BLAST，于 2009 年发布，非免费/开源）。这些软件包包括命令行工具 wu-blastall 和 ab-blastall，它们模仿了 NCBI “传统” BLAST 套件中的 blastall。Biopython 目前没有提供调用这些工具的包装器，但应该能够解析来自它们的任何与 NCBI 兼容的输出。
解析 BLAST 输出
如上所述，BLAST 可以生成多种格式的输出，例如 XML、HTML 和纯文本。最初，Biopython 具有用于解析 BLAST 纯文本和 HTML 输出的解析器，因为这些是当时提供的唯一输出格式。这些解析器现在已从 Biopython 中删除，因为这些格式的 BLAST 输出一直在发生变化，每次都会破坏 Biopython 解析器。如今，Biopython 可以解析 XML 格式、XML2 格式和表格格式的 BLAST 输出。本章使用 Bio.Blast.parse 函数描述了用于解析 XML 和 XML2 格式的 BLAST 输出的解析器。此函数会自动检测 XML 文件是采用 XML 格式还是 XML2 格式。表格格式的 BLAST 输出可以作为对齐使用 Bio.Align.parse 函数解析（见 来自 BLAST 或 FASTA 的表格输出 节）。
您可以通过多种方式获取 XML 格式的 BLAST 输出。对于解析器来说，输出是如何生成的并不重要，只要它是 XML 格式即可。
您可以使用 Biopython 通过互联网运行 BLAST，如第 通过互联网运行 BLAST 节所述。
您可以使用 Biopython 在本地运行 BLAST，如第 在本地运行 BLAST 节所述。
您可以在 NCBI 网站上通过网络浏览器自己执行 BLAST 搜索，然后保存结果。您需要选择 XML 作为接收结果的格式，并将最终获得的 BLAST 页面（您知道的，就是包含所有有趣结果的那个页面！）保存到文件中。
您还可以不使用 Biopython 在本地运行 BLAST，并将输出保存到文件中。同样，您需要选择 XML 作为接收结果的格式。
重要的是，您不必使用 Biopython 脚本获取数据才能解析它。通过这些方法之一执行操作后，您需要为结果获取一个类文件对象。在 Python 中，类文件对象或句柄仅仅是对任何信息源输入的一种很好的通用描述方式，以便可以使用 read() 和 readline() 函数检索信息（见第 什么是句柄？ 节）。
如果您按照上面通过脚本与 BLAST 交互的代码操作，那么您已经拥有了 result_stream，即指向 BLAST 结果的类文件对象。例如，使用 GI 编号进行在线搜索
>>> from Bio import Blast
>>> result_stream = Blast.qblast("blastn", "nt", "8332116")
如果您通过其他方式运行 BLAST，并在文件 my_blast.xml 中拥有 BLAST 输出（采用 XML 格式），您只需要以可读方式（作为 bytes）打开该文件即可
>>> result_stream = open("my_blast.xml", "rb")
现在我们有了数据流，就可以准备解析输出。解析它的代码非常简单。如果您期望单个 BLAST 结果（即，您使用了一个查询）
>>> from Bio import Blast
>>> blast_record = Blast.read(result_stream)
或者，如果您有大量结果（即，多个查询序列）
>>> from Bio import Blast
>>> blast_records = Blast.parse(result_stream)
就像 Bio.SeqIO 和 Bio.Align（见第 序列输入/输出 章和第 序列对齐 章）一样，我们有一对输入函数，read 和 parse，其中 read 用于您只有一个对象时，而 parse 是一个迭代器，用于您可能拥有大量对象时 - 但我们得到的是 BLAST 记录对象，而不是 SeqRecord 或 Alignment 对象。
为了能够处理 BLAST 文件可能很大，包含数千个结果的情况，Blast.parse() 返回一个迭代器。通俗地说，迭代器允许您逐步遍历 BLAST 输出，为每个 BLAST 搜索结果逐个检索 BLAST 记录。
>>> from Bio import Blast
>>> blast_records = Blast.parse(result_stream)
>>> blast_record = next(blast_records)
# ... do something with blast_record
>>> blast_record = next(blast_records)
# ... do something with blast_record
>>> blast_record = next(blast_records)
# ... do something with blast_record
>>> blast_record = next(blast_records)
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
StopIteration
# No further records
或者，您可以使用 for 循环
>>> for blast_record in blast_records:
... pass # Do something with blast_record
...
但请注意，您只能遍历 BLAST 记录一次。通常，您会从每个 BLAST 记录中保存您感兴趣的信息。
或者，您可以将 blast_records 作为列表使用，例如通过索引提取一条记录，或者通过对 blast_records 调用 len 或 print。解析器随后会自动迭代记录并将它们存储起来。
>>> from Bio import Blast
>>> blast_records = Blast.parse("xml_2222_blastx_001.xml")
>>> len(blast_records) # this causes the parser to iterate over all records
7
>>> blast_records[2].query.description
'gi|5690369|gb|AF158246.1|AF158246 Cricetulus griseus glucose phosphate isomerase (GPI) gene, partial intron sequence'
但是，如果您的 BLAST 文件很大，您可能会遇到将所有记录保存在列表中时出现的内存问题。
如果您在将 blast_records 用作列表之前就开始迭代记录，解析器会首先将文件流重置到数据的开头，以确保所有记录都被读取。请注意，如果流无法重置到开头，例如，如果数据是从远程读取的（例如，通过 qblast；见 通过互联网运行 BLAST 小节），则这将失败。在这种情况下，您可以在迭代记录之前通过调用 blast_records = blast_records[:] 明确地将记录读入列表中。读取记录后，就可以安全地迭代它们或将它们用作列表。
您无需自己打开文件，只需提供文件名即可
>>> from Bio import Blast
>>> with Blast.parse("my_blast.xml") as blast_records:
... for blast_record in blast_records:
... pass # Do something with blast_record
...
在这种情况下，Biopython 会为您打开文件，并在不再需要文件时立即关闭它（虽然可以使用 blast_records = Blast.parse("my_blast.xml")，但它有一个缺点，即文件可能会比严格需要的更长时间保持打开状态，从而浪费资源）。
您可以 print 记录以快速了解其内容
>>> from Bio import Blast
>>> with Blast.parse("my_blast.xml") as blast_records:
... print(blast_records)
...
Program: BLASTN 2.2.27+
db: refseq_rna
Query: 42291 (length=61)
mystery_seq
Hits: ---- ----- ----------------------------------------------------------
# # HSP ID + description
---- ----- ----------------------------------------------------------
0 1 gi|262205317|ref|NR_030195.1| Homo sapiens microRNA 52...
1 1 gi|301171311|ref|NR_035856.1| Pan troglodytes microRNA...
2 1 gi|270133242|ref|NR_032573.1| Macaca mulatta microRNA ...
3 2 gi|301171322|ref|NR_035857.1| Pan troglodytes microRNA...
4 1 gi|301171267|ref|NR_035851.1| Pan troglodytes microRNA...
5 2 gi|262205330|ref|NR_030198.1| Homo sapiens microRNA 52...
6 1 gi|262205302|ref|NR_030191.1| Homo sapiens microRNA 51...
7 1 gi|301171259|ref|NR_035850.1| Pan troglodytes microRNA...
8 1 gi|262205451|ref|NR_030222.1| Homo sapiens microRNA 51...
9 2 gi|301171447|ref|NR_035871.1| Pan troglodytes microRNA...
10 1 gi|301171276|ref|NR_035852.1| Pan troglodytes microRNA...
11 1 gi|262205290|ref|NR_030188.1| Homo sapiens microRNA 51...
12 1 gi|301171354|ref|NR_035860.1| Pan troglodytes microRNA...
13 1 gi|262205281|ref|NR_030186.1| Homo sapiens microRNA 52...
14 2 gi|262205298|ref|NR_030190.1| Homo sapiens microRNA 52...
15 1 gi|301171394|ref|NR_035865.1| Pan troglodytes microRNA...
16 1 gi|262205429|ref|NR_030218.1| Homo sapiens microRNA 51...
17 1 gi|262205423|ref|NR_030217.1| Homo sapiens microRNA 52...
18 1 gi|301171401|ref|NR_035866.1| Pan troglodytes microRNA...
19 1 gi|270133247|ref|NR_032574.1| Macaca mulatta microRNA ...
20 1 gi|262205309|ref|NR_030193.1| Homo sapiens microRNA 52...
21 2 gi|270132717|ref|NR_032716.1| Macaca mulatta microRNA ...
22 2 gi|301171437|ref|NR_035870.1| Pan troglodytes microRNA...
23 2 gi|270133306|ref|NR_032587.1| Macaca mulatta microRNA ...
24 2 gi|301171428|ref|NR_035869.1| Pan troglodytes microRNA...
25 1 gi|301171211|ref|NR_035845.1| Pan troglodytes microRNA...
26 2 gi|301171153|ref|NR_035838.1| Pan troglodytes microRNA...
27 2 gi|301171146|ref|NR_035837.1| Pan troglodytes microRNA...
28 2 gi|270133254|ref|NR_032575.1| Macaca mulatta microRNA ...
29 2 gi|262205445|ref|NR_030221.1| Homo sapiens microRNA 51...
~~~
97 1 gi|356517317|ref|XM_003527287.1| PREDICTED: Glycine ma...
98 1 gi|297814701|ref|XM_002875188.1| Arabidopsis lyrata su...
99 1 gi|397513516|ref|XM_003827011.1| PREDICTED: Pan panisc...
通常，您一次只运行一个 BLAST 搜索。然后，您只需从 blast_records 中取出第一个（也是唯一一个）BLAST 记录
>>> from Bio import Blast
>>> blast_records = Blast.parse("my_blast.xml")
>>> blast_record = next(blast_records)
或者更简洁地
>>> from Bio import Blast
>>> blast_record = Blast.read(result_stream)
或者，等效地
>>> from Bio import Blast
>>> blast_record = Blast.read("my_blast.xml")
（在这里，您不需要使用 with 块，因为 Blast.read 会读取整个文件并在之后立即关闭它）。
我想您现在想知道 BLAST 记录中包含什么内容。
BLAST Records、Record 和 Hit 类
BLAST Records 类
单个 BLAST 输出文件可以包含来自多个 BLAST 查询的输出。在 Biopython 中，BLAST 输出文件中的信息存储在一个 Bio.Blast.Records 对象中。这是一个迭代器，它为每个查询返回一个 Bio.Blast.Record 对象（见 BLAST Record 类 小节）。Bio.Blast.Records 对象具有以下描述 BLAST 运行的属性
source: 用于构建Bio.Blast.Records对象的输入数据（可以是文件名或路径，也可以是类似文件的对象）。program: 使用的特定 BLAST 程序（例如，“blastn”）。version: BLAST 程序的版本（例如，“BLASTN 2.2.27+”）。reference: BLAST 发表的文献引用。db: 用于运行查询的 BLAST 数据库（例如，“nr”）。query: 一个SeqRecord对象，可能包含以下部分或全部信息query.id: 查询的 SeqId；query.description: 查询的定义行；query.seq: 查询序列。
param: 包含 BLAST 运行中使用的参数的字典。您可能会在该字典中找到以下键'matrix': BLAST 运行中使用的评分矩阵（例如，“BLOSUM62”）（字符串）；'expect': 随机匹配预期数量的阈值（浮点数）；'include': psiblast 中多通道模型中包含的 e 值阈值（浮点数）；'sc-match': 匹配核苷酸的得分（整数）；'sc-mismatch': 不匹配核苷酸的得分（整数；'gap-open': 间隙打开成本（整数）；'gap-extend': 间隙扩展成本（整数）；'filter': BLAST 运行中应用的过滤选项（字符串）；'pattern': PHI-BLAST 模式（字符串）；'entrez-query': 对 Entrez 查询的请求限制（字符串）。
mbstat: 包含 Mega BLAST 搜索统计信息的字典。有关该字典中项目的描述，请参阅下面Record.stat属性的描述（在 BLAST 记录类 小节中）。只有旧版本的 Mega BLAST 会存储此信息。由于它存储在 BLAST 输出文件的末尾，因此只有在文件完全读取后（通过遍历记录直到发出StopIteration）才能访问该属性。
对于我们的示例，我们发现
>>> blast_records
<Bio.Blast.Records source='my_blast.xml' program='blastn' version='BLASTN 2.2.27+' db='refseq_rna'>
>>> blast_records.source
'my_blast.xml'
>>> blast_records.program
'blastn'
>>> blast_records.version
'BLASTN 2.2.27+'
>>> blast_records.reference
'Stephen F. Altschul, Thomas L. Madden, Alejandro A. Schäffer, Jinghui Zhang, Zheng Zhang, Webb Miller, and David J. Lipman (1997), "Gapped BLAST and PSI-BLAST: a new generation of protein database search programs", Nucleic Acids Res. 25:3389-3402.'
>>> blast_records.db
'refseq_rna'
>>> blast_records.param
{'expect': 10.0, 'sc-match': 2, 'sc-mismatch': -3, 'gap-open': 5, 'gap-extend': 2, 'filter': 'L;m;'}
>>> print(blast_records)
Program: BLASTN 2.2.27+
db: refseq_rna
Query: 42291 (length=61)
mystery_seq
Hits: ---- ----- ----------------------------------------------------------
# # HSP ID + description
---- ----- ----------------------------------------------------------
0 1 gi|262205317|ref|NR_030195.1| Homo sapiens microRNA 52...
1 1 gi|301171311|ref|NR_035856.1| Pan troglodytes microRNA...
2 1 gi|270133242|ref|NR_032573.1| Macaca mulatta microRNA ...
3 2 gi|301171322|ref|NR_035857.1| Pan troglodytes microRNA...
...
BLAST 记录类
一个 Bio.Blast.Record 对象存储 BLAST 为单个查询提供的的信息。 Bio.Blast.Record 类继承自 list，本质上是一个 Bio.Blast.Hit 对象列表（请参阅 BLAST 命中类 部分）。 Bio.Blast.Record 对象具有以下两个属性
query: 一个SeqRecord对象，可能包含以下部分或全部信息query.id: 查询的 SeqId；query.description: 查询的定义行；query.seq: 查询序列。
stat: 包含 BLAST 命中统计数据的字典。您可能会在该字典中找到以下键'db-num': BLAST 数据库中的序列数量（整数）；'db-len': BLAST 数据库的长度（整数）；'hsp-len': 有效 HSP（高得分对）长度（整数）；'eff-space': 有效搜索空间（浮点数）；'kappa': Karlin-Altschul 参数 K（浮点数）；'lambda': Karlin-Altschul 参数 Lambda（浮点数）；'entropy': Karlin-Altschul 参数 H（浮点数）
message: 一些（错误？）信息。
继续我们的示例，
>>> blast_record
<Bio.Blast.Record query.id='42291'; 100 hits>
>>> blast_record.query
SeqRecord(seq=Seq(None, length=61), id='42291', name='<unknown name>', description='mystery_seq', dbxrefs=[])
>>> blast_record.stat
{'db-num': 3056429, 'db-len': 673143725, 'hsp-len': 0, 'eff-space': 0, 'kappa': 0.41, 'lambda': 0.625, 'entropy': 0.78}
>>> print(blast_record)
Query: 42291 (length=61)
mystery_seq
Hits: ---- ----- ----------------------------------------------------------
# # HSP ID + description
---- ----- ----------------------------------------------------------
0 1 gi|262205317|ref|NR_030195.1| Homo sapiens microRNA 52...
1 1 gi|301171311|ref|NR_035856.1| Pan troglodytes microRNA...
2 1 gi|270133242|ref|NR_032573.1| Macaca mulatta microRNA ...
3 2 gi|301171322|ref|NR_035857.1| Pan troglodytes microRNA...
...
由于 Bio.Blast.Record 类继承自 list，因此您可以像使用它一样使用它。例如，您可以遍历记录
>>> for hit in blast_record:
... hit
...
<Bio.Blast.Hit target.id='gi|262205317|ref|NR_030195.1|' query.id='42291'; 1 HSP>
<Bio.Blast.Hit target.id='gi|301171311|ref|NR_035856.1|' query.id='42291'; 1 HSP>
<Bio.Blast.Hit target.id='gi|270133242|ref|NR_032573.1|' query.id='42291'; 1 HSP>
<Bio.Blast.Hit target.id='gi|301171322|ref|NR_035857.1|' query.id='42291'; 2 HSPs>
<Bio.Blast.Hit target.id='gi|301171267|ref|NR_035851.1|' query.id='42291'; 1 HSP>
...
要检查 blast_record 有多少个命中，您只需调用 Python 的 len 函数
>>> len(blast_record)
100
与 Python 列表类似，您可以使用索引从 Bio.Blast.Record 中检索命中
>>> blast_record[0] # retrieves the top hit
<Bio.Blast.Hit target.id='gi|262205317|ref|NR_030195.1|' query.id='42291'; 1 HSP>
>>> blast_record[-1] # retrieves the last hit
<Bio.Blast.Hit target.id='gi|397513516|ref|XM_003827011.1|' query.id='42291'; 1 HSP>
要从 Bio.Blast.Record 中检索多个命中，您可以使用切片表示法。这将返回一个新的 Bio.Blast.Record 对象，该对象仅包含切片后的命中
>>> blast_slice = blast_record[:3] # slices the first three hits
>>> print(blast_slice)
Query: 42291 (length=61)
mystery_seq
Hits: ---- ----- ----------------------------------------------------------
# # HSP ID + description
---- ----- ----------------------------------------------------------
0 1 gi|262205317|ref|NR_030195.1| Homo sapiens microRNA 52...
1 1 gi|301171311|ref|NR_035856.1| Pan troglodytes microRNA...
2 1 gi|270133242|ref|NR_032573.1| Macaca mulatta microRNA ...
要创建 Bio.Blast.Record 的副本，请获取完整的切片
>>> blast_record_copy = blast_record[:]
>>> type(blast_record_copy)
<class 'Bio.Blast.Record'>
>>> blast_record_copy # list of all hits
<Bio.Blast.Record query.id='42291'; 100 hits>
这在您要对 BLAST 记录进行排序或过滤（请参阅 排序和过滤 BLAST 输出）但又希望保留原始 BLAST 输出的副本时特别有用。
您还可以将 blast_record 作为 Python 字典访问，并使用命中的 ID 作为键检索命中
>>> blast_record["gi|262205317|ref|NR_030195.1|"]
<Bio.Blast.Hit target.id='gi|262205317|ref|NR_030195.1|' query.id='42291'; 1 HSP>
如果在 blast_record 中找不到 ID，则会引发 KeyError
>>> blast_record["unicorn_gene"]
Traceback (most recent call last):
...
KeyError: 'unicorn_gene'
您可以通过通常使用 .keys() 来获取完整键列表
>>> blast_record.keys()
['gi|262205317|ref|NR_030195.1|', 'gi|301171311|ref|NR_035856.1|', 'gi|270133242|ref|NR_032573.1|', ...]
如果您只想检查特定命中是否出现在查询结果中怎么办？您可以使用 in 关键字进行简单的 Python 成员资格测试
>>> "gi|262205317|ref|NR_030195.1|" in blast_record
True
>>> "gi|262205317|ref|NR_030194.1|" in blast_record
False
有时，仅仅知道命中是否存在还不够；您还需要知道命中的排名。在这里， index 方法可以帮到您
>>> blast_record.index("gi|301171437|ref|NR_035870.1|")
22
请记住，Python 使用基于零的索引，因此第一个命中将位于索引 0。
BLAST 命中类
blast_record 列表中的每个 Bio.Blast.Hit 对象表示查询与目标的一次 BLAST 命中。
>>> hit = blast_record[0]
>>> hit
<Bio.Blast.Hit target.id='gi|262205317|ref|NR_030195.1|' query.id='42291'; 1 HSP>
>>> hit.target
SeqRecord(seq=Seq(None, length=61), id='gi|262205317|ref|NR_030195.1|', name='NR_030195', description='Homo sapiens microRNA 520b (MIR520B), microRNA', dbxrefs=[])
我们可以通过打印它来获得 hit 的摘要
>>> print(blast_record[3])
Query: 42291
mystery_seq
Hit: gi|301171322|ref|NR_035857.1| (length=86)
Pan troglodytes microRNA mir-520c (MIR520C), microRNA
HSPs: ---- -------- --------- ------ --------------- ---------------------
# E-value Bit score Span Query range Hit range
---- -------- --------- ------ --------------- ---------------------
0 8.9e-20 100.47 60 [1:61] [13:73]
1 3.3e-06 55.39 60 [0:60] [73:13]
您会看到我们已经涵盖了基本内容
命中始终针对一个查询；查询 ID 和描述显示在摘要的顶部。
命中包含查询与一个目标序列的一个或多个比对。目标信息显示在摘要中的下一个位置。如上所示，目标可以通过命中的
target属性访问。最后，还有一个表格包含有关每个命中包含的比对的快速信息。在 BLAST 行话中，这些比对被称为“高得分片段对”，或 HSP（请参阅 BLAST HSP 类 部分）。表中的每一行都总结了一个 HSP，包括 HSP 索引、e 值、位得分、跨度（包括间隙的比对长度）、查询坐标和目标坐标。
Bio.Blast.Hit 类是 Bio.Align.Alignments（复数；请参阅 Alignments 类 部分）的子类，因此本质上是一个 Bio.Align.Alignment（单数；请参阅 Alignment 对象 部分）对象列表。特别是当将核苷酸序列比对到基因组时，如果特定查询比对到染色体的多个区域，则 Bio.Blast.Hit 对象可能包含多个 Bio.Align.Alignment。对于蛋白质比对，通常一个命中只包含一个比对，尤其是对于高度同源序列的比对。
>>> type(hit)
<class 'Bio.Blast.Hit'>
>>> from Bio.Align import Alignments
>>> isinstance(hit, Alignments)
True
>>> len(hit)
1
对于 XML2 格式的 BLAST 输出，命中可能具有多个具有相同序列但不同序列 ID 和描述的目标。这些目标可以作为 hit.targets 属性访问。在大多数情况下， hit.targets 的长度为 1，并且仅包含 hit.target
>>> from Bio import Blast
>>> blast_record = Blast.read("xml_2900_blastx_001_v2.xml")
>>> for hit in blast_record:
... print(len(hit.targets))
...
1
1
2
1
1
1
1
1
1
1
但是，正如您在上面的输出中看到的，第三个命中具有多个目标。
>>> hit = blast_record[2]
>>> hit.targets[0].seq
Seq(None, length=246)
>>> hit.targets[1].seq
Seq(None, length=246)
>>> hit.targets[0].id
'gi|684409690|ref|XP_009175831.1|'
>>> hit.targets[1].id
'gi|663044098|gb|KER20427.1|'
>>> hit.targets[0].name
'XP_009175831'
>>> hit.targets[1].name
'KER20427'
>>> hit.targets[0].description
'hypothetical protein T265_11027 [Opisthorchis viverrini]'
>>> hit.targets[1].description
'hypothetical protein T265_11027 [Opisthorchis viverrini]'
由于两个目标的序列内容彼此相同，因此它们的序列比对也相同。因此，此命中的比对仅引用 hit.targets[0]（与 hit.target 相同），因为 hit.targets[1] 的比对无论如何都会相同。
BLAST HSP 类
让我们回到我们的主要示例，看看第一个命中中的第一个（也是唯一一个）比对。此比对是 Bio.Blast.HSP 类的实例，它是 Bio.Align 中的 Alignment 类的子类
>>> from Bio import Blast
>>> blast_record = Blast.read("my_blast.xml")
>>> hit = blast_record[0]
>>> len(hit)
1
>>> alignment = hit[0]
>>> alignment
<Bio.Blast.HSP target.id='gi|262205317|ref|NR_030195.1|' query.id='42291'; 2 rows x 61 columns>
>>> type(alignment)
<class 'Bio.Blast.HSP'>
>>> from Bio.Align import Alignment
>>> isinstance(alignment, Alignment)
True
alignment 对象具有指向目标和查询序列的属性，以及一个 coordinates 属性，用于描述序列比对。
>>> alignment.target
SeqRecord(seq=Seq('CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCCTTT...GGG'), id='gi|262205317|ref|NR_030195.1|', name='NR_030195', description='Homo sapiens microRNA 520b (MIR520B), microRNA', dbxrefs=[])
>>> alignment.query
SeqRecord(seq=Seq('CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCCTTT...GGG'), id='42291', name='<unknown name>', description='mystery_seq', dbxrefs=[])
>>> alignment.target
SeqRecord(seq=Seq('CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCCTTT...GGG'), id='gi|262205317|ref|NR_030195.1|', name='NR_030195', description='Homo sapiens microRNA 520b (MIR520B), microRNA', dbxrefs=[])
>>> alignment.query
SeqRecord(seq=Seq('CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCCTTT...GGG'), id='42291', name='<unknown name>', description='mystery_seq', dbxrefs=[])
>>> print(alignment.coordinates)
[[ 0 61]
[ 0 61]]
对于翻译的 BLAST 搜索，目标或查询的 features 属性可能包含一个类型为 CDS 的 SeqFeature，该属性存储氨基酸序列区域。此特征的 qualifiers 属性是一个字典，其中只有一个键 'coded_by'；相应的 value 指定编码氨基酸序列的核苷酸序列区域，以 GenBank 风格的字符串形式表示，并带有基于 1 的坐标。
每个 Alignment 对象都有以下附加属性
score: 高得分对 (HSP) 的得分；annotations: 可能包含以下键的字典'bit score': HSP 的得分（以比特表示）（浮点数）；'evalue': HSP 的 E 值（浮点数）；'identity'’: HSP 中的相同性数量（整数）；'positive': HSP 中的正值数量（整数）；'gaps': HSP 中的间隙数量（整数）；'midline': 格式化中间线。
通常的 Alignment 方法（参见第 对齐对象 节）可以应用于 alignment。例如，我们可以打印对齐
>>> print(alignment)
Query : 42291 Length: 61 Strand: Plus
mystery_seq
Target: gi|262205317|ref|NR_030195.1| Length: 61 Strand: Plus
Homo sapiens microRNA 520b (MIR520B), microRNA
Score:111 bits(122), Expect:5e-23,
Identities:61/61(100%), Positives:61/61(100%), Gaps:0.61(0%)
gi|262205 0 CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCCTTTTAGAGG
0 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
42291 0 CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCCTTTTAGAGG
gi|262205 60 G 61
60 | 61
42291 60 G 61
我们只打印 BLAST 记录中所有大于特定阈值的命中的一些摘要信息
>>> E_VALUE_THRESH = 0.04
>>> for alignments in blast_record:
... for alignment in alignments:
... if alignment.annotations["evalue"] < E_VALUE_THRESH:
... print("****Alignment****")
... print("sequence:", alignment.target.id, alignment.target.description)
... print("length:", len(alignment.target))
... print("score:", alignment.score)
... print("e value:", alignment.annotations["evalue"])
... print(alignment[:, :50])
...
****Alignment****
sequence: gi|262205317|ref|NR_030195.1| Homo sapiens microRNA 520b (MIR520B), microRNA
length: 61
score: 122.0
e value: 4.91307e-23
gi|262205 0 CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTC 50
0 |||||||||||||||||||||||||||||||||||||||||||||||||| 50
42291 0 CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTC 50
****Alignment****
sequence: gi|301171311|ref|NR_035856.1| Pan troglodytes microRNA mir-520b (MIR520B), microRNA
length: 60
score: 120.0
e value: 1.71483e-22
gi|301171 0 CCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCC 50
0 |||||||||||||||||||||||||||||||||||||||||||||||||| 50
42291 1 CCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCC 51
****Alignment****
sequence: gi|270133242|ref|NR_032573.1| Macaca mulatta microRNA mir-519a (MIR519A), microRNA
length: 85
score: 112.0
e value: 2.54503e-20
gi|270133 12 CCCTCTAGAGGGAAGCGCTTTCTGTGGTCTGAAAGAAAAGAAAGTGCTTC 62
0 |||||||.|||||||||||||||||.|||||||||||||||||||||||| 50
42291 0 CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTC 50
...
对 BLAST 输出进行排序和过滤
如果 BLAST 输出文件中的命中顺序不符合您的口味，您可以使用 sort 方法重新排列 Bio.Blast.Record 对象中的命中。例如，这里我们根据每个目标的序列长度对命中进行排序，将 reverse 标志设置为 True，以便我们按降序排序。
>>> for hit in blast_record[:5]:
... print(f"{hit.target.id} {len(hit.target)}")
...
gi|262205317|ref|NR_030195.1| 61
gi|301171311|ref|NR_035856.1| 60
gi|270133242|ref|NR_032573.1| 85
gi|301171322|ref|NR_035857.1| 86
gi|301171267|ref|NR_035851.1| 80
>>> sort_key = lambda hit: len(hit.target)
>>> blast_record.sort(key=sort_key, reverse=True)
>>> for hit in blast_recordt[:5]:
... print(f"{hit.target.id} {len(hit.target)}")
...
gi|397513516|ref|XM_003827011.1| 6002
gi|390332045|ref|XM_776818.2| 4082
gi|390332043|ref|XM_003723358.1| 4079
gi|356517317|ref|XM_003527287.1| 3251
gi|356543101|ref|XM_003539954.1| 2936
这将在内部对 blast_record 进行排序。如果要保留原始的、未排序的 BLAST 输出的副本，请在排序之前使用 original_blast_record = blast_record[:]。
要根据 BLAST 命中的属性过滤 BLAST 命中，您可以使用 Python 内置的 filter 以及相应的回调函数来评估每个命中。回调函数必须以单个 Hit 对象作为其参数，并返回 True 或 False。以下是一个示例，我们从中过滤掉只有单个 HSP 的 Hit 对象
>>> filter_func = lambda hit: len(hit) > 1 # the callback function
>>> len(blast_record) # no. of hits before filtering
100
>>> blast_record[:] = filter(filter_func, blast_record)
>>> len(blast_record) # no. of hits after filtering
37
>>> for hit in blast_record[:5]: # quick check for the hit lengths
... print(f"{hit.target.id} {len(hit)}")
...
gi|301171322|ref|NR_035857.1| 2
gi|262205330|ref|NR_030198.1| 2
gi|301171447|ref|NR_035871.1| 2
gi|262205298|ref|NR_030190.1| 2
gi|270132717|ref|NR_032716.1| 2
类似地，您可以过滤每个命中中的 HSP，例如根据它们的 E 值
>>> filter_func = lambda hsp: hsp.annotations["evalue"] < 1.0e-12
>>> for hit in blast_record:
... hit[:] = filter(filter_func, hit)
...
可能您希望在之后删除所有没有剩余 HSP 的命中
>>> filter_func = lambda hit: len(hit) > 0
>>> blast_record[:] = filter(filter_func, blast_record)
>>> len(blast_record)
16
使用 Python 内置的 map 函数修改 BLAST 记录中的命中或 HSP。 map 函数接受一个回调函数，该函数返回修改后的命中对象。例如，我们可以使用 map 重命名命中 ID
>>> for hit in blast_record[:5]:
... print(hit.target.id)
...
gi|301171322|ref|NR_035857.1|
gi|262205330|ref|NR_030198.1|
gi|301171447|ref|NR_035871.1|
gi|262205298|ref|NR_030190.1|
gi|270132717|ref|NR_032716.1|
>>> import copy
>>> original_blast_record = copy.deepcopy(blast_record)
>>> def map_func(hit):
... # renames "gi|301171322|ref|NR_035857.1|" to "NR_035857.1"
... hit.target.id = hit.target.id.split("|")[3]
... return hit
...
>>> blast_record[:] = map(map_func, blast_record)
>>> for hit in blast_record[:5]:
... print(hit.target.id)
...
NR_035857.1
NR_030198.1
NR_035871.1
NR_030190.1
NR_032716.1
>>> for hit in original_blast_record[:5]:
... print(hit.target.id)
...
gi|301171322|ref|NR_035857.1|
gi|262205330|ref|NR_030198.1|
gi|301171447|ref|NR_035871.1|
gi|262205298|ref|NR_030190.1|
gi|270132717|ref|NR_032716.1|
请注意，在本例中，map_func 在内部修改了命中。与排序和过滤（参见上文）相比，使用 original_blast_record = blast_record[:] 对于保留未修改的 BLAST 记录的副本并不足够，因为它创建了 BLAST 记录的浅层副本，该副本包含指向相同 Hit 对象的指针。相反，我们使用 copy.deepcopy 创建 BLAST 记录的副本，其中每个 Hit 对象都被复制了。
写入 BLAST 记录
使用 Bio.Blast 中的 write 函数将 BLAST 记录保存为 XML 文件。默认情况下，使用（基于 DTD 的）XML 格式；您也可以通过使用 fmt="XML2" 参数为 write 函数，以（基于模式的）XML2 格式保存 BLAST 记录。
>>> from Bio import Blast
>>> stream = Blast.qblast("blastn", "nt", "8332116")
>>> records = Blast.parse(stream)
>>> Blast.write(records, "my_qblast_output.xml")
or
>>> Blast.write(records, "my_qblast_output.xml", fmt="XML2")
在本例中，我们可以直接将 Blast.qblast 返回的数据保存到 XML 文件中（参见第 保存 BLAST 结果 节）。但是，通过将 qblast 返回的数据解析为记录，我们可以在保存 BLAST 记录之前对它们进行排序或过滤。例如，我们可能只对正得分至少为 400 的 BLAST HSP 感兴趣
>>> filter_func = lambda hsp: hsp.annotations["positive"] >= 400
>>> for hit in records[0]:
... hit[:] = filter(filter_func, hit)
...
>>> Blast.write(records, "my_qblast_output_selected.xml")
除了文件名之外，Blast.write 的第二个参数还可以是文件流。在这种情况下，该流必须以二进制格式打开以供写入
>>> with open("my_qblast_output.xml", "wb") as stream:
... Blast.write(records, stream)
...
处理 PSI-BLAST
您可以直接从命令行或使用 python 的 subprocess 模块运行 PSI-BLAST 的独立版本（psiblast）。
在撰写本文时，NCBI 似乎不支持通过互联网运行 PSI-BLAST 搜索的工具。
请注意，Bio.Blast 解析器可以读取当前版本的 PSI-BLAST 的 XML 输出，但是像每次迭代中哪些序列是新的或重复使用的信息在 XML 文件中并不存在。
处理 RPS-BLAST
您可以直接从命令行或使用 python 的 subprocess 模块运行 RPS-BLAST 的独立版本（rpsblast）。
在撰写本文时，NCBI 似乎不支持通过互联网运行 RPS-BLAST 搜索的工具。
您可以使用 Bio.Blast 解析器读取当前版本的 RPS-BLAST 的 XML 输出。
BLAST（旧版）
嗨，每个人都喜欢 BLAST，对吧？我的意思是，天哪，还能有什么比将你的序列与已知世界中的所有其他序列进行比较更容易的事了？但是，当然，这一部分不是关于 BLAST 有多酷，因为我们已经知道这一点。而是关于 BLAST 的问题——处理大型运行产生的海量数据，以及自动化 BLAST 运行，可能非常困难。
幸运的是，Biopython 的开发人员非常了解这一点，因此他们开发了许多用于处理 BLAST 并使其变得更轻松的工具。本节详细介绍了如何使用这些工具以及如何用它们做一些有用的事情。
处理 BLAST 可以分为两个步骤，这两个步骤都可以在 Biopython 中完成。首先，为你的查询序列运行 BLAST，并获取一些输出。其次，在 Python 中解析 BLAST 输出以进行进一步分析。
您第一次接触到运行 BLAST 可能是通过 NCBI 网页服务。事实上，您可以通过多种方式运行 BLAST，这些方式可以分为几种类别。最重要的区别是在本地（在您自己的机器上）运行 BLAST，以及在远程（在另一台机器上，通常是 NCBI 服务器）运行 BLAST。我们将在本章开始时从 Python 脚本中调用 NCBI 在线 BLAST 服务。
注意：以下第 BLAST 和其他序列搜索工具 章介绍了 Bio.SearchIO。我们希望它最终能够取代旧的 Bio.Blast 模块，因为它提供了一个更通用的框架，可以处理其他相关的序列搜索工具。但是，目前您可以使用该模块或旧的 Bio.Blast 模块来处理 NCBI BLAST。
通过互联网运行 BLAST
我们使用 Bio.Blast.NCBIWWW 模块中的 qblast() 函数调用 BLAST 的在线版本。它有三个非可选参数
第一个参数是用于搜索的 BLAST 程序，以小写字符串形式表示。您可以在 https://blast.ncbi.nlm.nih.gov/Blast.cgi 中查看程序的选项和说明。目前
qblast仅适用于 blastn、blastp、blastx、tblast 和 tblastx。第二个参数指定要搜索的数据库。同样，您可以在 NCBI BLAST 指南中查看此选项的列表：https://blast.ncbi.nlm.nih.gov/doc/blast-help/.
第三个参数是一个字符串，包含你的查询序列。它可以是序列本身、fasta 格式的序列，或类似 GI 号码的标识符。
NCBI 指南中的 https://blast.ncbi.nlm.nih.gov/doc/blast-help/developerinfo.html#developerinfo 指出
不要比每 10 秒一次更频繁地联系服务器。
不要比每分钟一次更频繁地轮询任何单个 RID。
使用 URL 参数 email 和 tool，以便 NCBI 在出现问题时可以联系你。
如果要提交 50 次以上搜索，请在周末或工作日晚上 9 点到次日凌晨 5 点之间运行脚本。
为了满足第三点，可以设置 NCBIWWW.email 变量。
>>> from Bio.Blast import NCBIWWW
>>> NCBIWWW.email = "[email protected]"
qblast 函数还接受许多其他可选参数，这些参数基本上类似于你可以在 BLAST 网页上设置的不同参数。我们这里只重点介绍其中一些。
url_base参数设置用于通过互联网运行 BLAST 的基本 URL。默认情况下，它连接到 NCBI，但你可以使用它来连接到云中运行的 NCBI BLAST 实例。有关更多详细信息，请参阅qblast函数的文档。qblast函数可以以多种格式返回 BLAST 结果，您可以使用可选的format_type关键字选择："HTML"、"Text"、"ASN.1"或"XML"。默认值为"XML"，因为这是解析器期望的格式，在第 解析 BLAST 输出 节中进行了描述。expect参数设置期望值或 e 值阈值。
有关可选 BLAST 参数的更多信息，请参阅 NCBI 自身的文档，或者 Biopython 中内置的文档。
>>> from Bio.Blast import NCBIWWW
>>> help(NCBIWWW.qblast)
请注意，NCBI BLAST 网站上的默认设置与 QBLAST 上的默认设置略有不同。如果得到不同的结果，则需要检查参数（例如，期望值阈值和间隙值）。
例如，如果你有一个核苷酸序列，你想要使用 BLASTN 对其进行搜索，并且你了解你的查询序列的 GI 号码，你可以使用
>>> from Bio.Blast import NCBIWWW
>>> result_handle = NCBIWWW.qblast("blastn", "nt", "8332116")
或者，如果我们已经将查询序列保存在 FASTA 格式的文件中，我们只需要打开文件并将此记录读入为字符串，并将其用作查询参数即可
>>> from Bio.Blast import NCBIWWW
>>> fasta_string = open("m_cold.fasta").read()
>>> result_handle = NCBIWWW.qblast("blastn", "nt", fasta_string)
我们也可以将 FASTA 文件读入为 SeqRecord，然后仅提供序列本身
>>> from Bio.Blast import NCBIWWW
>>> from Bio import SeqIO
>>> record = SeqIO.read("m_cold.fasta", format="fasta")
>>> result_handle = NCBIWWW.qblast("blastn", "nt", record.seq)
只提供序列意味着 BLAST 会自动为您的序列分配一个标识符。您可能更喜欢使用 SeqRecord 对象的格式方法来生成 FASTA 字符串（其中将包括现有的标识符）
>>> from Bio.Blast import NCBIWWW
>>> from Bio import SeqIO
>>> record = SeqIO.read("m_cold.fasta", format="fasta")
>>> result_handle = NCBIWWW.qblast("blastn", "nt", record.format("fasta"))
如果你将序列保存在非 FASTA 文件格式中，你可以使用 Bio.SeqIO（参见第 序列输入/输出 章）提取这些序列，这种方法更有意义。
无论您为 qblast() 函数提供什么参数，您都应该在句柄对象中获得结果（默认情况下为 XML 格式）。下一步是将 XML 输出解析为表示搜索结果的 Python 对象（第 解析 BLAST 输出 节），但您可能希望首先保存输出文件的本地副本。我发现这在调试从 BLAST 结果中提取信息的代码时特别有用（因为重新运行在线搜索很慢，并且会浪费 NCBI 的计算时间）。
我们需要小心一点，因为我们只能使用 result_handle.read() 读取 BLAST 输出一次——再次调用 result_handle.read() 将返回空字符串。
>>> with open("my_blast.xml", "w") as out_handle:
... out_handle.write(result_handle.read())
...
>>> result_handle.close()
完成此操作后，结果将存储在文件 my_blast.xml 中，并且原始句柄已提取了所有数据（因此我们关闭了它）。但是，BLAST 解析器的 parse 函数（在 解析 BLAST 输出 中进行了描述）接受类似文件句柄的对象，因此我们可以直接打开保存的文件以供输入
>>> result_handle = open("my_blast.xml")
现在我们已将 BLAST 结果重新保存到句柄中，我们已准备好对其进行操作，因此这将我们直接引导到解析部分（参见下面的第 解析 BLAST 输出 节）。您现在可能希望直接跳到该部分……。
在本地运行 BLAST
简介
在本地运行 BLAST（而不是通过互联网，参见第 通过互联网运行 BLAST 节）至少有两个主要优势
本地 BLAST 可能比通过互联网运行 BLAST 更快；
本地 BLAST 允许你创建自己的数据库来搜索序列。
在本地运行 BLAST 的另一个原因可能是处理专有或未公开的序列数据。您可能不被允许重新发布这些序列，因此将它们提交给 NCBI 作为 BLAST 查询将不可行。
不幸的是，也有一些主要的缺点 - 安装所有组件并使其正确设置需要一些努力。
本地 BLAST 需要安装命令行工具。
本地 BLAST 需要设置（并可能保持更新）大型 BLAST 数据库。
为了进一步说明问题，有几个不同的 BLAST 软件包可用，并且还有一些其他工具可以生成类似 BLAST 输出文件的输出，例如 BLAT。
独立的 NCBI BLAST+
“新的”NCBI BLAST+ 套件于 2009 年发布。它取代了旧的 NCBI “遗留” BLAST 软件包（参见下文）。
本节将简要介绍如何从 Python 中使用这些工具。如果您已经阅读或尝试过第 对齐工具 节中的对齐工具示例，那么这一切都应该看起来相当简单。首先，我们构建一个命令行字符串（就像您在运行独立 BLAST 时手动在命令行提示符中输入的那样）。然后，我们可以在 Python 中执行此命令。
例如，如果您有一个包含基因核苷酸序列的 FASTA 文件，您可能想对非冗余 (NR) 蛋白质数据库运行 BLASTX（翻译）搜索。假设您（或您的系统管理员）已下载并安装了 NR 数据库，您可以运行
$ blastx -query opuntia.fasta -db nr -out opuntia.xml -evalue 0.001 -outfmt 5
这应该对 NR 数据库运行 BLASTX，使用 \(0.001\) 的期望截止值，并将 XML 输出生成到指定的输出文件（然后我们可以解析该文件）。在我的电脑上，这大约需要六分钟 - 这说明保存输出到文件是一个好主意，这样您就可以根据需要重复任何分析。
在 Python 中，我们可以使用 subprocess 模块来构建命令行字符串并运行它
>>> import subprocess
>>> cmd = "blastx -query opuntia.fasta -db nr -out opuntia.xml"
>>> cmd += " -evalue 0.001 -outfmt 5"
>>> subprocess.run(cmd, shell=True)
在这个例子中，BLASTX 应该不会向终端输出任何内容。您可能需要检查输出文件 opuntia.xml 是否已创建。
正如你在本教程之前示例中所看到的，opuntia.fasta 文件包含七个序列，因此 BLAST XML 输出应该包含多个结果。因此，使用 Bio.Blast.NCBIXML.parse() 解析它，如下面第 解析 BLAST 输出 节所述。
其他版本的 BLAST
NCBI BLAST+（用 C++ 编写）于 2009 年首次发布，用于替换原始的 NCBI “传统” BLAST（用 C 编写），该版本不再更新。发生了很多变化——旧版本有一个单独的核心命令行工具 blastall，涵盖了多种不同的 BLAST 搜索类型（现在在 BLAST+ 中是独立的命令），所有命令行选项都被重新命名。Biopython 对 NCBI “传统” BLAST 工具的封装已被弃用，将在未来版本中删除。为了避免混淆，本教程不涵盖从 Biopython 调用这些旧工具的内容。
你可能还会遇到 华盛顿大学 BLAST (WU-BLAST) 及其继任者，高级生物计算 BLAST (AB-BLAST，于 2009 年发布，不是免费/开源软件)。这些软件包包括命令行工具 wu-blastall 和 ab-blastall，它们模仿了 NCBI “传统” BLAST 套件中的 blastall。Biopython 目前没有提供用于调用这些工具的封装，但应该能够解析来自它们的任何与 NCBI 兼容的输出。
解析 BLAST 输出
如上所述，BLAST 可以生成多种格式的输出，例如 XML、HTML 和纯文本。最初，Biopython 具有用于解析 BLAST 纯文本和 HTML 输出的解析器，因为这些是当时提供的唯一输出格式。不幸的是，这些格式中的 BLAST 输出不断发生变化，每次都会破坏 Biopython 解析器。我们的 HTML BLAST 解析器已被删除，而已弃用的纯文本 BLAST 解析器现在只能通过 Bio.SearchIO 使用。请自行承担风险使用它，它可能有效，也可能无效，具体取决于你使用的 BLAST 版本。
由于跟上 BLAST 的变化越来越困难，特别是对于运行不同 BLAST 版本的用户，我们现在建议解析 XML 格式的输出，该格式可以由最近的 BLAST 版本生成。XML 输出不仅比纯文本和 HTML 输出更稳定，而且更容易自动解析，从而使 Biopython 更加稳定。
您可以通过多种方式获取 XML 格式的 BLAST 输出。对于解析器来说，输出是如何生成的并不重要，只要它是 XML 格式即可。
您可以使用 Biopython 通过互联网运行 BLAST，如第 通过互联网运行 BLAST 节所述。
您可以使用 Biopython 在本地运行 BLAST，如第 在本地运行 BLAST 节所述。
您可以在 NCBI 网站上通过网络浏览器自己执行 BLAST 搜索，然后保存结果。您需要选择 XML 作为接收结果的格式，并将最终获得的 BLAST 页面（您知道的，就是包含所有有趣结果的那个页面！）保存到文件中。
您还可以不使用 Biopython 在本地运行 BLAST，并将输出保存到文件中。同样，您需要选择 XML 作为接收结果的格式。
重要的是，你不必使用 Biopython 脚本获取数据才能解析它。通过以下方式之一进行操作后，你需要获取结果的句柄。在 Python 中，句柄只是一个很好的通用术语，用于描述任何信息源的输入，以便可以使用 read() 和 readline() 函数（见第 什么是句柄？ 节）检索信息。
如果你按照上面代码通过脚本与 BLAST 交互，那么你已经有了 result_handle，它是 BLAST 结果的句柄。例如，使用 GI 号进行在线搜索
>>> from Bio.Blast import NCBIWWW
>>> result_handle = NCBIWWW.qblast("blastn", "nt", "8332116")
如果你以其他方式运行 BLAST，并且在文件 my_blast.xml 中拥有 BLAST 输出（XML 格式），那么你只需打开该文件进行读取即可
>>> result_handle = open("my_blast.xml")
现在我们有了句柄，就可以解析输出。解析它的代码非常小。如果你期望单个 BLAST 结果（即，你使用了一个查询）
>>> from Bio.Blast import NCBIXML
>>> blast_record = NCBIXML.read(result_handle)
或者，如果您有大量结果（即，多个查询序列）
>>> from Bio.Blast import NCBIXML
>>> blast_records = NCBIXML.parse(result_handle)
就像 Bio.SeqIO 和 Bio.Align（见第 序列输入/输出 和 序列比对 章），我们有一对输入函数，read 和 parse，其中 read 用于你只有一个对象的情况，而 parse 是一个迭代器，用于你可能有很多对象的情况——但我们获得的不是 SeqRecord 或 MultipleSeqAlignment 对象，而是 BLAST 记录对象。
为了能够处理 BLAST 文件可能很大，包含数千个结果的情况，NCBIXML.parse() 返回一个迭代器。简单来说，迭代器允许你逐步遍历 BLAST 输出，逐一检索每个 BLAST 搜索结果的 BLAST 记录
>>> from Bio.Blast import NCBIXML
>>> blast_records = NCBIXML.parse(result_handle)
>>> blast_record = next(blast_records)
# ... do something with blast_record
>>> blast_record = next(blast_records)
# ... do something with blast_record
>>> blast_record = next(blast_records)
# ... do something with blast_record
>>> blast_record = next(blast_records)
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
StopIteration
# No further records
或者，您可以使用 for 循环
>>> for blast_record in blast_records:
... pass # Do something with blast_record
...
但请注意，你只能遍历 BLAST 记录一次。通常，你会从每个 BLAST 记录中保存你感兴趣的信息。如果你想保存所有返回的 BLAST 记录，你可以将迭代器转换为列表
>>> blast_records = list(blast_records)
现在，你可以使用索引像往常一样访问列表中的每个 BLAST 记录。但是，如果你的 BLAST 文件很大，那么尝试将所有记录保存到列表中可能会遇到内存问题。
通常，您一次只运行一个 BLAST 搜索。然后，您只需从 blast_records 中取出第一个（也是唯一一个）BLAST 记录
>>> from Bio.Blast import NCBIXML
>>> blast_records = NCBIXML.parse(result_handle)
>>> blast_record = next(blast_records)
或者更简洁地
>>> from Bio.Blast import NCBIXML
>>> blast_record = NCBIXML.read(result_handle)
我想您现在想知道 BLAST 记录中包含什么内容。
BLAST 记录类
BLAST 记录包含你可能想从 BLAST 输出中提取的所有内容。现在，我们只展示如何从 BLAST 报告中提取一些信息的示例，但如果你想要这里没有描述的特定内容，请详细查看有关记录类的信息，并查看代码或自动生成的文档——文档字符串包含大量关于每个信息片段中存储内容的信息。
为了继续我们的示例，让我们只打印出 BLAST 报告中所有高于特定阈值的命中信息的摘要信息。以下代码执行此操作
>>> E_VALUE_THRESH = 0.04
>>> for alignment in blast_record.alignments:
... for hsp in alignment.hsps:
... if hsp.expect < E_VALUE_THRESH:
... print("****Alignment****")
... print("sequence:", alignment.title)
... print("length:", alignment.length)
... print("e value:", hsp.expect)
... print(hsp.query[0:75] + "...")
... print(hsp.match[0:75] + "...")
... print(hsp.sbjct[0:75] + "...")
...
这将打印出以下摘要报告
****Alignment****
sequence: >gb|AF283004.1|AF283004 Arabidopsis thaliana cold acclimation protein WCOR413-like protein
alpha form mRNA, complete cds
length: 783
e value: 0.034
tacttgttgatattggatcgaacaaactggagaaccaacatgctcacgtcacttttagtcccttacatattcctc...
||||||||| | ||||||||||| || |||| || || |||||||| |||||| | | |||||||| ||| ||...
tacttgttggtgttggatcgaaccaattggaagacgaatatgctcacatcacttctcattccttacatcttcttc...
基本上，一旦你解析了 BLAST 报告中的信息，你可以对它做任何你想做的事情。这当然取决于你想要用它做什么，但希望这能帮助你开始做你需要做的事情！
从 BLAST 报告中提取信息的一个重要考虑因素是信息存储在何种类型的对象中。在 Biopython 中，解析器返回 Record 对象，具体取决于你正在解析的内容，可以是 Blast 或 PSIBlast。这些对象在 Bio.Blast.Record 中定义，并且非常完整。
图 表示 BLAST 报告的 Blast 记录类的类图。 和 PSIBlast 记录类的类图。 是我对 Blast 和 PSIBlast 记录类的 UML 类图的尝试。PSIBlast 记录对象类似，但支持 PSIBlast 迭代步骤中使用的轮次。
图 1 表示 BLAST 报告的 Blast 记录类的类图。
图 2 PSIBlast 记录类的类图。
如果你精通 UML 并看到可以改进的地方，请告诉我。
处理 PSI-BLAST
你可以直接从命令行或使用 python 的 subprocess 模块运行 PSI-BLAST 的独立版本（传统 NCBI 命令行工具 blastpgp 或其替代工具 psiblast）。
在撰写本文时，NCBI 似乎不支持通过互联网运行 PSI-BLAST 搜索的工具。
请注意，Bio.Blast.NCBIXML 解析器可以读取当前版本的 PSI-BLAST 的 XML 输出，但诸如每次迭代中哪些序列是新的或被重用的信息不在 XML 文件中。
处理 RPS-BLAST
你可以直接从命令行或使用 python 的 subprocess 模块运行 RPS-BLAST 的独立版本（传统 NCBI 命令行工具 rpsblast 或其同名替代工具）。
在撰写本文时，NCBI 似乎不支持通过互联网运行 RPS-BLAST 搜索的工具。
你可以使用 Bio.Blast.NCBIXML 解析器读取当前版本的 RPS-BLAST 的 XML 输出。