食譜

檢索不匹配的 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 解析器填充命中列表，如果遇到非常大的檔案，每次將結果檔案中的每個記錄與字典逐一檢查，是否會是一種記憶體使用較少的方法？