生物信息学入门——“山东大学生物信息学”学习笔记3序列比较

提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档

---

前言

本章节主要介绍什么是序列，序列比对的意义，怎么进行序列比对，序列比对的工具。

一、序列是什么？

序列的英文名字是 sequence，s 就是一条序列：s=abcdefghijklmnopqrstuvwxyz。序列就是字符串，
由字符穿成的串。通常用$s_i$ 代表序列 $s$ 的第 i 个字符，比如$s_4$是 s 的第 4 个字符，也就是 $d$。当然你的字符串是从 0 位开始计数还是从 1 位开始，自己根据需要定！还有，习惯上用 $s^{\prime }$代表 $s$ 的子序列。这些都是大家上计算机课已经熟知的了。
我们研究的序列主要是指蛋白质序列和核酸序列。蛋白质序列是由 20 个不同的字母，也就是 20 种不同的氨基酸排列组合而成。核酸序列是由 4 种不同的字母，也就是 4 种不同的碱基排列组合而成。核酸序列又分为 DNA 序列和 RNA 序列。
生物序列有自己的书写格式，最常用的是 FASTA 格式：第一行是一个大于号“>”开头，后面紧接注释信息；第二行开始就是纯序列部分，除了序列什么都不要写。

二、序列相似性

比较序列相似性的意义：对于一个已知序列但未知结构和功能的蛋白质，如果与它序列相似的某些蛋白质的结构和功能已知，那么就可以推测出这个未知结构和功能的蛋白质的结构和功能。
如何比较：用一致度（identity）和相似度（similarity）这两个指标来定量描述序列有多相似。
一致度：如果两个序列长度相同，那么它们的一致度可以暂时定义为它们对应位置上相同的残基数目占总长度的百分数。
相似度：如果两个序列长度相同，那么它们的相似度可以暂时定义为他们对应位置上相似的残基与相同的残基的数目和占总长度的百分数。

三、替换记分矩阵

替换记分矩阵是反映残基之间相互替换率的矩阵，描述了残基两两相似的量化关系。
替换记分矩阵有很多种，例如， DNA 序列的替换记分矩阵，蛋白质序列的替换记分矩阵。
三种DNA 序列的替换记分矩阵：等价矩阵(unitary matrix)、转换-颠换矩阵(transition-transversion matrix)、BLAST矩阵。
等价矩阵：相同核苷酸之间的匹配得分为 1，不同核苷酸间的替换得分为 0。由于不含有碱基的理化信息和不区别对待不同的替换，在实际的序列比较中很少使用，一般只用于理论计算。
转换-颠换矩阵（转换和颠换的概念：核酸的碱基按照环结构特征被划分为两类，一类是嘌呤，包括腺嘌呤 A 和鸟嘌呤 G，它们都有两个环；另一类是嘧啶，包括胞嘧啶 C 和胸腺嘧啶 T，它们只有一个环。如果 DNA 碱基的替换保持环数不变，则称为转换，比如腺嘌呤 A 替换为 鸟嘌呤 G、或者胞嘧啶 C 替换为胸腺嘧啶 T，也就是嘌呤变嘌呤，嘧啶变嘧啶；如果环数发生变化，则称为颠换，比如腺嘌呤 A 替换为胞嘧啶 C、或者胸腺嘧啶 T 替换为鸟嘌呤 G，也就是嘌呤变嘧啶，或者嘧啶变嘌呤）：在进化过程中，转换发生的频率远比颠换高。转换-颠换矩阵中，转换的得分比颠换要高为-1 分，而颠换的得分为-5 分。
BLAST 矩阵：令被比对的两个核苷酸相同时得分为+5 分，不相同为-4 分，比对效果最好。这个矩阵广泛地被 DNA 序列比较所采用。结论实践经验所得。因为这个矩阵最早应用于BLAST 工具，因此得名 BLAST 矩阵。

五种蛋白质序列的替换记分矩阵：等价矩阵(unitary matrix)、PAM 矩阵(Dayhoff突变数据矩阵)、BLOSUM 矩阵（blocks substitution matrix）、遗传密码矩阵(genetic code matrix, GCM)、疏水矩阵。
等价矩阵：相同得 1 分，不同得 0 分。
PAM 矩阵：PAM 矩阵基于进化原理。如果两种氨基酸替换频繁，说明自然界容易接受这种替换，那么这一对氨基酸替换的得分就应该高。PAM 矩阵是目前蛋白质序列比较中最广泛使用的记分方法之一。基础的 PAM-1矩阵反应的是进化产生的每一百个氨基酸平均发生一个突变的量值，由统计方法得到。PAM-1 自乘 n 次，可以得到 PAM-n ，表示发生了更多次突变。需要根据要比较的序列之间的亲缘关系远近，选择适合的 PAM 矩阵。如果序列亲缘关系远，即序列间会有很多突变，那就选 PAM 后面跟一个大数字的矩阵。如果亲缘关系近，也就是突变比较少，序列间大多数地方都是一样的，那就选 PAM 后面跟一个小数字的矩阵。
BLOSUM 矩阵：BLOSUM 矩阵后面也带有一个编号，有很多种 BLOSUM 矩阵。不同的是，BLOSUM 矩阵都是通过对大量符合特定要求的序列计算而来的。这点和 PAM 矩阵不同的。PAM-1 矩阵是基于相似度大于85%的序列计算产生的，也就是通过关系较近的序列计算出来的。那些进化距离较远的矩阵，如 PAM-250，是通过 PAM-1 自乘得到的。也就是说，BLOSUM 矩阵的相似性是根据真实数据产生的，而 PAM 矩阵是通过矩阵自乘外推而来的。和 PAM 矩阵的另一个不同之处是BLOSUM 矩阵的编号。这些编号，比如 BLOSUM80 中的 80，代表这个矩阵是由一致度≥80%的序列计算而来的。同理BLOSUM62 是指这个矩阵是由一致度≥62%的序列计算而来的。
亲缘关系较近的序列之间的比较，用 PAM 数小的矩阵或
BLOSUM 数大的矩阵；而亲缘关系较远的序列之间的比较，用 PAM 数大的矩阵或 BLOSUM数小的矩阵。
对于关系较远的序列之间的比较，由于 PAM250 是通过矩阵自乘推算而来的，所以其准确度受到一定限制。相比之下BLOSUM 矩阵更具优势。对于关系较近的序列之间的比较，用 PAM 或 BLOSUM 矩阵做出的比对结果，差别不大。如果关于要比较的序列不知道亲缘关系远近，那么就用BLOSUM62 。
遗传密码矩阵：是通过计算一个氨基酸转换成另一个氨基酸所需的密码子变化的数目而得到的。矩阵的值对应为据此付出的代价。如果变化一个碱基就可以使一个氨基酸的密码子转换为另一个氨基酸的密码子，则这两个氨基酸的替换代价为 1；如果需要 2 个碱基的改变，则替换代价为 2；再比如从蛋氨酸（Met）到酪氨酸（Tyr）三个密码子都要变，则代价为 3。遗传密码矩阵常用于进化距离的计算，它的优点是计算结果可以直接用于绘制进化树，但是它在蛋白质序列比对，尤其是相似程度很低的蛋白质序列比对中，很少被使用。
疏水矩阵：是根据氨基酸残基替换前后疏水性的变化而得到的矩阵。若一次氨基酸替换导致疏水特性不发生太大的变化，则这种替换得分高，否则替换得分低。疏水矩阵物理意义明确，有一定的理化性质依据，适用于偏重蛋白质功能方面的序列比对。

四、序列两两比较

比较两个序列的方法：打点法、序列对比法

3.1 打点法

最简单的比较两个序列的方法，理论上可以用纸和笔来完成。
如果要比较下面这两条序列：

把序列 1 整齐的水平书写，把序列 2 整齐的竖直书写，然后依次横横竖竖的比较每一个位置上的残基。相同的话就在这个位置上打个点，不同话，什么也不干。一共需要比较序列 1 的长度乘以序列 2 的长度个位置，也就是进行 90 次比较，最终产生打点矩阵。

放眼全局，红色的线和蓝色的线加起来基本上构成了一条主对角线。由此可以得出结论：序列 1 和序列 2 是比较相似的两条序列。
用这种方法还可以快捷的发现序列中的串联重复序列以及重复的次数。只要数数在半个矩阵中包括主对角线在内的所有等距的平行线的个数，就可以知道重复的次数，而且最短的平行线对应的序列就是重复单元。短的串联复序列具有高度多态性，也就是说不同的个体间重复次数存在差异，而且这种差异在基因遗传过程中一般遵循孟德尔共显性遗传规律，所以快速查找某些特定的短的串联复序列的重复次数可以用于法医学的个体识别或亲子鉴定等领域。

自动打点的软件：Dotlet 软件https://dotlet.vital-it.ch/
默认的颜色方案是在越相似的地方打的点的颜色越浅，越不相似的地方颜色越深。

3.2 序列比对法

3.2.1 序列比对法

意义：比较两条序列可以通过打点法，但是用打点法只能让你大致了解两条序列是否相似，无法定量的描述。如果想要精确地知道两条序列到底有多相似，就需要使用序列比对法。
定义：运用特定的算法找出两个或者多个序列之间产生最大相似度
得分的空格插入和序列排列方案。
举个例子，比如序列 s 和序列 t：

这两条序列的比对就是把 s 和 t 这两个字符串上下排列起来，在某些位置插入空格，这些空格叫空位（gap）。然后依次比较它们在每一个位置上字符的匹配情况，从而使所有位置的得分之和尽可能的高。就是通过插入空位，让上下两行中尽可能多的一致的和相似的字符对在一起。这不是随便摆摆看看就能完成的，需要使用专门的序列比对算法。
根据比对序列的个数可以把序列比对分为双序列比对和多序列比对。双序列比对就是比 2 条，而多序列比对是比 2 条以上。此外根据序列比对的算法不同，双序列比对又分为全局比对和局部比对。全局比对就是全长比较，一个字符都不能落下。而局部比对是比较对得最好的局部，对得不好的部分会被忽略不计。

3.2.2 双序列全局比对及算法

全局比对(global alignment)：用于比较两个长度近似的序列。
经典的全局比对算法是 Needleman-Wunsch 算法。
用 Needleman-Wunsch 算法为序列 p 和序列 q 创建全局比对。
输入值两条序列、替换积分矩阵、空位罚分。
空位罚分就是当字母对空位的时候的得分。原则希望一致或相似的字母尽可能的对在一起，字母对空位的情况和不相似的字母对在一起的情况少出现，所以通常字母对空位会得到一个负分，这个负分就叫做空位罚分。让空位罚分，也就是 gap 分值为-5 分。在比对中没有空位对空位的情况。

创建一个得分矩阵，并根据公式把得分矩阵填满。填满后即得到全局比对。得分矩阵的第一行是序列 p，第一列是序列 q，注意，p 和 q 的前面各留一个空列和一个空行，也就是第 0 列和第 0 行。

给得分矩阵赋值。根据公式： s(0,0)是初始值 0。
第 0 行：$s(0,j)=gap\ast j$
j 从 1 到 m，m 是序列 p 的长度。也就是 $s(0,1)=gap\ast 1=-5$，$s(0,2)=gap\ast 2=-10$，依次类推。第 0 行实际是一种极端情况的假设。也就是当序列 p 全部对空位时的得分。A 对空位是 -5 分，AC 都对空位就累计到了-10 分，ACG 都对空位就累积到了-15 分，如果序列 p 全部对空位，最终的累积得分就是 -25 分。
第 0 列：$s(i,0)=gap\ast i$
第 0 列和第 0 行一样，也是反映了序列 q 如果全部对空位的累计得分。对一个空位累积$gap\ast 1=-5$分，对两个空位累积 $gap\ast 2=-10$ 分，对三个空位累积 $gap\ast 3=-15$ 分，对四个空位累积 $gap\ast 4=-20$ 分。

接下来填 $s(1,1)$。这个格里的值来源于三个值中的最大值。一个是上面格$s(0,1)$里的值加 gap，一个是左面格 $s(1,0)$里的值加 gap，还有一个是斜上格 $s(0,0)$里的值加当前这个位置字母对字母在替换记分矩阵里的分值 $w(i,j)$。就是累积到这个位置时，是字母对字母得分高，还是序列 p 的字母对空位得分高，还是序列 q 的字母对空位得分高？有且只有这三种情况，我们要的是得分最高的那种情况。逐个看一下，上面格 $s(0,1)+gap=-5+-5=-10$。左面格 $s(1,0)+gap=-5+-5=-10$。斜上格 $s(0,0)+w(1,1)=0+10=10$。$max(-10,-10,10)=10$。所以当前这个格 $s(1,1)$的分值就是 10。此外，需要用箭头记录一下这个 10 是从哪里来的。它是从斜上这个格来的，所以我们画一个指向斜上的箭头。


填满之后，右下角的分数就是整个全局比对最终的得分。然后从这个位置开始追溯箭头一直到左上角的零，并且把这些箭头标记出来。

3.2.3 双序列局部比对及算法

局部比对(local alignment)：用于比较一长一短两条序列。
Smith - Waterman算法
局部比对的算法和全局比对很相似，只是在选最大值时通过增加了第四个元素“0”，来达到比对局部的效果。序列 p 和序列 q，一长一短，其他输入值跟全局比对的一样。

局部比对的计算公式在全局比对的基础上增加了第四个元素“0”。得分矩阵初始值仍是0，但第一行和第一列与全局比对不同，全是 0。

从 $s(1,1)$开始要选择四个值中的最大值。除了上面格 $s(0,1)+gap=0+-5=-5$，左边格 $s(1,0)+gap=0+-5=-5$，斜上格 $s(0,0)+w(1,1)=0+-3=-3$，还有一个 0。$max(-5,-5,-3,0)=0$。并且这个 0 既不是从上面格，也不是从左边格，以及斜上格三个方向来的，而是来自于公式里增加的“0”，所以不用画箭头。

按照公式，填充满整个得分矩阵。与全局比对不同，局部比对的得分不是在右下角，而是在整个矩阵中找最大值。这个最大值才是局部比对的最终得分，可能出现在任何一个位置。这次箭头追溯也不是从右下角到左上角，而是从刚刚找到的最大值开始追溯到没有箭头为止。追溯箭头终止的位置也可以是得分矩阵中的任何一个位置。

最后根据标记好的箭头写出比对结果。从左上到右下标记的红色箭头依次是：
斜箭头字母对字母，C 对 C；斜箭头字母对字母，G 对 G。相比这两条序列的全局比对结果，两端的空位在局部比对中就全部被忽略掉了。

四、一致度和相似度

不管两条序列长度是否相同，都要先对它们做全局比对。让两条序列先以最优的方式比对起来，再从全局比对中数出一致字符和相似字符的个数，除以全局比对的长度，来得到它们的一致度和相似度。

五、在线双序列比对工具

EMBL双序列比对工具 http://www.ebi.ac.uk/Tools/psa
全局双序列比对：调整 gap open 和 gap extend 以达到期望的比对结果
局部双序列比对：

六、BLAST搜索

BLAST(Basic Local Alignment SearchTool，基本局部比对搜索工具)是目前最常用的数据库搜索程序。
BLAST基本原理：要点是片段对的概念。所谓片段对是指两个给定序列中的一对子序列，它们的长度相等，且可形成无空位的完全匹配。

BLAST首先找出探测序列和目标序列间所有的匹配程度超过一定阀值的序列片段对，然后对片段对根据给定的相似性阀值进行延伸，得到一定长度的相似性片段，最后给出高分值片段对(high-scoring pairs,HSPs)。改进后的BLAST允许空位的插入。

BLAST种类：BLAST实际上是综合在一起的一组工具的统称，它不仅可用于直接对蛋白质序列数据库和核酸序列数据库进行搜索，而且可以将带搜索的核酸序列翻译成蛋白质序列后再进行搜索，或反之，以提高搜索效率。
Blastp：用蛋白质序列搜索蛋白质序列数据库
Blastn：用核酸序列搜索核酸序列数据库
Blastx：将核酸序列按6条链翻译成蛋白质序列后搜索蛋白质序列数据库
tblastn：用蛋白质序列搜索核酸序列数据库，数据库中的核酸序列要按6条链翻译成蛋白质序列后再搜索。
tblastx：将核酸序列按6条链翻译成蛋白质序列后搜索核酸序列数据库，数据库中的核酸序列要按6条链翻译成的蛋白质序列后再搜索。
根据搜索算法分：标准BLAST，PSI-BLAST，PHI-BLAST等。
标准BLAST：只能找到那些和探索序列十分相近的序列，而其他那些远亲就找不到了。
PSI-BLAST：每次用位置特异权重矩阵(Position-Specific ScoringMatrix，PSSM)搜索数据库后再利用搜索的结果重新构建PSSM，然后用新的PSSM再次搜
索数据库，如此反复(iteration)直至没有新的结果产生为止。
PHI-BLAST(Pattern-Hit InitiatedBLAST,模式识别BLAST):能找到与输入序列相似的并符合某种特定模式（pattern)的序列。

WU-BLAST-WU代表WashingtonUniversity.比NCBI-BLAST更灵敏，在插入空位的算法上更灵活。
Smith andWaterman(SSEARCH):有点儿慢，但是比BLAST更准确。
FASTA:有点儿慢，但是对于DNA序列的比较比BLAST更准确。
BLAT:用于小的序列（如cDNA等）在大基因组中的搜索。

七、多序列比对介绍

多序列比对(multiple alignment)，对两条以上的生物序列进行全局比对。
多序列比对的主要用途：
1.确认：一个未知的序列是否属于某个家族
2.建立：系统发生树，查看物种间或者序列间的关系
3.模式识别：一些特别保守的序列片段往往对应重要的功能区域，通过多序列比对，可以找到这些保守片段
4.已知推未知：把已知有特殊功能的序列片段通过多序列比对做成模型，然后根据该模型推测未知的序列片段是否也具有该功能。
5.其他：预测蛋白质/RNA二级结构等等。
多序列比对的算法：目前所有的多序列比对工具都不是完美的，它们都使用一种近似的算法。为了速度快牺牲了一定准确度。
对要进行多序列比对的序列的几点要求：
1.不支持数量太多的序列。一般10-15条序列，最好别超过50条。
2.不支持关系太远的序列。两两之间序列相似度低于30%的一组序列，作多序列比对会有麻烦。
3.不支持关系太近的序列。两两之间序列相似度大于90%的序列，有再多条都等于只有一条。
4.不支持短序列。多序列比对支持一组差不多长的序列，个别很短的序列属于捣乱分子
5.不支持有重复域的序列。如果序列里包含重复域，大多数多序列比对的程序都会出错，甚至崩溃。
序列的名字建议：
1.名字里不要有“空格”，用“_”代替“空格”。
2. 不要用特殊字符，(比如中文,@,#,&,^等)。
3.名字的长度不要超过15个字符。
4.一组序列里，不要有重名的序列。
5.如果不按上述几点建议命名的话，多序列比对的工具会在不告知你的情况下修改你的序列名称。

八、在线多序列比对工具

Tcoffee：多序列比对工具，算法上与Clustal系列类似，准确度上比Clustal系列
略高，但计算耗时也比Clustal系列略高。
http://www.tcoffee.org
http://tcoffee.crg.cat

多序列比对保存格式
选择保存格式之前，需要考虑：
·大多数软件都支持我要选的这个格式吗?
·我的同事们能使用我选的这个格式吗?
·这个格式能保存我所需要的所有信息吗?
·这个格式适合我进一步加工吗?
如果比对工具输出的格式里没有想要的，可以通过第三方软件格式转换，比如：
fmtseq :http://evol.mcmaster.ca/Pise/5.a/fmtseq.html
http://www.bioinformatics.org/JaMBW/1/2

九、多序列比对的编辑和发布

为了能对多序列比对的结果进行彩色显示和手工编辑，人们开发了多序列比对结果
编辑器。Jalview是一个特别常用的编辑器。
http://www.jalview.org

十、寻找保守区域

从多序列比对中寻找保守区域。
序列标识图(sequencelogo)是以图形的方式依次绘出序列比对中各个位置上出现的残基。每个位置上残基的累积可反应出该位置上残基的一致性。每个残基对应图形字符的大小与残基在该位置上出现的频率成正比。但图形字符的大小并不等于频率百分比，而是经过简单统计计算后转化的结果。
一款流行的创建序列标识图的软件：WebLogo3
http://weblogo.threeplusone.com/


MEME是一款可以自动从一组相关的DNA或蛋白质序列中发现序列基序的软件。http://meme-suite.orgl
在核酸/蛋白质序列中存在有特定模式(pattern)的序列片段，这些片段称为序列的基序(motif)。序列基序与生物功能密切相关。
例如，N-糖基化位点基序(N-glycosylationsite motif)总是符合以下特定模式：Asn开始，然后紧跟除了Pro之外的任何一个氨基酸，再紧跟Ser或者Thr，再紧跟除Pro外的任何一个氨基酸。
这个特定模式可通过正则表达式(Regular expression)来表述:N{P}[ST]{P}
其中N=Asn，P=PrO，S=Ser，T=Thr;(X}代表除X外的任意氨基酸;[XY］代表X或者Y。

PRINTS是蛋白质序列指纹图谱数据库，存储了目前已发现的绝大多数蛋白质家族的指纹图谱。对于一个陌生的蛋白质，只要看看它的序列是否符合某个家族的图谱就可以对它进行分类并预测它的功能。
PRINTS http://www.bioinf.manchester.ac.uk/dbbrowser/PRINTS/
一个蛋白质的指纹(Prints)就是一组保守的序列基序，用于刻画蛋白质家族的特征。这些基序由多序列比对结果获得，且他们在氨基酸序列上是不相邻的，但是在三维结构中，他们可能紧密结合在一起。

总结

以上就是今天要讲的内容，本文仅仅是对序列比对的简单介绍。