📚 课程目录第2章 BLAST结果格式解析（tabular, XML）

BLAST结果格式解析（tabular, XML）

拿到BLAST输出结果，第一件事是什么？

我个人习惯先看格式。BLAST默认给人看的格式是带标题、对齐块、打分条的“美观”版面，但机器读起来简直要命。你想想看，几千条比对记录，每一条都带一堆排版符号，写脚本来解析能把人逼疯。

还好NCBI给了我们更聪明的选项：outfmt 6（表格格式）和outfmt 5（XML格式）。这一节我带你把这些格式解剖清楚，顺便讲讲我在实际项目中踩过的坑。

outfmt 6：那个“偷懒”的表格格式

说白了，outfmt 6就是一行一条比对记录，用制表符分隔。格式简单到令人感动。

qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore

每条记录包含12个字段，下面是标准定义：

字段名	含义	类型
qseqid	查询序列ID	字符串
sseqid	目标序列ID	字符串
pident	序列一致性百分比	浮点数
length	比对长度	整数
mismatch	错配数	整数
gapopen	缺口开放数	整数
qstart	查询序列起始位置	整数
qend	查询序列结束位置	整数
sstart	目标序列起始位置	整数
send	目标序列结束位置	整数
evalue	期望值	浮点数（科学计数）
bitscore	比特分数	浮点数

经验提醒：我刚开始做特征工程时，以为evalue越低越好，就直接把它当负向特征用。后来发现有些假阳性的evalue也很低。建议你把bitscore和pident组合起来，比单用evalue靠谱得多。

运行BLAST时加这个参数就能输出表格格式：

blastn -query query.fasta -db nt -outfmt 6 -out result.tab

想要自定义字段顺序？可以用这个：

blastn -query query.fasta -db nt -outfmt "7 qseqid sseqid evalue bitscore" -out result_custom.tab

XML格式：那个“丰满”的备选方案

表格格式虽然简洁，但信息密度上比XML差远了。XML格式把每条比对的详细信息都包进去了——包括比对中间区域的逐位匹配状态、HSP分段、框架信息。

我一般在做深度特征提取时才会用XML，比如要分析gap具体在哪些位点出现、比对区域的碱基匹配度分布之类。

典型XML结构长这样：

<?xml version="1.0"?>
<!DOCTYPE BlastOutput PUBLIC "-//NCBI//NCBI BlastOutput/EN" "http://www.ncbi.nlm.nih.gov/dtd/NCBI_BlastOutput.dtd">
<BlastOutput>
  <BlastOutput_program>blastn</BlastOutput_program>
  <BlastOutput_db>nt</BlastOutput_db>
  <BlastOutput_iterations>
    <Iteration>
      <Iteration_iter-num>1</Iteration_iter-num>
      <Iteration_query-ID>Query_1</Iteration_query-ID>
      <Iteration_query-def>Sample_sequence</Iteration_query-def>
      <Iteration_hits>
        <Hit>
          <Hit_num>1</Hit_num>
          <Hit_id>gi|12345|ref|NM_001234.1|</Hit_id>
          <Hit_def>Species gene mRNA</Hit_def>
          <Hit_len>1850</Hit_len>
          <Hit_hsps>
            <Hsp>
              <Hsp_num>1</Hsp_num>
              <Hsp_bit-score>654.3</Hsp_bit-score>
              <Hsp_score>678</Hsp_score>
              <Hsp_evalue>1.2e-45</Hsp_evalue>
              <Hsp_query-from>1</Hsp_query-from>
              <Hsp_query-to>450</Hsp_query-to>
              <Hsp_hit-from>100</Hsp_hit-from>
              <Hsp_hit-to>549</Hsp_hit-to>
              <Hsp_qseq>ATCGACGTAGCTAGCTAGCTAGC...</Hsp_qseq>
              <Hsp_hseq>ATCGACGTAGCTAGCTAGCTAGC...</Hsp_hseq>
              <Hsp_midline>|||| ||| |||||||||||||||...</Hsp_midline>
            </Hsp>
          </Hit_hsps>
        </Hit>
      </Iteration_hits>
    </Iteration>
  </BlastOutput_iterations>
</BlastOutput>

关键节点：

顶部信息：程序名、数据库、版本号
迭代循环：每条查询序列是一个Iteration节点
Hit节点：每个比中的目标序列
Hsp节点：每个比中“片段”（一条序列可以有多个比中区域）
Hsp_midline：逐位反映匹配/错配/缺口的状态

两种格式的对比

对比维度	outfmt 6（表格）	outfmt 5（XML）
文件大小	小	大（通常3~5倍）
解析难度	简单（split了事）	中等（需要XML解析器）
信息密度	仅摘要字段	包含完整比对细节
特征提取支持度	基础特征（分数、长度、位置）	高级特征（gap分布、HSP分段、序列级相似度）
适用场景	快速筛选、初筛模型	精细化特征工程

你看这个对比就明白了：如果你的模型只需要evalue、bitscore、identity这些“粗粒度”特征，outfmt 6完全够用。但如果你要分析“比对区域内匹配分布的波动性”这类高级模式，就得用XML。

我曾经在一个抗性基因鉴定的项目里，仅靠outfmt 6的特征死活区分不开两个高度相似的家族。切换到XML之后，我提取了midline连续匹配片段的长度方差——竟然发现这是一个关键区分特征。有时候细节真的藏在全量信息里。

解析代码示例

先看解析outfmt 6有多简单：

import pandas as pd

columns = [
    "qseqid", "sseqid", "pident", "length",
    "mismatch", "gapopen", "qstart", "qend",
    "sstart", "send", "evalue", "bitscore"
]

df = pd.read_csv("result.tab", sep="\t", names=columns, header=None)
print(df.head())

搞定。一行不出错。

XML解析稍微复杂一点，我用ElementTree来实现：

import xml.etree.ElementTree as ET

tree = ET.parse("result.xml")
root = tree.getroot()

hits = []
for iteration in root.iter("Iteration"):
    query_id = iteration.find("Iteration_query-ID").text
    for hit in iteration.iter("Hit"):
        hit_id = hit.find("Hit_id").text
        for hsp in hit.iter("Hsp"):
            evalue = hsp.find("Hsp_evalue").text
            bitscore = hsp.find("Hsp_bit-score").text
            hits.append({
                "query_id": query_id,
                "hit_id": hit_id,
                "evalue": float(evalue),
                "bitscore": float(bitscore)
            })

print(f"共解析出 {len(hits)} 条HSP记录")

注意：XML文件如果特别大（比如上G级别的比对结果），直接全量解析到内存可能会导致OOM。我个人建议用xml.etree.ElementTree.iterparse做流式解析，逐条释放节点。我在处理全基因组比对时就遇到过这个问题，内存爆了两次才学乖。

context = ET.iterparse("large_result.xml", events=("end",))
for event, elem in context:
    if elem.tag == "Hsp":
        bitscore = elem.find("Hsp_bit-score").text
        evalue = elem.find("Hsp_evalue").text
        # 处理该HSP
        elem.clear()  # 释放内存，重要！

解析流程知识图谱

把上面讲的内容理一理，核心逻辑其实就是一条流水线：从原始输出到特征能用，中间需要过滤、解析、精炼三步。

这个图的核心信息：

BLAST输出可以有两个流向：表格和XML
表格格式解析快，适合大规模初筛
XML格式信息多，适合精细化特征工程
最终都汇入特征提取阶段，但提取的“原料”丰富度不同

决定之前，想清楚一个问题

到底选outfmt 6还是outfmt 5？我个人判断标准很简单：

选outfmt 6：数据量超过10万条比对、特征维度在20维以内、需要快速出基线模型
选outfmt 5：模型精度卡在瓶颈、需要从比对细节里挖掘新特征、数据量在万级以下（或有大内存服务器）

好的数据集决定模型上限，特征工程只是逼近这个上限。解析格式这一步，其实就是决定你“能看到多少信息”的起点。