Biopython调用BLAST比对16S rRNA时得到的e值全部为0.0是什么原因？

为什么e值全为0.0 #

• 这不是代码或者BLAST返回错误：16S rRNA序列长度在1500bp左右，当比对的序列同源性较高时，期望e值会远小于双精度浮点数的最小可表示正数值（约2.2e-308），Python解析时会直接截断显示为0.0，和序列是否100%匹配没有直接关系。
• 肉眼观测到的非100%匹配结果排序非常靠后，原代码没有做结果条数限制，默认返回的前数百条高同源性16S序列的e值都会低于浮点数下限，显示为0.0。

修复代码实现需求（取前5条最佳比对） #

调整后的代码如下：

from Bio import SeqIO
from Bio import Entrez
from Bio.Blast import NCBIWWW
from Bio.Blast import NCBIXML

Entrez.email = "obsto123@gmail.com"
handle = Entrez.efetch(db="nucleotide",id="NC_000913.1",rettype="gb",retmode="text")
seq_record = SeqIO.read(handle,"genbank")

records = []
for feature in seq_record.features:
    if feature.type == "rRNA" and "16S" in feature.qualifiers["product"][0]:
        records.append(seq_record.seq[feature.location.start.position:feature.location.end.position])
rRNA = records[0]

print("Now doing blast search...")
# 直接在qblast里指定返回前5条结果，减少请求和解析耗时
blast_data = NCBIWWW.qblast("blastn","nt",rRNA, hitlist_size=5)

# 保存xml文件
with open("blast_xml.xml","w") as f:
    f.write(blast_data.read())

# 解析并打印结果
with open("blast_xml.xml","r") as xml:
    blast_record = NCBIXML.read(xml)
    
for idx, alignment in enumerate(blast_record.alignments, 1):
    for hsp in alignment.hsps:
        print(f"===== 第{idx}条最佳比对 =====")
        print(f"序列名称: {alignment.title}")
        print(f"比对长度: {hsp.align_length}")
        print(f"一致性: {hsp.identities/hsp.align_length:.2%}")
        print(f"e值: {hsp.expect}")
        # 要查看真实的极小e值可以打印科学计数法格式
        print(f"e值（科学计数法）: {hsp.expect:.2e}")
        print("")

额外说明 #

• 如果需要对e值做排序或者阈值筛选，直接用0.0判断即可，所有显示为0.0的结果都属于同源性极高的可信比对结果。
• 如果需要查看非高同源的比对结果，可以在qblast参数中调整expect阈值，比如设置expect=1e-5来过滤低置信度结果，也可以调大hitlist_size返回更多结果。

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