从水果连连看到两条序列比对

一、序列比对 Sequence Alignment

序列比对（sequence alignment），目前是生物信息学的基本研究方法。

算法类似于连连看，规则是上下两个水果一样，就可以连起来，计如得分：

现在如果上下两行代表两条序列，把水果换成碱基，可消除的碱基中间连线，就像下面这样：

到目前为止，我们已经实现了一个简单的序列比对。

序列比对最终结果可以用比对得分来评估，然后通过统计学分析后，得到序列间的相似性与同源性，以及它们的显著性水平即可进行下一步生物信息分析。

在应用上，如果找到了不同序列之间的相似性，那就可以推断功能或建立进化关系，以此更好地了解基因的起源和功能。

反之，如果找到序列间的不相似性，就能推断插入，突变，缺失等生物学过程，比如推断新冠病毒的突变位点。

根据序列比对范围和目的，分为两种：

1、全局比对 Global Alignment

顾名思义，就是对两条序列的全长都进行比对

当然有时候序列比对时会不尽人意，类似于这样：

细心的小伙伴可能会发现只有在其中空一格就会“连”到更多的碱基

其中引入的空格，也叫空位（Gap），在生物学中也有依据：DNA 序列在进化过程中会发生的碱基删除事件。但更多的是一种算法规则，即在计算打分时，需要遵循以下规则：

• 碱基匹配加分：+1，也就是中间有连线的碱基对
• 碱基错配罚分：-1，也就是没有连线的碱基对
• 碱基空位罚分：-3，也就是空位，不组成碱基对

根据规则，上述的比对结果为：8-1-3=4

这种比对常常用于基因家族分析，系统发育树构建等

2、局部比对 Local Alignment

目的是在两条序列比对后，获取序列比对分数或置信度最高的匹配序列片段。我们拿刚刚的全局比对结果举例

最佳的片段为

根据规则计算比对分数为：7，显然这就是我们要找的最佳匹配序列。

这种比对常常用于功能域查找，转录因子足迹搜索等等。

！！！

但是有个坏消息是，现实中的序列是要长的多，比如癌基因 p53 的序列长度为 25760 个碱基。

而且也要复杂的多，比如氨基酸之间的错配，由于氨基酸之间的物化性质，虽然是不同的氨基酸，但是介于错配与匹配分数之间，这就不能用简单的规则来计算比对分数啦。

为了获得最佳的比对序列，就需要比较序列间的比对得分大小。那么现在有两个需要解决的问题：

• 设计一种规则，用于计算最真实的比对得分
• 设计一种算法，来快速精准的比对序列

这时，有大牛提出计分矩阵和最优比对算法来解决这两个问题。

这篇我们先来探讨比对的得分的计算，也就是计分矩阵的由来与计算方法：

二、计分矩阵 Scoring Matrix

在序列比对过程中，需要一个计分规则来对匹配到的每个位置的碱基，氨基酸，错配等进行打分，因此该矩阵也叫替换矩阵（substitution matrix）。

2.1 碱基计分矩阵

比如我们来计算下面两条 DNA 序列的分值：

一个常用与DNA序列的计分矩阵

这个矩阵的意思就是，只有匹配一样计分 0.9，但凡不一致扣分 0.1。根据这个计分规则，我们可以很骄傲地拿出我们幼儿园学到的数学，得出：

同源性：代表两条序列间有进化关系，也就是进化中的突变概率也会被考虑在内。

相似性：只代表两条序列的相似度

空位问题 Gap

对于序列在进化过程中，插入或缺失造成的序列空位，可能是一个或多个碱基，氨基酸，甚至功能域。

在进行序列比对时，就需要考虑到这些问题，一般用空位罚分（Gap penalty）来处理。

用公式表示：

这里的参数代表：

• a 空位设置 Gap opening：只要有一个空位出现，就以空位设置罚分
• b 空位扩展 Gap extension：任一空位的扩大，以空位扩展罚分，一般长度越大，罚分越重
• k 空位数

很容易可以得到

2.2 氨基酸计分矩阵

蛋白质序列的计分矩阵相较于只有 4 个碱基的 DNA 序列要复杂的多。对于蛋白质序列的分析可以避免在翻译时的简并性问题（几个三联体可能编码同一个氨基酸），而且也是最接近执行生物学功能和自然选择的分子，因此蛋白质分析在同源家族，基因进化等方面更具研究意义。

2.2.1 PAM 矩阵

各种氨基酸在进化过程中，由于其自身的物化性质，一种氨基酸替换为另一种氨基酸的概率并不一样。

1978年，以 Dayhoff 为首的科学家，对大量蛋白质家族进行统计学分析，观测到 1572 次氨基酸替换，构建了最原始的 PAM 矩阵（Percent Accepted Mutation，PAM），也叫 MDM 矩阵（Mutation Data Matrix）或 Dayhoff 矩阵。

根据该表可以计算突变概率矩阵，其中每个矩阵元素代表在进化过程中氨基酸之间的替换频率。

接下来将观测到的突变百分率作为一种进化时间单位，现在假设同一位点不会发生两次以上的突变，我们对 PAM 进行次方处理，比如 PAM 的 100 次方，意味着进化了 100 次，时间尺度也会更大。

在Dayhoff 和她的小伙伴研究过程中，发现将突变概率矩阵进行 250 次方处理后得到的 PAM 250，适合用于研究远缘蛋白质进化，换句话说这是一个研究这种蛋白质最合适的时间尺度。

然后再将 PAM 250 矩阵进行对数处理，得到 PAM250 的对数概率矩阵，该矩阵用于表示氨基酸间互相替换的观测规律。

经过前人的不懈努力，我们终于拿到了最终的计分矩阵，可以计算比对得分啦。

后来随着蛋白质序列的增加，有人扩大了统计样本，新构建了 JTT 矩阵等，但最终效果都与 PAM 类似。因此，目前使用最为广泛的还是 PAM。

不清楚选择哪种矩阵怎么办？

需要注意，由于不同的蛋白质家族进化速度并不相同，因此选用的 PAM 也会不一样。

总的来说，如果研究进化关系远的蛋白质序列，最好选 > 100 PAM。

如果蛋白质序列本身相似度高，PAM 的影响会比较小，默认 PAM250 就行。

2.2.2 BLOSUM 矩阵

1992 年，Steve Henikoff 和他的妻子 Jorja Henikoff 一起引入了 BLOSUM 替换矩阵。

BLOSUM(Blocks Substitution Matrix)，不同于 PAM 直接用全序列对统计氨基酸替换规律，Henikoff 先对区域保守的蛋白质家族进行局部对比，得到蛋白质序列的高度保守区域，也就是 Blocks，然后基于局部比对块获得每个位置的替换分数。

在计算时首先要构建一个蛋白质家族最保守区域的序列比对数据库，得到局部比对块，计算块中的氨基酸对。

现在计算每个氨基酸对的替换分数：

a. 计算观察概率

假设 代表 i，j 氨基酸对， 代表观察到的氨基酸频率:

b. 计算期望概率

在完全独立情况下， 代表 i，j 氨基酸频率，该氨基酸替换频率的期望值

其中，

c. 计算氨基酸替换分数

也就是， 代表观察到的可能性， 代表预料之中的可能性。每个氨基酸对的出现与该对出现的预期值的比率，再被四舍五入并用于替换矩阵中，得到这样一种矩阵，类似于 PAM 矩阵：

其中，

• 零分表示在数据库中发现给定的两个氨基酸比对的频率是偶然的
• 正分表示比对被发现的频率高于偶然
• 负分表示比对被发现的频率低于偶然

根据构建数据库时， Block 的最小相似比例，可以定义不同的 BLOSUM。

比如上图所示的是 BLOSUM 62， 代表使用局部序列比对相似度为至少为 62%的 Block 构建的替换矩阵。

不难理解，BOLSUM 后的数字越大，代表进化关系越近。

不清楚选择哪种矩阵怎么办？

BLOSUM 80：进化关系近的蛋白质

BLOSUM 62：大部分蛋白质，默认

BLOSUM 45：进化关系远的蛋白质

目前 BLOSUM 矩阵就是我们使用 BLAST 算法经常用的一种计分矩阵。

2.2.3 PSSM 位置特异性矩阵

位置特异性矩阵（PSSM，Position-Specific Scoring Matrix），计算每种碱基或氨基酸，在特定位置的频率矩阵。一般常用于保守序列的搜索，比如PSI-BLAST，DELTA-BLAST等。

计算举例：

a. 计算频数矩阵

对每条序列中的 ATCG 计算频数

b. 计算频率矩阵

c. 标准化矩阵

碱基背景频率/位置频率

d. 对数转换矩阵

，其中 代表标准化矩阵中对应的频率，如果 为 0，则取 伪值 0.1

实际应用中，可以用该 PSSM 矩阵作为打分矩阵，去搜索未知序列中与该 PSSM 相似的保守序列片段。当然也可以计算并绘制基序 Logo。

下一篇我们详细探究序列比对算法：

• Needleman-Wunsch 算法
• Smith-Waterman 算法