在我们身体里，存在着肉眼看不见的"生命代码"——DNA、RNA和蛋白质序列。它们由四种碱基或二十种氨基酸排列而成，看似是简单的字符串，却记录着生命的全部奥秘。
生命代码的两大主角
要比对生命序列，首先得搞清楚它们长什么样。DNA序列由A（腺嘌呤）、T（胸腺嘧啶）、C（胞嘧啶）、G（鸟嘌呤）四种碱基组成。RNA将T换成了U（尿嘧啶）。蛋白质序列更为复杂，由20种氨基酸用单字母表示，比如A是丙氨酸，L是亮氨酸。这些字符的排列顺序决定了基因的功能——就像词序决定了句子的意义。
序列比对的基本概念
拿到一段陌生的DNA序列，想知道它的作用？最直接的方法是序列比对（sequence alignment）——把它与已知序列对齐，看看两段序列"像不像"。
相似性背后可能藏着共同的祖先、相似的功能，或进化的痕迹。如果能在已知序列中找到"亲戚"，往往就能推断出它的功能。
两段序列相似，究竟有多相似？科学家用一套打分体系来量化：
- 匹配得分：两个位置的字符相同，得分。
- 错配惩罚：字符不同，扣分。
- 缺口惩罚：需要插入"空白"（gap）对齐时，扣更多分。
相似度（similarity）指的是两个序列在对应位置上，匹配字符的数目占总长度的百分比。例如两条序列经过比对后：
- Sequence1:
C V H K A T - Sequence2:
C I H K - T
有4处一致，1处相似，1处缺口。
- identity = (4/6)*100% = 67%
- similarity = (4+1/6)*100% = 83%
这套评分系统，让科学家能够快速评估比对结果的优劣。
打分矩阵的必要性
如果说序列比对是解读生命密码的艺术，打分矩阵（scoring matrix）就是它的调色盘。比对比的不仅是字符是否"对得上"，还需要考虑生物学意义——哪些替换是自然界常见的？哪些匹配最能反映进化关系？
序列比对中的"相似性"并非一概而论。以蛋白质为例，氨基酸的化学性质、功能、结构对其进化替换的频率有重要影响。比如，亮氨酸（L）替换为异亮氨酸（I）相对常见，因为都是疏水性氨基酸。而亮氨酸被替换为带正电的赖氨酸（K）就很少发生，可能破坏蛋白质功能。
打分矩阵的任务，正是为这些替换赋予合适的分值，使比对结果既符合生物学事实，又有助于后续分析。
常用的打分矩阵
PAM矩阵
PAM（Point Accepted Mutation）矩阵是最早的蛋白质打分矩阵之一。它以"单位进化时间"为基础——1 PAM相当于每100个氨基酸中约有1个发生了可接受的突变。PAM1自乘n次可以得到PAM-n，表示发生了更多次突变。
PAM矩阵的核心思想是：通过研究近亲物种的蛋白质序列，推测不同氨基酸在进化中的替换频率。
- PAM1：适用于非常相似的序列。
- PAM250：适用于亲缘关系较远的序列。
PAM矩阵像一台时间机器，帮我们回溯序列在进化中的替换轨迹。随着PAM数值增加，比对倾向于寻找远亲关系。
BLOSUM矩阵
BLOSUM（BLOcks SUbstitution Matrix）矩阵更关注当前的序列相似性，尤其适用于远亲序列的比对。它通过分析保守序列区块（blocks）中的氨基酸替换频率得出得分。
与PAM不同，BLOSUM不假设进化模型，而是基于实际观测数据生成。BLOSUM编号代表该矩阵由一致度≥该数值的序列计算而来。例如BLOSUM62是由一致度≥62%的序列计算得出。
- BLOSUM62：最常用，适用于大多数情况。
- BLOSUM80：更适合相似性较高的序列。
- BLOSUM45：适合亲缘关系更远的序列。
简单说，BLOSUM后面数字越小，适合比较的序列相似度越低；数字越大，适合比较的序列相似度越高。
选择矩阵的原则
如何选择合适的矩阵？关键在于序列的"亲疏"。
- 相似序列（近亲）：选择细致的矩阵，如PAM1或BLOSUM80，强调微小差异。
- 差异较大序列（远亲）：选择更粗略的矩阵，如PAM250或BLOSUM45，允许更多替换。
这一选择直接影响比对结果的敏感性与特异性。近亲需要精确，远亲需要包容。
如何解读打分矩阵
下表为PAM-250矩阵，对角线上的数值为匹配氨基酸的得分。在其他位置上，≥0的得分代表对应的一对氨基酸为相似氨基酸，<0的是不相似的氨基酸。
| A | B | C | D | E | F | G | H | I | K | L | M | N | P | Q | R | S | T | U | V | W | X | Y | Z | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| A | 2 | |||||||||||||||||||||||
| B | 0 | 3 | ||||||||||||||||||||||
| C | -2 | -4 | 12 | |||||||||||||||||||||
| D | 0 | 3 | -5 | 4 | ||||||||||||||||||||
| E | 0 | 3 | -5 | 3 | 4 | |||||||||||||||||||
| F | -3 | -4 | -4 | -6 | -5 | 9 | ||||||||||||||||||
| G | 1 | 0 | -3 | 1 | 0 | -5 | 5 | |||||||||||||||||
| H | -1 | 1 | -3 | 1 | 1 | -2 | -2 | 6 | ||||||||||||||||
| I | -1 | -2 | -2 | -2 | -2 | 1 | -3 | -2 | 5 | |||||||||||||||
| K | -1 | 1 | -5 | 0 | 0 | -5 | -2 | 0 | -2 | 5 | ||||||||||||||
| L | -2 | -3 | -6 | -4 | -3 | 2 | -4 | -2 | 2 | -3 | 6 | |||||||||||||
| M | -1 | -2 | -5 | -3 | -2 | 0 | -3 | -2 | 2 | 0 | 4 | 6 | ||||||||||||
| N | 0 | 2 | -4 | 2 | 1 | -3 | 0 | 2 | -2 | 1 | -3 | -2 | 2 | |||||||||||
| P | 1 | -1 | -3 | -1 | -1 | -5 | 0 | 0 | -2 | -1 | -3 | -2 | 0 | 6 | ||||||||||
| Q | 0 | 1 | -5 | 2 | 2 | -5 | -1 | 3 | -2 | 1 | -2 | -1 | 1 | 0 | 4 | |||||||||
| R | -2 | -1 | -4 | -1 | -1 | -4 | -3 | 2 | -2 | 3 | -3 | 0 | 0 | 0 | 1 | 6 | ||||||||
| S | 1 | 0 | 0 | 0 | 0 | -3 | 1 | -1 | -1 | 0 | -3 | -2 | 1 | 1 | -1 | 0 | 2 | |||||||
| T | 1 | 0 | -2 | 0 | 0 | -3 | 0 | -1 | 0 | 0 | -2 | -1 | 0 | 0 | -1 | -1 | 1 | 3 | ||||||
| U | 0 | -1 | -3 | -1 | -1 | -2 | -1 | -1 | -1 | -1 | -1 | -1 | 0 | -1 | -1 | -1 | 0 | 0 | -1 | |||||
| V | 0 | -2 | -2 | -2 | -2 | -1 | -1 | -2 | 4 | -2 | 2 | 2 | -2 | -1 | -2 | -2 | -1 | 0 | -1 | 4 | ||||
| W | -6 | -5 | -8 | -7 | -7 | 0 | -7 | -3 | -5 | -3 | -2 | -4 | -4 | -6 | -5 | 2 | -2 | -5 | -4 | -6 | 17 | |||
| X | 0 | -1 | -3 | -1 | -1 | -2 | -1 | -1 | -1 | -1 | -1 | -1 | 0 | -1 | -1 | -1 | 0 | 0 | -1 | -1 | -4 | -1 | ||
| Y | -3 | -3 | 0 | -4 | -4 | 7 | -5 | 0 | -1 | -4 | -1 | -2 | -2 | -5 | -4 | -4 | -3 | -3 | -2 | -2 | 0 | -2 | 10 | |
| Z | 0 | 2 | -5 | 3 | 3 | -5 | 0 | 2 | -2 | 0 | -3 | -2 | 1 | 0 | 3 | 0 | 0 | -1 | -1 | -2 | -6 | -1 | -4 | 3 |
打分矩阵的表格看似复杂，但很直观：
- 正分：替换较为常见，进化上"被接受"。
- 负分：替换较为罕见，可能对功能有害。
分值越高，替换越可能发生或越保守。
核苷酸序列的打分
相比蛋白质，核苷酸序列的打分体系简单许多。通常采用恒定的匹配与错配分值，如等价矩阵：相同核苷酸之间匹配得分为1，不同核苷酸间替换得分为0。
由于等价矩阵不含有碱基的理化信息，且不区别对待不同的替换，在实际序列比较中很少使用，一般只用于理论计算。这种简单模型适用于短序列的快速比对，但对于复杂比对，仍需要引入更多生物学背景，如转换和颠换等。
通过打分矩阵，科学家得以在序列的无尽组合中寻找进化的蛛丝马迹。矩阵的每一个分值，既是自然界千百万年进化的浓缩，又是我们理解生命复杂性的钥匙。
一切生命皆源于演化。基因序列的比对，正是解读生命演化密码的钥匙。比对算法的背后，是科学家对"相似性"的深刻思考与实践创新。接下来介绍从精准到高效的比对方法，看看科学家如何在浩瀚的序列海洋中寻觅隐藏的规律。
全局比对
设想两首长度相近的诗歌——一首是原作，另一首是略有修改的译本。要精确比对它们的异同，我们会从头到尾逐字逐句核对。这便是全局比对：用动态规划方法从起点到终点逐步匹配。
Needleman-Wunsch算法
Needleman-Wunsch算法是全局比对的奠基之作。它通过递归的动态规划，逐步为每种比对可能性打分，最终找到得分最高的路径。这过程如同精心修复残卷，每一步都兼顾对齐、插入和删除的代价，力求整体最优。
虽然精准，但随着比对长度增加，计算成本也急剧上升。全局比对适合短而相似的序列，例如两个同源基因或两个物种间的DNA片段。
局部比对
如果说全局比对是一场全程录制的舞台剧，局部比对则是剪辑精华片段。例如，我们更关心一首诗中某两句绝妙的句子是否在另一首诗中出现过——这便是局部比对的用武之地。
Smith-Waterman算法
局部比对允许我们聚焦于序列中最相似的部分。Smith-Waterman算法基于动态规划，与全局比对的不同之处在于允许得分从零开始，遇到负得分时即停止。这如同一次自由探索，随时放弃无意义的片段，只记录得分最高的"局部之美"。
这种方法对发现功能域、保守片段等局部相似性尤其有效。但计算成本依然较高，因此通常用于精细比对，而非大规模筛选。
启发式算法
科学研究常需在海量数据中快速找到感兴趣的片段。若将全局和局部比对比作精雕细琢的工匠技艺，启发式算法则更像一场大规模的寻宝行动，以速度见长，同时兼顾一定精度。
BLAST
BLAST（Basic Local Alignment Search Tool）是生物信息学领域的明星算法。它将序列比对拆解为两步：
- 种子匹配：将序列切割成小片段（种子），快速扫描数据库，找到能形成无空位完全匹配的区域。
- 延展优化：从这些种子开始，向两端扩展，形成局部比对。
BLAST从头至尾将两条序列扫描一遍找出所有种子，在允许阈值范围内对片段对进行延伸，最终找出高分值片段对（HSPs，high-scoring pairs）。
计算复杂度是n的一次方（n是序列长度）。相比双序列比对的n²复杂度，这大大节省了时间。这种方法极大加速了序列搜索，尤其适合基因组数据库中的比对任务。然而，快速的代价是对低复杂度区域和弱相似性的敏感性略逊。
FASTA
FASTA是BLAST的前辈，采用类似思路。尽管计算速度略慢，但在某些特定场景中敏感性较高。
敏感性与特异性
比对算法的选择，总在"敏感性"和"特异性"间权衡。
- 敏感性：发现真正相似性的能力，偏向"宁多勿漏"。Smith-Waterman和FASTA在这方面表现突出。
- 特异性：排除无意义匹配的能力，偏向"宁缺毋滥"。BLAST在优化参数后可提供较高特异性。
具体应用中，寻找稀有突变时更看重敏感性；进行大规模筛选时倾向于特异性，以减少假阳性干扰。
结语
序列比对不仅是技术，更是一门艺术。科学家从全局到局部，从动态规划到启发式算法，为我们打开了解读生命的窗口。这背后，始终贯穿的是对生物学意义的执着追寻。
无论是全局比对的严谨工匠精神，局部比对的片段发现之美，还是启发式算法的大规模探索，都指向同一个目标——从遗传密码中揭示生命的奥秘。
比对算法的发展告诉我们：科学需要精度，也需要效率。它们不仅让我们找到序列间的相似性，更让我们理解隐藏其中的生物学意义。这是序列比对的魅力所在，也是生命科学不断进步的源泉。