2. BLAST基础：序列比对原理与BLAST算法核心机制

在我们讨论GPU加速之前，我觉得有必要先把BLAST这个东西到底在干什么说清楚。很多初学者一上来就奔着调参和加速去了，结果连BLAST的算法骨架都没搞明白，遇到实际数据就抓瞎。这章我带你从序列比对原理出发，把BLAST的核心机制拆开来看。

2.1 序列比对：你到底在比什么？

序列比对说白了，就是看两条DNA或蛋白质序列有多像。你想想看，一条新测出来的序列扔到库里，怎么知道它是什么基因？怎么判断它和已知序列的亲缘关系？本质上就是在找“相似的区域”。

但这里面有个坑——相似不是只有完全匹配才算数。生物序列在演化过程中会变异，有的位置碱基换了（替换），有的多了（插入），有的少了（删除）。所以比对这活儿，本质是在容忍一定程度的错配的前提下，把两条序列对齐。

2.2 比对的数学模型：打分与回溯

怎么衡量“对齐得好不好”？业界统一的做法是打分。每条序列看成一个字符序列，比对结果就是给每对位置打个分：

• 匹配（Match）：两个字符相同，加分。
• 错配（Mismatch）：两个字符不同，减分。
• 空位（Gap）：一个序列比另一个多出字符，引入空位扣分。空位通常还分“开空位（Gap Open）”和“延伸空位（Gap Extend）”，开空位惩罚重，延伸轻一些。因为生物学上，连续插入几个碱基比频繁地插删更常见。

打分矩阵决定了具体的分值。蛋白质比对常用BLOSUM62或PAM250，DNA比对一般匹配+5、错配-4（或其派生变体）。我个人习惯在项目里先跑一次默认矩阵看效果，再根据数据特点微调空位罚分。

产生最高得分的算法，是经典的动态规划（Needleman-Wunsch做全局比对，Smith-Waterman做局部比对）。但动态规划要想做全局最优，复杂度是O(n*m)。一条查询序列几百万碱基，数据库一两个T的规模，直接用动态规划去搜……嗯，你算一下就知道，哪怕用泰坦显卡也扛不住。所以工业界不这么干。

2.3 BLAST的突破：不求最优，只求够用够快

BLAST（Basic Local Alignment Search Tool）在1990年由Altschul等人提出，它的核心思想是放弃全局最优，用启发式策略快速定位相似区域。说白了就是——我不跟你把整个比对图塞满，只挑可能有戏的局部玩一把。

BLAST的算法可以概括成三个步骤，我习惯叫它“seed-and-extend”策略：

1. Seeding（种子筛选）：把查询序列切成很多小片段（一般是3个氨基酸或11个碱基），这些小片段就是“种子”。去数据库里建一个哈希索引，查所有出现过的种子。匹配上的位置，都是候选的“高相似区域”。
2. Extending（延伸）：在每个种子匹配位置，向两边不断扩展比对，直到得分下降到阈值以下。这次扩展不要求全局最优，只算局部动态规划，保证快。
3. E-value评估（统计显著性检验）：扩展出来的比对片段（High-scoring Segment Pair, HSP）算一个E-value，即“随机出现这样高得分的概率有多大”。E-value越低，结果越可靠。通常E-value < 0.05才被认为是有统计意义的。

我最初读到这个设计时觉得它很“糙”，后来自己在项目中调优基因注释流程才意识到——大部分序列就是靠少量关键保守区域来鉴定的，全局动态规划算出来的结果和BLAST的top hit差别不大，但BLAST快了两个数量级。

2.4 BLAST家族：参数与变种

BLAST不是单一程序，而是一个工具族。不同的查询序列类型和数据库类型，BLAST会选不同的seed长度和打分策略：

工具名	用途	seed长度	备注
blastn	核酸序列 ≓ 核酸数据库	11 bp	默认seed，通用性好
blastp	蛋白质序列 ≓ 蛋白质数据库	3 aa	蛋白质保守性更强，seed可以更短
blastx	核酸序列 → 翻译 → 蛋白质数据库	3 aa（翻译后）	常用于物种间比对
tblastn	蛋白质序列 ≓ 核酸翻译库	3 aa	查找编码区时很有效
megablast	高相似度的核酸快速比对	16-28 bp	seed更长，避免短seed干扰

2.5 BLAST的编译与命令行基础

虽然本课程重点是GPU加速的BLAST，但命令行原版BLAST的基本用法你得先过一遍。下面是最常见的一组操作：

# 建索引（格式化成数据库）
makeblastdb -in reference.fasta -dbtype nucl -out ref_db

# 用blastn运行比对
blastn -query query.fasta -db ref_db -outfmt 6 -out results.tsv -evalue 1e-5 -num_threads 8

参数解释：

• -outfmt 6输出tab分隔的标准表格，方便下游分析。
• -evalue 1e-5控制显著性阈值，值越小越严格。
• -num_threads 8用8个CPU线程并行，这是CPU版BLAST唯一能利用多核的地方。

对比此时你会看到，BLAST在CPU上的并行粒度仅仅是线程级的，且数据库数据的读取和k-mer索引查询都跑在单机内存里，瓶颈很快会落在I/O和索引搜索上。

2.6 BLAST的困境：为什么需要GPU加速？

我们算一笔账：常规的blastn查询一条5000bp的序列在小规模数据库（100万个序列，2Gb）上大约需要几秒钟。但是如果你的查询序列有几十万条（比如宏基因组的组装bin），或者数据库是完整的NCBI NT库（50+Tb），CPU版本的BLAST基本要跑上数天甚至数周。而高吞吐测序数据一两天就产出几十条这样的任务。

这就是我关注GPU加速的初衷。BLAST的seed-and-extend流程中存在几个高度可并行的计算热点：

• 种子匹配阶段的hash索引搜索——单个查询的每个种子可以独立并行走。
• 延伸阶段的无依赖计算——每个seed的HSP不影响别的seed的HSP。
• 打分矩阵的访问——GPU的共享内存可以让多条线程同时读取同一评分矩阵。

把这三个环节抓到GPU上，理论上可以获得几十甚至上百倍的加速。这就是为什么我们在这个课程里要自己动手做GPU版BLAST。

2.7 BLAST核心流程结构图

嗯，这张图你多看两眼。整个BLAST就干了这三件事——其他都是细节参数和I/O包装。我们后续GPU加速的重点，就放在Step 1的哈希查询和Step 2的延伸计算上。

2.8 本章小结

序列比对本质上是个最优化打分问题，而BLAST用“先粗筛再细比对”的策略绕开了全局动态规划的高计算量。我建议新手在接触GPU加速BLAST之前，至少手动跑一次blastn，看一眼-outfmt 6的输出长什么样，理解默认E-value阈值怎么影响结果条数。概念清楚了，后面写CUDA代码时才不会迷失在索引和寄存器里。