📚 课程目录第4章提取BLAST特征

4. 提取BLAST特征：score归纳总结

做机器学习特征，最怕的就是原始特征太粗糙。

BLAST 输出的 score 本身就是一个极好的判别信号。但直接把它扔进模型，效果往往一般。你想想看，一个 score 是 100 的比对，和一个 score 是 85 的，它们真的只差了 15 分么？不一定。

我刚开始做 BLAST 特征提取时，也犯了懒——直接把 raw score 丢进去。结果模型在验证集上表现奇差，调参调了一周也没救。后来我重新审视数据，发现 score 的分布藏着很多信息。嗯，当时我就意识到，真正有用的不是单个 score 值，而是它在全局中的 "相对位置" 和 "统计特性"。

score 的 "骨架"：我们到底要提取什么？

BLAST 输出里，每个 HSP 都会带一串 score 相关的字段。但我们要归纳的，不仅是那个 "最高分"，而是一整套描述比对质量的统计指标。

我个人习惯把 score 特征分为三个层次：

原始值：bit score、raw score
归一化值：对长度、数据库规模做修正后的 score
衍生统计量：均值、方差、分位数、极差、偏度

  核心原则：不要只取 "最佳 hit" 的 score，而要抽取所有 hits 的 score 分布特征。这样做能捕捉到序列在数据库中的 "整体相似度轮廓"。

例如，一个 query 比对到 100 个 subject，score 集中在 40-50 之间；另一个 query 比对到同样多 subject，score 从 30 跨越到 200。它们代表的生物学意义完全不同。前者可能是一个保守结构域，后者可能是同一超家族内的远源同源——特征向量必须能编码这种差异。

常见的 score 特征及其计算方式

下面这张表，是我在项目中总结的一套 "score 归纳特征模板"。每次做 BLAST 特征工程，我基本都会先跑一遍。

特征名	计算方式	用途说明
max_bit_score	所有 HSP 中 bit score 的最大值	最强比对信号，区分度高
mean_bit_score	所有 HSP 的 bit score 算术平均	反映整体比对强度
std_bit_score	bit score 的标准差	衡量得分离散程度，过大可能意味着数据库污染
score_range	max_bit_score - min_bit_score	极端值跨度
median_bit_score	bit score 的中位数	抗噪性强，适合有长尾分布的数据
log_evalue_neg_min	-log10(最小的 e-value)	将 e-value 转化为正向特征，便于模型训练
hit_count	e-value < 阈值（如 0.001）的 hit 数量	近似同源序列丰富度
top5_score_slope	前 5 个 hit score 随 rank 的衰减斜率	描述比对强度随排名下降的速度

我的经验：不要忽视 "top5_score_slope"。在很多实际案例中，斜率陡峭的 query 往往是高度特异的单拷贝基因；而斜率平缓的 query 通常属于多基因家族或含有重复结构域——这个特征在分类任务中非常有用。

实践流程：从 BLAST 输出到 score 特征向量

我一般分四步走。每一步都很简单，但顺序不能乱。

解析 BLAST 输出：用 BioPython 的 NCBIXML.parse() 或自己写 regex 拉取所有 HSP 的 score 字段
构建 score 列表：对每个 query，将其所有 hits 的 bit score 收集到一个数组
计算统计量：用 numpy 或 pandas 快速计算上面表格中的特征
合并为特征向量：将所有统计量按固定顺序拼接成一行特征

我曾经在一个跨物种比较项目里，忘了做第 2 步——直接取了第一个 HSP 的 score。结果模型在某个目上表现极差，调了三天才发现，因为这个目的序列普遍存在多域结构，第一个 HSP 根本不代表整体。所以，一定要确保你收集的是全部 hits 的 score 分布，而不是被默认排序欺骗了。

代码实战：从 BLAST XML 提取 score 统计量

下面这段代码，是我自己项目里一直在用的模板。它读入一个 BLAST XML 文件，输出一个包含 8 个 score 特征的 DataFrame。

from Bio.Blast import NCBIXML
import numpy as np
import pandas as pd

def extract_score_features(xml_path):
    records = []
    with open(xml_path) as handle:
        for record in NCBIXML.parse(handle):
            scores = []
            for alignment in record.alignments:
                for hsp in alignment.hsps:
                    scores.append(hsp.bits)
            # 如果没有 hit，填充 NaN
            if not scores:
                records.append({
                    'query_id': record.query,
                    'max_bit_score': np.nan,
                    'mean_bit_score': np.nan,
                    'std_bit_score': np.nan,
                    'score_range': np.nan,
                    'median_bit_score': np.nan,
                    'log_evalue_neg_min': np.nan,
                    'hit_count': 0,
                    'top5_score_slope': np.nan
                })
                continue
            scores = np.array(scores)
            # 降序排列
            scores_sorted = np.sort(scores)[::-1]
            # 计算 e-value 正向化
            evalues = []
            for alignment in record.alignments:
                for hsp in alignment.hsps:
                    # 避免 log 0
                    evalues.append(max(hsp.expect, 1e-180))
            log_evalue = -np.log10(np.min(evalues)) if evalues else np.nan
            # 计算 top5 斜率
            if len(scores_sorted) >= 5:
                top5 = scores_sorted[:5]
                x = np.array([0,1,2,3,4])
                slope, _ = np.polyfit(x, top5, 1)
            else:
                slope = np.nan
            records.append({
                'query_id': record.query,
                'max_bit_score': np.max(scores),
                'mean_bit_score': np.mean(scores),
                'std_bit_score': np.std(scores),
                'score_range': np.ptp(scores),
                'median_bit_score': np.median(scores),
                'log_evalue_neg_min': log_evalue,
                'hit_count': len(scores_sorted),
                'top5_score_slope': slope
            })
    return pd.DataFrame(records)

# 使用示例
df_scores = extract_score_features('blast_results.xml')
print(df_scores.head())

我曾经踩过的坑：早期版本我直接用了 hsp.expect 取最小值，没做 log 转换。e-value 是一个非常 skewed 的分布（从 1 到 1e-180），原值输入模型会让线性模型完全失效。记得一定要先做 -log10() 变换，才能让模型有效捕捉到信号。

更复杂的场景：多数据库与多 BLAST 输出合并

有时你会对同一 query 做多个 BLAST 搜索（比如不同数据库、不同参数），这时 score 特征需要跨输出归纳。

我个人做法是把每个输出当作一个 "视角"，然后将视角内的统计量作为独立特征组。举个例子：

对 nr 库的 BLAST 结果，提取 score 的 max、mean、std
对 Swiss-Prot 库的 BLAST 结果，提取同样的统计量
对 Pfam 库的 HMMER 结果，提取 domain score 的分布特征

最后把这些特征组水平拼接。这样做，模型能看到 "在 nr 中普遍高分但在 Swiss-Prot 中低分" 这类跨库矛盾信号——往往与假阳性或错误注释高度相关。

SVG 流程图：score 归纳的完整链路

下面这张图，展示了我刚才讲的全部流程。你边看代码边对照这张图，会更容易理解。

避坑指南：我摔过的跤

列几条常见的坑，你在做的时候一定要绕开。

坑 1：有些 BLAST 输出的 HSP 数量是 0。这种情况，scores 列表是空的。如果直接计算 np.mean(scores) 会报错。记得提前判断，用 NaN 填充。

坑 2：e-value 不要直接用原值做特征。它的分布太偏了，模型几乎学不到东西。我习惯用 -log10 变换，把它变成 "越大越显著" 的正向指标。

坑 3：做跨数据库合并时，注意 score 的尺度可能不同。不同 BLAST 程序（blastn vs blastp）甚至不同数据库版本的 bit score 计算方式有细微差别。建议在特征工程阶段先做标准化，或至少确认每个特征的分布是否可比。

小结：score 归纳是特征工程的 "地基"

score 特征提取，说白了就是把 BLAST 输出的原始比对信息，转化为模型能理解的数值信号。这个步骤没做好，后面的模型再怎么调参都是白费功夫。

我个人的原则是：宁可多算十个统计量，也不要只用一个 max。因为你永远不知道哪个特征会在你的数据上产生奇效。