1. BLAST基础：算法原理与常用变体
各位好，我是老张。咱们开始第一讲，先把BLAST这个老朋友重新认识一下。很多做生信的同学天天用BLAST，但真问到它背后怎么跑的，有时候还真说不清楚。
BLAST的全称是Basic Local Alignment Search Tool，说白了就是“局部比对搜索工具”。它不是做全局比对，而是找局部相似片段。为什么要这样？因为两条同源序列可能只在某个结构域高度保守，其他部分早就漂变到妈都不认得了。全局比对在这里意义不大，局部比对才是正道。
1.1 BLAST的核心思想
我当年刚接触BLAST时有个困惑：为什么它比全库暴力Smith-Waterman快那么多？直到我看了算法原文，才明白它用了启发式策略——用命中种子词加速搜索。
核心流程分三步走：
- 种子词索引：把查询序列切成固定长度的小片段，比如3个氨基酸或11个碱基；同时把目标数据库也建成索引。
- 命中筛选：用查询种子在数据库中找完全或近似匹配的种子位点。
- 延伸与打分：对命中区域向两侧延伸，直到打分下降到阈值以下，输出高分片段对。
嗯，这里要注意：BLAST不会遍历所有可能性，而是先通过种子快速锁定候选区域。我对它的评价是——用可接受的敏感性损失，换来了几个数量级的加速。
关键数字：BLAST在常用数据库上（如nr、nt）的搜索速度通常是Smith-Waterman算法的10到50倍，而敏感度损失在5%以内。对于长序列比对，这个比例还会更优。
1.2 常用BLAST变体
BLAST不是一套工具，而是一族。不同的输入输出形式，对应不同的变体。我个人习惯按查询序列类型和目标数据库类型来选，这样不容易出错。
| 变体名称 | 查询序列 | 目标数据库 | 典型应用场景 |
|---|---|---|---|
| BLASTn | 核苷酸 | 核苷酸 | 核酸序列比对、引物设计验证、基因组注释 |
| BLASTp | 蛋白质 | 蛋白质 | 同源蛋白搜索、功能预测、结构域识别 |
| BLASTx | 核苷酸（翻译） | 蛋白质 | 新物种转录组鉴定、病毒序列功能标注 |
| tBLASTn | 蛋白质 | 核苷酸（翻译） | 在基因组中寻找未注释的编码区 |
| tBLASTx | 核苷酸（翻译） | 核苷酸（翻译） | 跨物种保守非编码区发现 |
我的选型建议：如果是已知基因的相似性搜索，直接上BLASTn，速度快。如果你是拿着一段cDNA序列要做功能预测，优先考虑BLASTx，因为它把翻译后的蛋白相似性映射回来，敏感性更高。
1.3 关键参数与算法细节
我记得有一次帮团队做大规模宏基因组注释，调参调了一个星期才找到合适的配置。BLAST的参数远比你想象的多，但几个核心参数必须懂：
- word size（种子长度）：BLASTn默认11，BLASTp默认3。这个是敏感性和速度的平衡杆。短种子更敏感，但慢；长种子快，但容易漏。我曾经在寻找高度变异病毒株时，把word size从11降到7，才把同源序列抓回来。
- e-value（期望值）：默认10。它表示在当前数据库大小下，随机得到当前或更好打分的次数。e-value越小越严格。1e-5是很多论文的默认分界点。
- gap penalty（空位罚分）：BLASTp默认存在空位罚-11，延伸罚-1。这个组合对绝大多数蛋白质结构域识别是合理的。
避坑指南：我遇到过有人拿着BLASTx跑一群细菌16S rRNA序列——这其实完全没有必要，因为16S是核酸序列比对，BLASTn更直接更快。更严重的是，BLASTx会把所有六帧翻译都跑一遍，数据库大小其实变成了6倍，输出结果也多了很多没有生物学意义的匹配。选错工具不仅浪费时间，还会引入噪声。
1.4 实操：命令行示例
说了这么多，咱们动手看看。BLAST+工具包安装基本已成标配，核心命令就是blastn、blastp、blastx这几个。
# BLASTn 示例：用核酸序列查询核酸数据库
blastn -query query.fasta -db nt -out result_blastn.txt -outfmt 6 -evalue 1e-5
# BLASTp 示例：用蛋白序列查询蛋白数据库
blastp -query query_prot.fasta -db nr -out result_blastp.txt -outfmt 6
# BLASTx 示例：用核酸序列查询蛋白数据库（自动六帧翻译）
blastx -query query_nuc.fasta -db nr -out result_blastx.txt -outfmt 6 -evalue 1e-3这里-outfmt 6 是制表符分隔的简洁输出，列顺序是：qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore。每条记录就是一条HSP。
小技巧：如果你的查询序列只有几百bp，建议加上 -max_target_seqs 20，避免输出几千条无用的低分匹配。我以前在服务器上被同事的BLAST作业弄炸过硬盘，就是因为没有限制输出数量。
1.5 数据库格式：从FASTA到可搜索索引
很多人忽略了一个步骤——数据库需要预处理。原始FASTA文件BLAST是无法直接搜索的，必须先用makeblastdb构建索引。这个过程会生成若干后缀的文件，比如对于核酸库：
# 构建核酸库
makeblastdb -in my_sequences.fasta -dbtype nucl -out my_db
# 构建蛋白库
makeblastdb -in my_proteins.fasta -dbtype prot -out my_prot_db生成的文件包括：.nsq/.nhr/.nin（核酸库索引）和 .psq/.phr/.pin（蛋白库索引）。后面章节我们会详细讲怎么把这些索引搬到GPU显存里做加速，但第一步是确保你手中的数据库格式正确。
到目前为止，我们已经把BLAST的核心概念、常用变体、关键参数和基本使用方法过了一遍。这一套是后续所有GPU加速优化的基础，跑偏了后面调参越调越痛苦。嗯，下节课咱们直接进入BLAST在GPU上的运行原理，看看把索引放进显存到底能快多少。