《生物信息学 第四版》读书笔记·ch03 序列比对原理
生信中，常用的研究手段之一是比对。例如双序列比对寻找相似区域、保守位点，分析进化关系，多序列比对寻找保守区域、位点、模式。
3.1 序列比对相关概念
序列比对（sequence alignment），指运用某种数学模型或算法，寻找多个序列之间的最大匹配碱基数。比对结果反映了算法在多大程度上提供序列之间的相似性关系及生物学特征，是序列分析、数据库搜索的基础。
序列比对类型。从比对序列的数目上分双序列比对、多序列比对。双序列比对从比对范围分全局比对（global alignment）、局部比对（local alignment）。
序列差异由突变引起，常见的突变有替换、插入、删除。
字符编辑操作（edit operation），通过编辑操作将一个序列转换为另一个新序列，用于解决插入、删除字符等问题。分成以下类型：
Match(a, a)，字符匹配。Delete(a, -)，删除。从序列1删除字符a或从序列2插入-。Replace(a, b)，替换。用序列2的字符b替换序列1的字符a。Insert(-, b)，插入。在序列1插入-或删除序列2的b。
编辑距离（edit distance），指通过编辑操作计算获得两条序列之间的距离。
代价函数或得分函数，用于评价比对结果，规则如下：
- 两条序列比对的得分或代价，等于所有编辑操作的得分或代价总和。
- 最优比对是总得分最高的比对。
- 最小编辑距离是在得分函数值最优时的距离。
双序列比对，指对两条序列进行编辑操作，使两条序列长度相等，编辑距离尽可能小，匹配字符尽可能多。
全局序列比对，指对给定序列全长进行比较，算法与最长公共子序列问题类似。给定序列$A, B$，用$S(i, j)$表示两序列所有比对的最好分数，可以用递归关系计算，见P47。主要优势在于对具有高度同源性的序列进行优化，例如根据已知三维结构的同源序列对未知结构的序列进行建模。
局部序列比对，指获取特定序列在数据库中配对最好的亚区，具体算法同全局比对，见P47。注意，全局比对要从序列前端开始，局部比对从任何一个地方都可以作为起点，适合具有局部同源的片段，例如搜索特定序列、结构域、重复序列等。
相似性（similarity），指两序列间直接的数量关系，例如部分相同、相似的百分比等度量。
同一性（identity），指两序列在同一位点核苷酸碱基或氨基酸残基完全相同的序列比例。
同源性（homology），指从共同祖先经过趋异进化形成的不同序列，即从数据中推断出两个基因在进化上具有共同祖先，属于质的判断。
高相似性不等于同源，但大多数情况下表明了同源性。
同源基因，指来自共同祖先的基因。通常，一个物种基因组中，两个基因或可读框在各自全长的60%以上范围，同一性大于等于30%，称为同源基因。包括：
- 直系同源基因（orthologous gene），指不同物种中具有相同功能的同源基因，在物种形成过程中形成。由共同祖先演化而来，具有序列相似性。
- 旁系同源基因（paralogous gene），指一个物种内的同源基因，是种内基因倍增的结果。
序列相似性高时，直系同源暗示功能性同源，旁系同源有相似但不相同的功能。实践中，单从序列区分直系同源、旁系同源较难。
3.2 序列比对打分方法
对于序列比对结果，通常用打分矩阵计算其分值，评价结果优劣。打分矩阵是序列比对的基础，选择不同打分矩阵得到的结果不同。通常保守替换比随机替换更能维持蛋白功能，因此理化性质相近的AA替换应该比理化性质远的AA替换得分高，保守的AA替换得分高于非保守AA替换。
构建相似性打分矩阵。基于远距离进化过程中观测的残基替换率，用不同分值表征不同残基之间的相似性程度。注意，打分矩阵有其固有的噪声，即对2个不同残基赋予某个分值时引入了比对过程的噪声。
DNA打分矩阵。常用的有单位矩阵、0.99, 0.0033组成的矩阵、0.99, 0.006, 0.002组成的矩阵（考虑转换、颠换）。
单位矩阵（稀疏矩阵，单一相似性打分矩阵）。只考虑残基的同一性，即相同残基为1，不同残基为0，可看作0和1的打分矩阵。
AA打分矩阵。蛋白质进化过程中，不同AA相互替换的概率不同。一种方式是，Dayhoff将20种AA分成6组，同一组的AA看作相同。另一种方法是统计自然界中AA的替换率，例如PAM矩阵、BLOMSUM矩阵。
PAM矩阵（Dayhoff突变数据矩阵）。基于进化原理，建立在进化单点可接受突变（point accepted mutation，PAM）模型基础上。PAM是进化的变异单位，1个PAM表示每100个残基中有1个可接受单点突变。Dayhoff用统计方法得到PAM 1矩阵，外插到PAM 250，见图3-3。即，基于相似性高（85%及以上）的序列比对获得初始模型，再计算进化距离远的矩阵，准确性受一定限制。
BLOSUM矩阵（blocks substitution matrix，模块替换矩阵）。Henikoff夫妇从蛋白质模块数据库BLOCKS中找出的一组替换矩阵。设置最小相同残基数百分比，将序列片段整合在一起，避免同一残基的重复计数。每一个片段中，计算每个残基位置的平均贡献，使整个片段可以看作单一序列。设置不同百分比，产生不同矩阵。例如80%以上相同残基的序列组成的序列模块产生了BLOSUM80矩阵。通过试验确认了BLOSUM 62矩阵的最佳性能，实践中，BLOSUM 62是许多蛋白质比对工具的标准。BLOSUM矩阵的相似性根据真实数据产生。
序列分析的难点是确定只有20%相似性的序列之间是否同源。根据残基差异百分率与进化距离PAM对照表（表3-1，P51），PAM 250突变数据矩阵在20%的水平上反映序列之间的相似性。实际序列比对时，应该选择不同相似性分数矩阵多次比对，比较比对结果。例如，针对不同演化距离，考虑PAM 100 到PAM 500的矩阵，关系近用PAM 100到PAM 150，关系远用更高的。
PAM矩阵和BLOSUM矩阵在高度相似的序列比对上整体差异不大，但在朦胧区（twilight zone，两条序列相似程度在20%左右），可能产生显著影响。此时，增强微弱信号以探测远距离相关变得很重要。
BLOSUM矩阵基于序列局部比对结果建立。考虑到保守区域的比对含有各种不同进化距离的序列，因此用BLOSUM矩阵分析高度保守的残基可能造成结果偏倚。
空位罚分。比对结果中，每插入1个gap就从总分值中减去一定的分值，分为空位起始罚分、空位延伸罚分。序列比对结果的分值是匹配残基的总分值与空位罚分总和。
- 起始空位，指序列比对时，在某条序列中插入1个空位使其有更好的匹配。
- 延伸空位，指引入空位后，继续引入下一个连续的空位使序列之间有更好的匹配。
空位罚分的方式有：
- 线性空位罚分（constant gap penalty），只对起始空位罚分，不考虑延伸空位。
- 仿射空位罚分（affine gap penalty），起始空位对应初始罚分，延伸空位对应较低的罚分，即$g(k) = a + b * k$，其中$a, b, k$分别表示初始罚分、延伸罚分、空位数目。
- 其他罚分系统。
3.3 序列比对算法
序列比对算法，用于从多个比对结果中获取合适的序列比对结果。分成dotplot算法、动态规划算法、blast算法。
dotplot算法。古老的序列比对算法，通过点阵图表示。算法步骤如下：
- 构建点阵矩阵。在矩阵中，将两条序列的碱基/残基分别沿x、y轴排列，依次比较两条序列的每个碱基，相同时则在矩阵中填充点，形成点阵矩阵。在线工具有Dotlet（https://dotlet.vital-it.ch/）。
- 获取相似性片段。将对角线方向上相邻的点连接，直线形成的矩形区域就是这两条序列的相似性片段。
dotplot算法获得的是相同片段，不能提供相似片段在统计学意义上的相似性。
动态规划算法（dynamic programming algorithm），最精确的一种序列比对算法，分成以下两种：
- 全局动态规划算法。Needleman-Wunsch算法，用于发现两条序列全局的相似性。
- 局部动态规划算法。Smith-Waterman算法，用于发现两条序列局部水平上的相似性。
动态规划算法步骤如下：
- 计算得分矩阵。使用迭代方法计算两个序列的相似分值，存储于得分矩阵中。公式见P55。
- 寻找最优比对序列。从最佳路径的终点利用回溯法寻找最优路径，表示两序列的最优比对结果。注意，全局动态规划算法中，最佳路径终点在最后1行、最后 1列的位置。局部动态规划算法中，最佳路径终点在元素值最大的位置。
动态规划算法比对非常精确，但运行时间较长，不适合大量数据的数据库搜索。
例子：见P56。
BLAST算法。Altschul等1990年提出，采用短片段匹配、统计模型寻找目的序列与数据库之间的最佳局部比对效果。基本思想为通过产生数量更少但质量更好的增强点提高速度。
BLAST算法步骤如下：
- 根据查询序列，编译长度固定的字段编译列表。
- 在数据库中搜索与字段匹配的记录。
- 以字段为中心向两端延伸，寻找超过阈值分数S的高分值片段对（high-scoring segment pair，HSP）。
期望值（E值），指在1次数据库搜索中，随机条件下期望发生的得分大于S的不同比对的数目，可用于估计BLAST搜索中的假阳性结果。公式为：$E = Kmne^{-\lambda S}$，其中$m$为待查询序列长度，$n$为整个数据库长度，$S$为比对原始分数，$K$和$\lambda$为Karlin-Altschul统计量。
blast算法是一种近似算法，速度快、比较精确。注意，动态规划算法适合较少序列，blast算法适合从一组大量序列中搜索与查询相似序列。
3.4 序列比对工具
序列比对工具，即序列比对数据库搜索工具，有EBI FASTA、NCBI BLAST。
FASTA工具。FAST-ALL的缩写，算法由Pearon、Lipman在1988年提出，可用于核酸、蛋白质序列快速序列比对数据库搜索工具，包括fasta, fastx/y, fastf, fasts, tfastx/y。搜索速度不如blast，逐步被取代。
BLAST工具。利用blast算法，从核酸、蛋白质序列数据库中找出与目标序列具有一定程度相似性的序列。包括以下工具：
blastn。将目标DNA序列的两条链与DNA数据库比较。blastp。将目标蛋白质序列与蛋白质数据库比较。blastx。将目标DNA序列翻译成6个蛋白质，逐个与蛋白质数据库比较。tblastn。将DNA数据库中的每条序列翻译成6种蛋白质，与目标蛋白质序列比较。tblastx。将目标DNA序列、DNA数据库都翻译成6种蛋白质，进行36次蛋白-蛋白数据库搜索。- PSI-BLAST（position specific iterated BLAST），位点特异性迭代BLAST，用于寻找远缘相关蛋白质序列，比常规
blastp敏感。步骤为：用blastp搜索。对结果进行多序列比对，构建位点特异性矩阵PSSM。根据矩阵PSSM再次搜索目标数据库。位点特异性反复比对后，用缺失比对的参数检验每个匹配的统计显著性。重复2-4步骤，直至不出现新的结果，一般需要重复5次。 - PHI-BLAST（pattern hit initiated BLAST），模式识别BLAST。用于查询与目标序列相似的复合某种模式（pattern）的蛋白质序列。
- MEGABLAST。快速的局部核酸序列比对工具，适合基因预测、单核苷酸多态性等，可有效识别相似性高的序列，例如相似性95%以上的序列，比
blastn快速、准确。
3.5 多序列比对
目的。对3条以上序列进行下比对，发现构成同一基因家族的序列的共性，有助于研究分子结构、功能、进化。例如发现结构域或功能保守片段，发现蛋白质序列的系统发育关系。
定义。在多条序列中插入空位，使全局比对结果具有相同的长度，但结果中不能出现一列全为空位的列。
应用。发现新序列与已知序列家族的同源性，预测蛋白质序列的二级、三级结构，发现蛋白的系统发生关系，分析蛋白质家族结构、功能相似片段等。
多序列比对算法，如下：
- 动态规划算法。类似于双序列比对动态规划算法，只是维数升级为三维。
- 渐进式算法。基于相似序列具有进化相关性假设，步骤为：执行双序列比对，构建距离矩阵，产生指导树（guide tree，GT），执行渐进式比对，对关系密切的序列执行加权。
- 迭代算法。核心是用比对记分函数反复添加一个附加的序列到已知的比对中。首先从双序列比对中找到距离最小的一组，形成最优比对。反复寻找与最优比对距离最小的序列，更新结果。
- 统计概率算法。常用方法为HMM（隐马尔可夫模型，hidden Markov model，HMM），描述大量相互联系的状态之间发生相互转换概率的模型，本质上是一条表示匹配、插入、缺失 状态的链。用于多序列比对时可用于检测序列比对结果中的保守区，需要把所有位置表示为三种状态之一。
多序列比对工具，如下：
- ClustalX/W。使用广泛的工具，采用渐进式算法。
- Clustal Omega。类似于X/W，但提高了可扩展性，支持多线程，比对结果更好。可在几小时完成数十万条序列比对。
- T-Coffee。核酸或AA多序列比对工具，更适合蛋白质序列对比。算法由Jaap Heringa在2000年设计，步骤包括生成基本信息库、扩展库、生成指导树、渐进式比对。比Clustal准确率、敏感性更高，但速度慢。在相似性较高的相关序列比对中，T-Coffee的快速模式速度快。
- MultAlin。（http://multalin.toulouse.inra.fr/multalin/）Florence Cprpet开发，核酸、蛋白多序列比对工具。算法的基本思想为启发式聚类。将序列行双序列比对，进行分层次聚类，再进行多序列比对，根据多序列比对中的双序列比对分值建立指导树。循环直至分值上升，结束。支持多种序列格式，但运行时间长、运算空间要求高，适合序列少、长度短的情况。
- MAFFT（multiple alignment using fast fourier transform），快速的多序列比对工具，算法基于快速傅立叶变换。将序列表示为向量序列，把序列信息看作信号，通过信号处理获得比对结果。特点在于速度快，特别是高度保守序列，准确度方面与T-Coffee相当。支持下载、在线使用。
参考
《生物信息学》第四版 陈铭 chapter03