食谱

检索不匹配的 BLAST 查询

问题

NCBI 的独立 BLAST 程序的 XML 输出不包含目标数据库中“无命中”的查询序列信息。有时你想知道哪些序列没有匹配到数据库，并相应地进一步分析/注释它们。有多种方法可以做到这一点，一种方法是使用 SeqIO 的方法 .index() 将查询文件转换为字典，然后解析结果文件以获取与字典匹配的序列。然后，可以使用 Python 的 set() 算术创建一个查询文件中但结果文件中不存在的序列列表，这些列表可以用作键来检索每个“无命中”查询的完整 SeqRecord。懂了吗？嗯，也许做起来更容易

解决方案

假设你设置了一个 BLAST 运行，其中一个名为 queries.fasta 的文件中的序列被搜索到数据库中，结果保存在 BLAST_RESULTS.xml 中

from Bio import SeqIO
from Bio.Blast import NCBIXML

# Build an index, but we don't need to parse the record
q_dict = SeqIO.index("queries.fasta", "fasta")

hits = []
for record in NCBIXML.parse(open("BLAST_RESULTS.xml")):
    # As of BLAST 2.2.19 the xml output for multiple-query blast searches
    # skips queries with no hits so we could just append the ID of every blast
    # record to our 'hit list'. It's possible that the NCBI will change this
    # behaviour in the future so let's make the appending conditional on there
    # being some hits (ie, a non-empty alignments list) recorded in the blast
    # record

    if record.alignments:
        # The blast record's 'query' contains the sequences description as a
        # string. We used the ID as the key in our dictionary so we'll need to
        # split the string and get the first field to remove the right entries
        hits.append(record.query.split()[0])

misses = set(q_dict.keys()) - set(hits)
orphan_records = [q_dict[name] for name in misses]

我们可以进行一个小小的健全性检查以确保一切正常

>>> print(
...     "found %i records in query, %i have hits, making %i misses"
...     % (len(q_dict), len(hits), len(misses))
... )
found 11955 records in query, 2802 have hits, making 9153 misses

很好，现在你有了一个包含所有“非命中”序列的列表 (orphan_records)，其中包含 SeqRecord 对象，你可以根据需要对其进行注释/分析，或者使用 SeqIO.write() 来创建一个仅包含这些序列的新文件，这些文件可以被输入到另一个程序中。

讨论

如上所述，大多数情况下，每次运行中都会花费大量时间填充 BLAST 解析器中的命中列表，如果对结果文件中的每个记录进行一次检查以查看它是否在字典中，这是否是一种更节省内存的方法，尤其是在处理非常大的文件时？