使用Biopython处理Sanger测序序列后解析NCBIWWW.qblast返回XML文件问题

BLAST运行慢的原因 #

当前使用的NCBIWWW.qblast是调用NCBI公共在线BLAST服务，速度受网络质量、NCBI服务器负载、同时提交任务数量影响，确实普遍较慢。如果需要大幅提升速度，可本地安装NCBI BLAST+工具包，下载nt数据库到本地后运行本地比对，速度会提升数倍到数十倍。

现有代码优化点 #

现有代码存在2个可修正的问题，避免后续运行出错：

• 字符串strip方法使用错误：filename.strip("_download.fasta")是按单个字符匹配删除首尾内容，不是匹配完整后缀，会导致文件名包含_、d、o、w等字符时被误删。Python 3.9及以上版本可替换为filename.removesuffix("_download.fasta")，全版本兼容的写法为filename.rsplit("_download.fasta", 1)[0]，第二个后缀处理同理。
• 修剪后序列输出无需手动按50bp换行，直接调用SeqIO写入即可保证FASTA格式规范，避免手动拼接字符串出现格式错误。

BLAST XML结果解析代码 #

在现有代码末尾添加以下代码，即可提取每条序列top hit的物种及比对信息：

from Bio.Blast import NCBIXML
import glob

for xml_file in glob.iglob("*trim.xml"):
    print(f"解析文件：{xml_file}")
    with open(xml_file, "r") as handle:
        # 遍历XML中所有序列的比对结果
        for blast_rec in NCBIXML.parse(handle):
            # 处理无比对结果的情况
            if not blast_rec.alignments:
                print(f"序列{blast_rec.query}：无匹配结果")
                print("-"*40)
                continue
            # 取top1比对结果
            top_align = blast_rec.alignments[0]
            top_hsp = top_align.hsps[0]
            # 提取物种信息（NCBI返回的hit标题中物种名包含在方括号内）
            if "[" in top_align.title and "]" in top_align.title:
                species = top_align.title.split("[")[-1].split("]")[0].strip()
            else:
                species = "物种信息未标注"
            # 输出结果
            print(f"序列ID：{blast_rec.query}")
            print(f"匹配物种：{species}")
            print(f"序列相似度：{top_hsp.identities / top_hsp.align_length * 100:.2f}%")
            print(f"比对长度：{top_hsp.align_length}")
            print("-"*40)

内容的提问来源于stack exchange，提问作者allisonrs