2、BLAST原理与数据库

做基因组注释，第一关往往就是BLAST。你拿到一条新序列，不知道它是什么，第一反应就是“拿去BLAST一下”。我刚开始做生物信息时，也是这么干的——把序列往NCBI网页里一贴，点下按钮，然后等着结果出来。后来发现，如果不理解BLAST背后怎么算的，你很可能被一堆低置信度的hit带偏方向。

这一节咱们就把它拆开看看：BLAST的算法到底怎么跑？NCBI那几个核心数据库（nr、nt、refseq）各有什么脾气？五种BLAST类型到底该用哪个？

2.1 BLAST算法核心思想

BLAST全称是Basic Local Alignment Search Tool。说白了，它不追求全局比对，只找局部的相似片段。为什么？因为生物学上真正保守的功能区域往往是局部的——比如一个蛋白的激酶结构域，长度也就一两百个氨基酸。你拿全长去比，反而把噪声引入进来了。

核心算法分三步：种子生成、无空位扩展、带空位扩展 + 打分。

1. 种子生成（Seeding）：把查询序列切成固定长度的小片段（比如11个核苷酸或3个氨基酸），称为“word”。然后在数据库里快速匹配这些word。我记得早期版本默认word长度是11bp（核酸）和3aa（蛋白），后来可以调参数。
2. 无空位扩展（Ungapped Extension）：找到种子匹配后，向左右两侧延伸，不引入空位，直到打分低于某个阈值。这一步为了快速筛掉大量无意义匹配。
3. 带空位扩展（Gapped Extension）：对通过前两步的候选区域，用动态规划（Smith-Waterman变形）做带空位的精确比对，给出最终分数。

空位罚分在第二步里特别关键。我曾经调过一次BLASTX的-罚分参数，结果出来的比对质量天差地别。原则是：空位开放罚分（Gap Open）要比扩展罚分（Gap Extend）高很多，因为引入一个空位是大事，但一旦有了，后续延续相对便宜。

2.2 NCBI数据库介绍

NCBI维护着几十个序列数据库，做基因组注释最常用的就三个：nr、nt、refseq。它们之间的区别，说白了就是“冗余程度”和“注释可信度”的取舍。

数据库	数据类型	特点	我常用的场景
nr	非冗余蛋白库	整合GenBank、RefSeq、Swiss-Prot、PDB等，去冗余但不绝对	新基因功能注释，尤其是未知蛋白
nt	非冗余核苷酸库	类似nr，但针对核酸序列，数据量极大	鉴定物种来源或找同源基因
refseq	参考序列库	经过人工审核，非冗余，质量高，但覆盖不全面	构建训练集、做标准注释

我建议：日常快速初步注释用nr，因为数据量大，容易找到相似序列。做论文或正式结果用refseq，每条序列都有唯一ID，可溯源。至于nt，我很少直接用它BLAST，除非我要确认一条序列是不是已知的某个物种的转录本。

举个例子：我在做水稻基因组注释时，先拿基因模型去nr库里BLAST，看能不能找到水稻已知蛋白的hit。如果nr里都是拟南芥的，而且一致性很低，那我就知道这个基因模型可能预测错了。

2.3 BLAST类型对比

BLAST有五个主要“变体”：blastn、blastp、blastx、tblastn、tblastx。每个对应不同的输入序列和数据库类型。我刚学的时候老搞混，后来总结了一张表：

类型	查询序列	目标数据库	典型用途
blastn	核酸	核酸	找同源DNA、EST比对
blastp	蛋白	蛋白	蛋白功能注释，找同源蛋白
blastx	核酸（翻译）	蛋白	预测新基因：将转录本翻译成蛋白再比对
tblastn	蛋白	核酸（翻译）	在基因组DNA中发现潜在编码区域
tblastx	核酸（翻译）	核酸（翻译）	高度保守的核酸区域（成本高，少用）

我做环境宏基因组时最常用blastx，因为测出来的是DNA片段，直接比对到蛋白库可以绕过内含子问题。相反，做进化分析时常用tblastn：拿一个已知蛋白，在近缘种的基因组里找可能的直系同源。

有一点要小心：tblastx会把两边的核酸都翻译成6种框架，计算量爆炸。我最多只用它做几Kb的小片段比对，一旦超过10Kb，跑着跑着内存就满了。

# 典型的命令行用法
# blastx：查询核酸，比对到蛋白库
blastx -query transcript.fa -db nr -out blastx_results.txt -outfmt 6 -evalue 1e-5

# tblastn：查询蛋白，比对到核酸库
tblastn -query protein.fa -db refseq_genomic -out tblastn_results.txt -outfmt 6 -evalue 1e-3

嗯，现在你应该对BLAST的“基因”有个底了。下面这张图是我画的BLAST知识体系，帮你把算法、数据库、类型串起来：