📚 课程目录第1章 BLAST简介与应用场景

1. BLAST 简介与应用场景

1.1 什么是 BLAST？

BLAST，全称 Basic Local Alignment Search Tool。说直白点，它就是做序列比对的。

给你的 DNA 或蛋白质序列找个亲戚。看看它在数据库里跟谁长得像。我早年写序列比对算法时，BLAST 是绕不开的基准。它的核心思想其实很简单——用短种子快速命中，再延伸比对。速度比老式的 Needleman-Wunsch 快几个数量级。

你可能会问：为什么不用 Smith-Waterman 做局部比对？嗯，问题在于速度。Smith-Waterman 精确但太慢。BLAST 牺牲了一点精确度，换来了实用的大规模搜索能力。我在项目中处理过几万条转录组序列，如果不用 BLAST，光比对就得跑一周。

1.2 BLAST 能解决什么问题？

说白了，BLAST 帮你回答三个问题：

这条序列是什么物种的？
它跟已知基因同源吗？
它的功能可能是什么？

但我今天要讲的，不是怎么调 BLAST 参数。我想说一个更重要的东西——BLAST 结果如何变成机器学习的输入特征。

这几年我做生物信息学与机器学习交叉的项目，发现一个大坑：很多人在跑完 BLAST 后，直接把 top 1 命中当答案，然后就丢给模型训练。这是巨大的信息浪费。BLAST 输出里的每个字段——e-value、bitscore、覆盖度、一致度——其实都是宝贵的特征。

个人经验：我曾在一个宏基因组项目里，光靠 BLAST 字段的组合特征，就把物种分类准确率从 82% 提到了 94%。把 BLAST 当成黑盒是浪费，你得把它当成特征生成器。

1.3 BLAST 的三种常见用法

我习惯把 BLAST 的应用场景分成三类。下表总结了一下：

应用类型	查询序列	目标数据库	典型目标
物种鉴定	未知序列	NT / NR	找到最近亲缘物种
功能注释	预测蛋白序列	Swiss-Prot / Pfam	推断分子功能
同源搜索	感兴趣基因	专用基因组库	发现新同源物

你看，不同场景下，我们需要关注的 BLAST 字段是不一样的。物种鉴定更看重覆盖度和一致度；功能注释更关注保守结构域的 hits 分布；同源搜索则看重 e-value 的排序和 bitscore 的比值。

我做过一个抗性基因挖掘的项目。当时用 BLAST 比对完，发现一个问题：单纯靠 e-value 排序，会把一些短序列的假阳性排在前面。后来我改用 bitscore 除以序列长度作为归一化特征，才稳定下来。

注意：e-value 是依赖数据库大小的。同一比对结果，数据库变大，e-value 会变差（变大）。所以跨数据库比较时，优先用 bitscore 而不是 e-value。

1.4 BLAST 输出的核心字段

当你跑 BLAST 时，输出格式可以配成表格模式的 -outfmt 6。我常用的是这个：

# BLAST 默认的 outfmt 6 字段顺序
1. qseqid     — 查询序列 ID
2. sseqid     — 目标序列 ID
3. pident     — 一致百分比
4. length     — 比对长度
5. mismatch   — 错配数
6. gapopen    — 缺口数
7. qstart     — 查询起始位
8. qend       — 查询结束位
9. sstart     — 目标起始位
10. send      — 目标结束位
11. evalue    — e 值
12. bitscore  — 比特分数

这些字段就是原始特征池。我个人习惯会用它们做衍生特征，比如：

覆盖度 = (qend - qstart + 1) / 查询序列全长
归一化 bitscore = bitscore / 比对长度
gap 密度 = gapopen / length
top hits 数量：同一个 query 有多少个显著 hits

你别小看这些衍生特征。我在一个酶活性预测模型里，只用 pident 和 evalue 两个原始特征，AUC 只有 0.71。加上 覆盖度 和 归一化 bitscore 后，AUC 直接跳到了 0.85。

1.5 BLAST 结果特征构建的整体流程

我们把这个过程画成一张图，你看一眼就明白了：

流程不复杂。关键在“特征抽取”那个环节。很多人跳过了这一步，直接用 BLAST 的输出字段扔给模型——这等于把原始矿石直接拿去冶炼，杂质太多。

1.6 常见应用场景举例

讲几个我用过的实际场景，方便你对号入座：

宏基因组物种分类：把每条 contig 的 BLAST 结果做成 20 维特征向量，再用随机森林分类器区分属级别。我试过，准确率比基于 k-mer 的方法高 10 个百分点。
抗菌肽活性预测：抗菌肽很短（10-50 aa），BLAST 比对通常覆盖不全。这时候我用 bitscore/query_length 做归一化特征，模型泛化能力明显提升。
酶功能亚类区分：同属同一个超家族的酶，活性位点不同。只用 e-value 基本分不开。但我把 BLAST 输出中 qstart~qend 的覆盖位置做成分布直方图特征，模型就能捕捉到保守基序的差异。

项目实例：有次做肠道微生物功能 profiling，跑完 BLAST 后我尝试用 top 10 hits 的 bitscore 分布替代简单的 top 1 e-value。结果逻辑回归的 F1 从 0.68 升到 0.83。你想想看，改的不是模型，是用什么特征。

1.7 一个简单的命令示例

不过多的讲基础，我顺手给你一个跑 BLAST 的例子。这在后面的章节会详细展开：

# 用蛋白质序列搜索 NR 数据库
blastp -query my_proteins.fa \
       -db nr \
       -out my_blast_results.tsv \
       -outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore" \
       -evalue 1e-5 \
       -num_threads 8

输出的 my_blast_results.tsv 就是你后续特征工程的原料。别忘了加上 -max_target_seqs 10 或更大的值，保证每个 query 有足够多的 hits。

1.8 从 BLAST 结果到特征矩阵：最简代码

我习惯用 Python 的 BioPython 解析 BLAST 输出。一段极简的特征提取代码：

from Bio.Blast import NCBIXML

# 解析 XML 格式的 BLAST 输出
for record in NCBIXML.parse(open("blast_results.xml")):
    query_len = record.query_length
    for align in record.alignments[:10]:   # 只取前 10 个 hits
        hsp = align.hsps[0]
        coverage = (hsp.query_end - hsp.query_start + 1) / query_len
        norm_bitscore = hsp.bitscore / hsp.align_length

        # 这就是你的特征
        features = [
            hsp.identities / hsp.align_length,  # 一致度
            hsp.expect,                          # e-value
            hsp.bitscore,                        # 原始 bitscore
            coverage,                            # 覆盖度
            norm_bitscore                        # 归一化 bitscore
        ]

这个循环跑完以后，你就有了一个特征矩阵。每行是一个 query 与其某个 hit 的组合特征。下一章我们会讲如何把这些特征整理成适合模型训练的格式。

嗯，本章先到这里。主要是让你对 BLAST 的特征潜力有个整体认识。别急，后面我们逐步拆解每一个特征的物理意义和建模价值。