Biopython Blast模块指定物种查询同源DNA序列无结果求助

需要在弗格森埃希氏菌（E. fergusoni，taxid:564）中筛选与给定大肠杆菌编码序列同源、全片段相似度在80%-98%之间的DNA序列，用于构建大肠杆菌基因进化树外类群。在线NCBI Blast指定该物种可正常返回结果，但批量处理大量基因时，使用Biopython的NCBIWWW.qblast函数通过entrez_query参数限定目标物种，尝试多种语法均无结果且请求耗时长达20分钟，需寻求正确的查询语法及优化方案。

运行环境：Python 3.11、Biopython 1.83、VS Code 1.89.1（Jupyter Notebook）

测试代码 #

Blast.email = my_email@email.com
from Bio.Blast import NCBIWWW
from Bio.Blast import NCBIXML

sequence_data = open("gene0.fasta").read()
result_handle = NCBIWWW.qblast(
    program="blastn",
    database="nt",
    sequence=sequence_data,
    # 尝试过的查询语法
    # entrez_query = "txid564[Organism]",
    # entrez_query = "Escherichia fergusonii (taxid:564)[Organism]",
    # entrez_query = "txid564[ORGN]",
    # entrez_query = "564[Taxid]",
    format_type="XML",
)

with open("gene0_results.xml", "w") as out_handle:
    out_handle.write(result_handle.read())

# 结果筛选逻辑
for record in NCBIXML.parse(open("gene0_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.identities/len(hsp.query) < 0.98) and (hsp.identities/len(hsp.query) >= 0.8)  : 
                    print("match: %s " % align.title[:100])
                    print(f"length: {align.length}")
                    print(f"e value: {hsp.expect}")
                    print(f"identity : {hsp.identities/len(hsp.query)} ({hsp.identities} nucleotides over {len(hsp.query)})")
                    print(hsp.query[0:75] + "...")
                    print(hsp.match[0:75] + "...")
                    print(hsp.sbjct[0:75] + "...")

注：未指定entrez_query时，1分钟内即可返回100%相似度的大肠杆菌序列；指定后无结果且请求超时。

1. 正确的entrez_query语法 #

"txid564[ORGN]"是NCBI Entrez查询物种的标准有效语法，之前尝试该语法无结果可能是因为：

• 未开启megablast模式，普通blastn对近缘物种搜索效率低，且易遗漏匹配序列
• 数据库返回结果过多导致超时，需限制返回数量

若需进一步缩小范围，可使用组合筛选语法：

entrez_query="txid564[ORGN] AND refseq[filter]"

该语法仅从RefSeq数据库中筛选弗格森埃希氏菌序列，减少非目标序列干扰。

2. 请求性能优化 #

针对批量处理场景，需调整qblast参数提升效率：

• 添加megablast=True：针对近缘物种同源序列，megablast模式比普通blastn速度更快，更适合高相似度匹配需求
• 设置hitlist_size=N：限制返回结果数量（如hitlist_size=50），减少数据传输量
• 批量处理时添加请求间隔：NCBI对API请求频率有限制，每两次请求间隔至少1秒，避免被限流导致超时

优化后的qblast调用代码：

result_handle = NCBIWWW.qblast(
    program="blastn",
    database="nt",
    sequence=sequence_data,
    entrez_query="txid564[ORGN]",
    megablast=True,
    hitlist_size=50,
    format_type="XML",
)

3. 结果筛选逻辑优化 #

原代码基于单个HSP的相似度计算可能不准确，建议改为基于整个比对序列的整体相似度：

# 替换原结果筛选循环
for record in NCBIXML.parse(open("gene0_results.xml")): 
    if record.alignments: 
        print("\n") 
        print("query: %s" % record.query[:100]) 
        for align in record.alignments: 
            total_identity = 0
            total_query_length = 0
            for hsp in align.hsps:
                total_identity += hsp.identities
                total_query_length += len(hsp.query)
            if total_query_length == 0:
                continue
            overall_identity = total_identity / total_query_length
            if 0.8 <= overall_identity < 0.98:
                print("match: %s " % align.title[:100])
                print(f"length: {align.length}")
                print(f"overall identity: {overall_identity:.2f}")
                # 可选：输出第一个HSP的比对片段
                if align.hsps:
                    hsp = align.hsps[0]
                    print(hsp.query[0:75] + "...")
                    print(hsp.match[0:75] + "...")
                    print(hsp.sbjct[0:75] + "...")

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