第三章 序列比较.pdf

第三章第三章第三章第三章 序列比较序列比较序列比较序列比较 序列比较是生物信息学中最基本 最重要的操作 通过序列比对可以发现生物序列中的功 能 结构和进化的信息 序列比较的根本任务是 通过比较生物分子序列 发现它们的相似性 找出序列之间共同的区域 同时辨别序列之间的差异 在分子生物学中 DNA 或蛋白质的相似 性是多方面的 可能是核酸或氨基酸序列的相似 可能是结构的相似 也可能是功能的相似 一 个普遍的规律是序列决定结构 结构决定功能 研究序列相似性的目的之一是 通过相似的序列 得到相似的结构或相似的功能 这种方法在大多数情况下是成功的 当然 也存在着这样的情况 即两条序列几乎没有相似之处 但分子却折叠成相同的空间形状 并具有相同的功能 这里先不 考虑空间结构或功能的相似性 仅研究序列的相似性 研究序列相似性的另一个目的是通过序列 的相似性 判别序列之间的同源性 推测序列之间的进化关系 这里 将序列看成由基本字符组 成的字符串 无论核酸序列还是蛋白质序列 都是特殊的字符串 本章着重介绍通用的序列比较 方法 3 13 13 13 1序列的相似性序列的相似性序列的相似性序列的相似性 序列的相似性可以是定量的数值 也可以是定性的描述 相似度是一个数值 反映两条序 列的相似程度 关于两条序列之间的关系 有许多名词 如相同 相似 同源 同功 直向同源 共生同源等 在进行序列比较时经常使用 同源 homology 和 相似 similarity 这两个概 念 这是两个经常容易被混淆的不同概念 两条序列同源是指它们具有共同的祖先 在这个意义 上 无所谓同源的程度 两条序列要么同源 要么不同源 而相似则是有程度的差别 如两条序 列的相似程度达到30 或60 一般来说 相似性很高的两条序列往往具有同源关系 但也有例 外 即两条序列的相似性很高 但它们可能并不是同源序列 这两条序列的相似性可能是由随机 因素所产生的 这在进化上称为 趋同 convergence 这样一对序列可称为同功序列 直向同 源 orthologous 序列是来自于不同的种属同源序列 而共生同源 paralogous 序列则是来自 于同一种属的序列 它是由进化过程中的序列复制而产生的 Edited by Foxit Reader Copyright C by Foxit Software Company 2005 2008 For Evaluation Only 序列比较的基本操作是比对 align 两条序列的比对 alignment 是指这两条序列中各个 字符的一种一一对应关系 或字符对比排列 序列的比对是一种关于序列相似性的定性描述 它 反映在什么部位两条序列相似 在什么部位两条序列存在差别 最优比对揭示两条序列的最大相 似程度 指出序列之间的根本差异 3 1 13 1 1 字母表和序列字母表和序列 在生物分子信息处理过程中 将生物分子序列抽象为字符串 其中的字符取自特定的字母 表 字母表是一组符号或字符 字母表中的元素组成序列 一些重要的字母表有 1 4字符 DNA 字母表 A C G T 2 扩展的遗传学字母表或 IUPAC 编码 见表2 3 3 单字母氨基酸编码 见表2 1 4 上述字母表形成的子集 下面所讨论的内容独立于特定的字母表 首先规定一些特定的符号 字母表 A 由字母表 A 中字符所形成的一系列有限长度序列或字符串的集合 a b c 单独的字符 s t u v x A 中的序列 s 序列 s 的长度 Edited by Foxit Reader Copyright C by Foxit Software Company 2005 2008 For Evaluation Only 为了说明序列 s 的子序列和 s 中单个字符 我们在 s 中各字符之间用数字标明分割边界 例 如 设 s ACCACGTA 则 s 可表示为 0A1C2C3A4C5G6T7A8 i s j 指明第 i 位或第 j 位之间的子序列 当然 0 i j s 子序列0 s i 称为前缀 即 prefix s i 而子序列 i s s 称为后缀 suffix s s i 1 有两种特殊的情况 即 i j 或 i j 1 i s i 表示空序列 j 1 s j表示 s 中的第j个字符 简记为 sj 一般认为 子序列与计算机算法中子串的概念相当 但是 严格地讲 子序列与子串的概 念是有区别的 子串是子序列 而子序列不一定是子串 可以通过选取 s 中的某些字符 或删除 s 中的某些字符 而形成 s 的子序列 例如 TTT 是ATATAT的子序列 而 s 的子串则是由 s 中相 继的字符所组成 例如 TAC 是 AGTACA 的子串 但不是 TTGAC 的子串 如果 t 是 s 的子串 则称 s 是 t 的超串 子串也可以称为连续子序列 两条序列 s 和 t 的连接用 s t 来表示 如 ACC CTA ACCCTA 字符串操作除连接操作之外 另有一个 k 操作 即删除一个字符串两端的字符 其定义如 下 prefix s l sk s l suffix s l k s ls i s j ki 1sk s j 序列比较可以分为四种基本情况 具体任务和应用说明如下 1 假设有两条长度相近的 来自同一个字母表的序列 它们之间非常相似 仅 仅是有一些细微的差别 例如字符的插入 字符的删除和字符替换 要求找出这两条序列的差别 这种操作实际应用比较多 例如 有两个实验室同时测定某个基因的 DNA 序列 其结果可能不 一样 需要通过序列比较来比较实验结果 2 假设有两条序列 要求判断是否有一条序列的前缀与另一条序列的后缀相似 如果是 则分别取出前缀和后缀 该操作常用于大规模 DNA 测序中序列片段的组装 3 假设有两条序列 要求判断其中的一条序列是否是另一条序列的子序列 这 种操作常用于搜索特定的序列模式 4 假设有两条序列 要求判断这两条序列中是否有非常相似的子序列 这种操 作可用于分析保守序列 当然 进行序列比较时 往往还需要说明是采取全局比较 还是采取局部比较 全局比较 是比较两条完整的序列 而局部比较是找出最大相似的子序列 3 1 23 1 2 编辑距离 编辑距离 EditEdit DistanceDistance 观察这样两条 DNA 序列 GCATGACGAATCAG 和 TATGACAAACAGC 一眼看上去 这 两条序列并没有什么相似之处 然而如果将第二条序列错移一位 并对比排列起来以后 就可以 发现它们的相似性 如果进一步在第二条序列中加上一条短横线 就会发现原来这两条序列有更多的相似之处 上面是两条序列相似性的一种定性表示方法 为了说明两条序列的相似程度 还需要定量 计算 有两种方法可用于量化两条序列的相似程度 一为相似度 它是两条序列的函数 其值越 大 表示两条序列越相似 与相似度对应的另一个概念是两条序列之间的距离 距离越大 则两 条序列的相似度就越小 在大多数情况下 相似度和距离可以交互使用 并且距离越大 相似度 越小 反之亦然 但一般而言 相似度使用得较多 并且灵活多变 最简单的距离就是海明 Hamming 距离 对于两条长度相等的序列 海明距离等于对应 位置字符不同的个数 例如 图3 1是3组序列海明距离的计算结果 使用距离来计算不够灵活 这是因为序列可能具有不同的长度 两条序列中各位置上的字 符并不一定是真正的对应关系 例如 在 DNA 复制的过程中 可能会发生像删除或插入一个碱 基这样的错误 虽然两条序列的其他部分相同 但由于位置的移动导致海明距离的失真 就图3 1 中例子最右边的情况 海明距离为6 简单地从海明距离来看 两条序列差别很大 整个序列的 长度只有8bp 但是 如果从 s 中删除 G 从 t 中删除 T 则两条序列都成为 ACACACA 这说 明两条序列仅仅相差两个字符 实际上 在许多情况下 直接运用海明距离来衡量两条序列的相 似程度是不合理的 为了解决字符插入和删除问题 引入字符 编辑操作 Edit Operation 的概念 通过编辑操 作将一个序列转化为一个新序列 用一个新的字符 代表空位 或空缺 Space 并定义下述 字符编辑操作 Match a a 字符匹配 Delete a 从第一条序列删除一个字符 或在第二条序列相应的位置插入 空白字符 Replace a b 以第二条序列中的字符 b 替换第一条序列中的字符 a a b Insert b 在第一条序列插入空位字符 或删除第二条序列中的对应字符 b 很显然 在比较两条序列 s 和 t 时 在 s 中的一个删除操作等价于在 t 中对应位置上的一个 插入操作 反之亦然 需要注意的是 两个空位字符不能匹配 因为这样的操作没有意义 引入 上述编辑操作后 重新计算两条序列的距离 就成为编辑距离 以上的操作仅仅是关于序列的常用操作 在实际应用中还可以引入复杂的序列操作 下面 是两条序列的一种比对 上述比对不能反映两条序列的本质关系 但是 如果将第二条序列头尾倒置 可以发现两 条序列惊人的相似 再比如 下面两条序列有什么关系 如果将其中一条序列中的碱基替换为其互补碱基 就 会发现其中的关系 CTAGTCGAGGCAATCT GAACAGCTTCGTTAGT 3 1 33 1 3 通过点矩阵分析两条序列的相似之处通过点矩阵分析两条序列的相似之处 进行序列比较的一个简单的方法是 矩阵作图法 或 对角线作图 这种方法是由 Gibb 首先 提出的 将两条待比较的序列分别放在矩阵的两个轴上 一条在 X 轴上 从左到右 一条在 Y 轴上 从下往上 如图3 2所示 当对应的行与列的序列字符匹配时 则在矩阵对应的位置作出 点 标记 逐个比较所有的字符对 最终形成点矩阵 图3 2序列比较矩阵标记图 显然 如果两条序列完全相同 则在点矩阵主对角线的位置都有标记 如果两条序列存在 相同的子串 则对于每一个相同的子串对 有一条与对角线平行的由标记点所组成的斜线 如图 3 3中的斜线代表相同的子串 ATCC 而对于两条互为反向的序列 则在反对角线方向上有标记 点组成的斜线 如图3 4所示 图3 3 相同子串矩阵标记图 图3 4 反向序列矩阵标记图 对于矩阵标记图中非重叠的与对角线平行斜线 可以组合起来 形成两条序列的一种比对 在两条子序列的中间可以插入符号 表示插入空位字符 在这种对比之下分析两条序列的 相似性 如图3 5所示 找两条序列的最佳比对 对应位置等同字符最多 实际上就是在矩阵标 记图中找非重叠平行斜线最长的组合 图3 5 多个相同连续子序列矩阵标记图 除非已经知道待比较的序列非常相似 一般先用点矩阵方法比较 因为这种方法可以通过 观察矩阵的对角线迅速发现可能的序列比对 实例实例1 1 1 1 实例实例2 2 2 2 两条序列中有很多匹配的字符对 因而在点矩阵中会形成很多点标记 当对比较长的序列 进行比较时 这样的点阵图很快会变得非常复杂和模糊 使用滑动窗口代替一次一个位点的比较 是解决这个问题的有效方法 假设窗口大小为10 相似度阈值为8 首先 将 X 轴序列的第1 10 个字符与 Y 轴序列的第1 10个字符进行比较 如果在第一次比较中 这10个字符中有8个或者8 个以上相同 那么就在点阵空间 1 1 的位置画上点标记 然后窗口沿 X 轴向右移动一个字符 的位置 比较 X 轴序列的第2 11个字符与 Y 轴序列的第1 10个字符 不断重复这个过程 直到 X 轴上所有长度为10的子串都与 Y 轴第1 10个字符组成的子串比较过为止 然后 将 Y 轴的窗口 向上移动一个字符的位置 重复以上过程 直到两条序列中所有长度为10的子串都被两两比较过 为止 基于滑动窗口的点矩阵方法可以明显地降低点阵图的噪声 并且可以明确地指出两条序列 间具有显著相似性的区域 3 1 43 1 4 序列的两两比对序列的两两比对 序列的两两比对 Pairwise SequenceAlignment 就是对两条序列进行编辑操作 通过字符匹 配和替换 或者插入和删除字符 使得两条序列达到一样的长度 并使两条序列中相同的字符尽 可能地一一对应 设两条序列分别是 s 和 t 在 s 或 t 中插入空位符号 使 s 和 t 达到一样的长度 图3 6是对序列 AGCACACA 和 ACACACTA 的两种比对结果以及对应的字符编辑操作 下面就不同类型的编辑操作定义函数 w 它表示 代价 cost 或 权重 weight 对字母 表 A 中的任意字符 a b 定义 这是一种简单的代价定义 在实际应用中还需使用更复杂的代价模型 一方面 可以改变各 编辑操作的代价值 例如 在蛋白质序列比较时 用理化性质相近的氨基酸进行替换的代价应该 比完全不同的氨基酸替换代价小 另一方面 也可以使用得分 score 函数来评价编辑操作 下面给出一种基本的得分函数 在进行序列比对时 可根据实际情况选用代价函数或得分函数 即选用 3 1 式或 3 2 式 下面给出在进行序列比对时常用的概念 1 两条序列 s 和 t 的比对的得分 或代价 等于将 s 转化为 t 所用的所有编辑操作的得分 或代价 总和 2 s 和 t 的最优比对是所有可能的比对中得分最高 或代价最小 的一个比对 3 s 和 t 的真实距离应该是在得分函数 p 值 或代价函数 w 值 最优时的距离 使用前面代价函数 w 的定义 可以得到下列比对的代价 s AGCACAC A t A CACACTA cost s t 2 而使用得分函数 p 的定义 可以得到下列比对的得分 s AGCACAC A t A CACACTA score s t 5 进行序列比对的目的是寻找一个得分最高 或代价最小 的比对 3 1 53 1 5 用于序列相似性的打分矩阵 用于序列相似性的打分矩阵 scoringscoring matrixmatrix 无论是3 1式还是3 2式 都是简单相似性评价模型 在计算比对的代价或得分时 对字符替 换操作只进行统一的处理 没有考虑 同类字符 替换与 非同类字符 替换的差别 实际上 不同类型的字符替换 其代价或得分是不一样的 特别是对于蛋白质序列 某些氨基酸可以很容 易地相互取代而不用改变它们的理化性质 例如 考虑这样两条蛋白质序列 其中一条在某一位 置上是丙氨酸 如果该位点被替换成另一个较小且疏水的氨基酸 比如缬氨酸 那么对蛋白质功 能的影响可能较小 如果被替换成较大且带电的残基 比如赖氨酸 那么对蛋白功能的影响可能 就要比前者大 直观地讲 比较保守的替换比起较随机替换更可能维持蛋白质的功能 且更不容 易被淘汰 因此 在为比对打分时 我们可能更倾向对丙氨酸与缬氨酸的比对位点多些奖励 而 对于丙氨酸与那些大而带电氨基酸 比如赖氨酸 的比对位点则相反 理化性质相近的氨基酸残 基之间替换的代价显然应该比理化性质相差甚远的氨基酸残基替换得分高 或者代价小 同样 保守的氨基酸替换得分应该高于非保守的氨基酸替换 这样的打分方法在比对非常相近的序列以 及差异极大的序列时 会得出不同的分值 这就是提出打分矩阵 或者称为取代矩阵 的原由 在打分矩阵中 详细地列出各种字符替换的得分 从而使得计算序列之间的相似度更为合理 在 比较蛋白质时 我们可以用打分矩阵来增强序列比对的敏感性 打分矩阵是序列比较的基础 选 择不同的打分矩阵将得到不同的比较结果 而了解打分矩阵的理论依据将有助于在实际应用中选 择合适的打分矩阵 以下介绍一些常用的打分矩阵或代价矩阵 3 1 5 13 1 5 1 核酸打分矩阵核酸打分矩阵 设核酸序列所用的字母表为 A C G T 1 等价矩阵 等价矩阵 见表3 1 是最简单的一种打分矩阵 其中 相同核苷酸匹配的得分为 1 而 不同核苷酸的替换得分为 0 没有得分 2 BLAST 矩阵 BLAST 是目前最流行的核酸序列比较程序 表3 2是其打分矩阵 这也是一个非常简单的矩 阵 如果被比的两个核苷酸相同 则得分为 5 反之得分为 4 3 转换 颠换矩阵 核酸的碱基按照环结构分为两类 一类是嘌呤 腺嘌呤 A 鸟嘌呤 G 它们有两个环 另 一类是嘧啶 胞嘧啶 C 胸腺嘧啶 T 它们的碱基只有一个环 如果 DNA 碱基的变化 碱基替 换 保持环数不变 则称为转换 transition 如 A G C T 如果环数发生变化 则称为颠 换 transversion 如 A C A T 等 在进化过程中 转换发生的频率远比颠换高 而表3 3 所示的矩阵正好反映了这种情况 其中转换的得分为 1 而颠换的得分为 5 3 1 5 23 1 5 2 蛋白质打分矩阵蛋白质打分矩阵 设蛋白质的字母表为表2 1 见第二章 1 等价矩阵 3 3 其中 Rij代表打分矩阵元素 i j 分别代表字母表第 i 个和第 j 个字符 2 遗传密码矩阵 GCM GCM 矩阵通过计算一个氨基酸残基转变到另一个氨基酸残基所需的密码子变化数目而得 到 矩阵元素的值对应于代价 如果变化一个碱基 就可以使一个氨基酸的密码子改变为另一个 氨基酸的密码子 则这两个氨基酸的替换代价为1 如果需要2个碱基的改变 则替换代价为2 以此类推 见表3 4 注意 Met 到 Tyr 的转变是仅有的密码子三个位置都发生变化的转换 在 表3 4中 Glx 代表 Gly Gln 或 Glu 而 Asx 则代表 Asn 或 Asp X 代表任意氨基酸 GCM 常用 于进化距离的计算 其优点是计算结果可以直接用于绘制进化树 但是它在蛋白质序列比对尤其 是相似程度很低的序列比对中很少被使用 表3 4遗传密码矩阵 ASGLKVTPEDNIQRFYCHMWZBX Ala A01122111112222222222222 Ser S10112211221121111221222 Gly G11022122112221221221222 Leu L21202121222111122111222 Lys K22220212121111222212122 Val V12112022112122122212222 Thr T11221201221121222212222 Pro P11212210222211222122222 Glu E12121122012212222222122 Asp D12122122101222212122212 Asn N21221212210122212122212 Ile I21211112221021122212222 Gln Q22211221122201222122122 Arg R21111211222110221111222 Phe F21212122222122011222222 Tyr Y21222222211222101132212 Cys C21122222222221110221222 His H22212221211211212022212 Met M22211112222121232202222 Trp W21112222222221221220222 Glx Z22221222122212222222122 Asx B22222222211222212122212 X 22222222222222222222222 3 疏水矩阵 该矩阵 见表3 5 是根据氨基酸残基替换前后疏水性的变化而得到得分矩阵 若一次氨基 酸替换疏水特性不发生太大的变化 则这种替换得分高 否则替换得分低 表3 5蛋白质疏水矩阵 RKDEBZSNQGXTHACMPVLIYFW Arg R10 10 998866655555433333210 Lys K10 10 998866655555433333210 Asp D9910 10 8876665555544433321 Glu E9910 10 8876665555544433321 Asx B888810 10 88887777666555443 Glx Z888810 10 88887777666555443 Ser S66778810 10 10 10 9999887777664 Asn N66668810 10 10 10 9999888777664 Gln Q66668810 10 10 10 9999888777664 Gly G55668810 10 10 10 9999888877665 X555577999910 10 10 10 998888775 Thr T555577999910 10 10 10 998888775 His H555577999910 10 10 10 999888775 Ala A555577999910 10 10 10 999888775 Cys C4455668888999910 10 9999885 Met M3344668888999910 10 10 10 99887 Pro P33446678888899910 10 10 99987 Val V33445577788888910 10 10 10 10 987 Leu L3333557777888899910 10 10 998 Ile I3333557777888899910 10 10 998 Tyr Y2233446666777788999910 10 8 Phe F1122446666777788889910 10 9 Trp W001133444555556777888910 4 PAM 矩阵 为了得到打分矩阵 更常用的方法是统计自然界中各种氨基酸残基的相互替换率 如果两 种特定的氨基酸之间替换发生得比较频繁 那么这一对氨基酸在打分矩阵中的互换得分就比较 高 PAM 矩阵就是这样一种打分矩阵 PAM 矩阵是第一个广泛使用的最优矩阵 它是基于进化 原理的 建立在进化的点接受突变模型 PAM PointAccepted Mutation 基础上 通过统计相似 序列比对中的各种氨基酸替换发生率而得到该矩阵 Dayhoff 和她的同事们研究了71个相关蛋白 质家族的1572个突变 发现蛋白质家族中氨基酸的替换并不是随机的 由此 断言一些氨基酸的 替换比其他替换更容易发生 其主要原因是这些替换不会对蛋白质的结构和功能产生太大的影 响 如果氨基酸的替换是随机的 那么 每一种可能的取代频率仅仅取决于不同氨基酸出现的背 景频率 然而 在相关蛋白中 存在取代频率大大地倾向于那些不影响蛋白质功能的取代 换句 话说 这些点突变已经被进化所接受 这意味着 在进化历程上 相关的蛋白质在某些位置上可 以出现不同的氨基酸 一个 PAM 就是一个进化的变异单位 即1 的氨基酸改变 但是 这并不意味着经过100次 PAM 后 每个氨基酸都发生变化 因为其中一些位置可能会经过多次改变 甚至可能变回到原 先的氨基酸 因此 另外一些氨基酸可能不发生改变 PAM 有一系列的替换矩阵 每个矩阵用 于比较具有特定进化距离的两条序列 例如 PAM 120矩阵用于比较相距120个 PAM 单位的序 列 一个 PAM N 矩阵元素 i j 的值反映两条相距 N 个 PAM 单位的序列中第 i 种氨基酸替换 第 j 种氨基酸的概率 从理论上讲 PAM 0是一个单位矩阵 主对角线上的元素值为1 其它矩 阵元素的值为0 其他 PAM N 矩阵可以通过统计计算而得到 首先针对那些确信是相距一个PAM 单位的序列进行统计分析 得到 PAM 1矩阵 PAM 1矩阵对角线上的元素值接近于1 而其它矩 阵元素值接近于0 例如 可以按下述方法构建 PAM 1矩阵 首先 构建一个序列间相似度很高 通常大于85 的比对 接着 计算每个氨基酸 j 的相对突变率 mj 相对突变率就是某种氨基 酸被其它任意氨基酸替换的次数 比如 丙氨酸的相对突变率是通过计算丙氨酸与非丙氨酸残基 比对的次数来得到 然后 针对每个氨基酸对 i 和j 计算氨基酸 j 被氨基酸 i 替换的次数 最后 将以上替换次数除以对应的相对替换率 利用每个氨基酸出现的频度对其进行标准化 并将以上 计算结果取常用对数 于是得到了 PAM 1矩阵中的元素 PAM 1 i j 这种矩阵被称作对数几 率矩阵 log odds matrix 因为其中的元素是根据每个氨基酸替换率的对数值来得到的 将 PAM 1自乘 N 次 可以得到矩阵 PAM N 虽然 Dayhoff 等人只发表了PAM 250 但潜 在的突变数据可以外推至其他 PAM 值 产生一组矩阵 可以根据待比较序列的长度以及序列间 的先验相似程度来选用特定的 PAM 矩阵 以发现最适合的序列比对 一般 在比较差异极大的 序列时 通常在较高的 PAM 值处得到最佳结果 比如在 PAM 200到 PAM 250之间 而较低值 的 PAM 矩阵一般用于高度相似的序列 实践中用得最多的且比较折衷的矩阵是 PAM 250 5 BLOSUM 矩阵 BLOSUM矩阵是由 Henikoff 首先提出的另一种氨基酸替换矩阵 它也是通过统计 相似蛋白质序列的替换率而得到的 PAM 矩阵是从蛋白质序列的全局比对结果推导出来的 而 BLOSUM 矩阵则是从蛋白质序列块 短序列 比对而推导出来的 但在评估氨基酸替换 频率时 应用了不同的策略 基本数据来源于 BLOCKS 数据库 其中包括了局部多重比对 包含较远的相关序列 与在 PAM 中使用较近的相关序列相反 虽然在这种情况下没有 用进化模型 但它的优点在于可以通过直接观察而不是通过外推获得数据 同 PAM 模型一 样 也有一系列的 BLOSUM 矩阵 可以根据亲缘关系的不同来选择不同的 BLOSUM 矩阵 进行序列比较 然而 BLOSUM 矩阵阶数的意义与 PAM 矩阵正好相反 低阶 PAM 矩阵适 合用来比较亲缘较近的序列 而低阶 BLOSUM 矩阵更多是用来比较亲缘较远的序列 一般 来说 BLOSUM 62矩阵适于用来比较大约具有62 相似度的序列 而 BLOSUM 80矩阵更 适合于相似度为80 左右的序列 第三章第三章第三章第三章 序列比较序列比较序列比较序列比较 3 23 23 23 2两两比对算法两两比对算法两两比对算法两两比对算法 进行序列的两两比对最直接的方法就是生成两条序列所有可能的比对 分别计算得分 或 代价 函数 然后挑选一个得分最高 或代价最小 的比对作为最终结果 但是 两条序列可能 的比对数非常多 是序列长度的指数函数 随着序列长度的增长 计算量呈指数增长 从算法时 间复杂性的角度来看 这种比对方法显然不合适 用上一节中所介绍的点矩阵分析方法 在寻找 斜线及斜线组合时 仍然需要较大的运算量 因此 必须设计高效的算法以找出最优的比对 著 名的Needleman Wunsch 算法 就是针对寻求最佳序列比对这一问题所设计的动态规划寻优策略 下面 首先介绍基本的序列比较算法 即动态规划算法 Dynamic Programming 该算法把一个 问题分解成计算量合理的子问题 并使用这些子问题的结果来计算最终答案 这个算法是生物信 息学的基本算法之一 动态规划是一种常用的规划方法 往往用于在一个复杂的空间中寻找一条最优路径 对于 一个具体的问题 如果该问题可以被抽象为一个对应的图论问题 并且问题的解对应于图中从起 点到终点的最短距离 那么就可以通过动态规划算法解决这个问题 在运用动态规划时 有以下 几个要求 1 首先 搜索问题能够划分成一系列相继的阶段 2 起始阶段包含基本子问题的 解 3 在后续阶段中 能够按递归方式逐步计算前面阶段的每个局部解 4 最后阶段包含全 局解 3 2 13 2 1 序列两两比对基本算法序列两两比对基本算法 设序列 s t 的长度分别为 m 和 n 考虑两个前缀 0 s i和0 t j i j 1 假如已知序列0 s i和0 t j 所有较短的连续子序列的最优比对 即已知 1 0 s i 1 和0 t j 1 的最优比对 2 0 s i 1 和0 t j的最优比对 3 0 s i和0 t j 1 的最优比对 则0 s i和0 t j的最优比对一定是上述三种情况之一的扩展 即 1 替换 si tj 或匹配 si tj 这取决于 si是否等于 tj 2 删除 si 3 插入 tj 令S 0 s i 0 t j 为序列0 s i和与序列0 t j比对的得分 可根据下列递归算式计算最大值 3 4 其初值为 3 5 按照这种方法 对于给定的打分函数p si tj 两条序列所有前缀的比对得分值定义了一个 m 1 n 1 的得分矩阵 D di j 其中 di j S 0 s i 0 t j 对于一个长度为 n 的序列 有 n 1个前缀 包括一个空序列 所以 得 分矩阵的大小为 m 1 n 1 其中 矩阵的纵轴方向自上而下对应于第一条序列 s 横轴方 向从左到右对应于第二条序列 t 矩阵横向移动表示在纵轴序列中加入一个空位 纵向的移动 表示在横轴序列中加入一个空位 而斜对角向的移动表示两序列各自相应的字符进行比对 注意 各轴第一个元素的索引下标为0 di j的计算公式如下 3 6 di j最大值的三种选择决定了各矩阵元素之间的关系 用图3 7表示 矩阵右下角元素即为期望的结果 dm n S 0 s m 0 t n S s t 首先初始化得分矩阵 D 然后计算 D 的其它元素 计算过程从d0 0开始 可以是按行计算 每行从左到右 也可以是按列计算 每列从上到下 当然 任何计算过程 只要满足在计算di j 时di 1 j di 1 j 1和di j 1都已经被计算这个条件即可 在计算di j后 需要保存di j是从di 1 j d i 1 j 1或di j 1中的哪一个推进的 或保存计算的路径 以便于后续处理 上述计算过程到dm n结 束 与计算过程相反 求最优路径或最优比对时 从dm n开始 反向前推 假设在反推时到达 di j 现在要根据保存的计算路径判断di j究竟是根据di 1 j di 1 j 1和di j 1中的哪一个计算而得 到的 找到这个点以后 再从此点出发 一直到d0 0为止 走过的这条路径就是最优路径 即得 分最大路径 其对应于两条序列的最优比对 假设我们采用公式 3 2 的打分函数 即序列字符匹配得分为 1 失配得分为0 空位罚分 为 1 并且假设待比较的两条核苷酸序列分别是 s AGCACACA 和 t ACACACTA 从公式 3 4 或公式 3 6 来看 动态规划算法的执行过程实际上是逐步计算得分矩阵 D 每一个元素 的过程 首先按照公式 3 5 对矩阵 D 进行初始化 即首先计算矩阵第0行和第0列的元素值 其中 d0 0代表序列 s 和 t 两个空前缀比对的得分 其值显然为0 第0行的其它元素d0 j表示序 列 s 空前缀与序列 t 前面连续 j 个字符组成的前缀比对的得分 相当于在序列 s 前面插入了 j 个 空位 而每个空位的罚分为 1 所以d0 j j 即矩阵第0行元素值依次减1 同理 矩阵第0列元 素值也应该依次减1 矩阵初始化的结果见图3 8 t t s s A AC CA AC CA AC CT TA A 0 1 2 3 4 5 6 7 8 A A 1 G G 2 C C 3 A A 4 C C 5 A A 6 C C 7 A A 8 图3 8 序列 AGCACACA 和 ACACACTA 比对的初始化得分矩阵 初始化以后 假设逐行计算得分矩阵 D 的其它元素值 假设现在计算矩阵第1行 第1列的 元素值 即计算矩阵 1 1 的元素值 可以通过三种途径到达该位置 1 从 0 0 出发 经过两条序列第一个字符的比对 A A 2 从 0 1 出发 在第二条序列加上的空位 A 3 从 1 0 出发 在第一条序列加上的空位 A 因此 矩阵 1 1 的元素值来自于下列三个值中最大的一个 1 左上方 0 0 位置的值加上匹配 A A 的得分1 和为1 2 上方 0 1 位置的值加上空位罚分 1 和为 2 3 左边 1 0 位置的值加上空位罚分 1 和为 2 最终 矩阵 1 1 的元素值等于1 当所有元素值计算完以后 形成图3 9所示的得分矩阵 计算完得分矩阵 D 后 从元素 8 8 所在的位置反推最优路径 在图3 9中画出了一条穿 越得分矩阵的路径 该路径表明如何通过合理的比对达到最大的得分 其中 斜线表示匹配或替 换 垂直线表示删除 而水平线则表示插入 由该路径可以得到下面这种序列比对 从图3 9可知 总的比对得分为5 值得注意的一点是 在有些情况下 最优比对并非唯一 亦 即存在几条最优路径 以上计算是在打分函数的基础上进行的 得分值表示序列之间的相似程度 在实际应用中 也可以利用 代价 函数进行计算 两条序列比对的代价越低 序列就越相似 比对的代价越高 序列的差异就越大 因为计算方式刚好与得分函数相反 所以 具体计算时应求出最小代价所对 应的路径 一般来讲 由于比对的得分可正可负 使用起来就更加灵活 所以大量的序列比对实 用程序在计算时都采用得分的概念 在计算序列相似程度时还应该考虑序列长度的影响 令 S s t 表示两条长度分别为 m 和 n 的序列的相似性得分 假设 S s t 99 两条序列有99个字符一致 如果 m n 100 则可以说 这两条序列非常相似 几乎一样 但是 如果 m n 1000 则仅有10 的字符相同 所以 在实 际序列比较时 使用相对长度的得分就更有意义 定义 3 7 以 sim s t 作为衡量序列相似性的指标 从所使用的数据结构 di j本身及其计算过程来看 序列两两比对基本算法的空间复杂度和时 间复杂度都是 O mn 如果两条序列的长度相近 空间和时间复杂度就为 O n2 3 2 23 2 2 子序列与完整序列的比对子序列与完整序列的比对 序列比较的基本动态规划算法找出两条序列的最佳比对 而不在乎它是否具有生物学意义 有些同源序列虽然全序列的相似性很小 但是存在高度相似的局部区域 如果在进行序列比对时 注重序列的局部相似性 则可能会发现重要的比对 因此 不能仅仅只关注全局最佳的那一个 Smith 和 Waterman 在 Needleman Wunsch 算法的基础上进行改进 提出序列局部比对算法 后来 其他人又进一步改进 形成改良 Smith Waterman 算法 该算法将寻找多种最好的但不相互交叉 的比对方式作为结果 下面分别介绍各种局部比对方法 在有些情况下 我们需要将一个较短的序列 或探测序列 或模式序列 与一个较长的完 整序列比较 试图找出局部的最优匹配 假设我们希望在较长序列 ATGCAGCTGCTT 中搜寻短 序列 AGCT 在所有可能的序列比对中 我们感兴趣的是 这之所以是我们最感兴趣的比对 是因为它表明了 较短的序列完整地出现在较长的序列 之中 我们有时希望避免对序列一端或两端出现的空位进行罚分 例如 在寻找一条短序列和整 个基因组的最佳比对时就希望这样 这是一种局部的比对 local alignment 其定义如下 给定两条序列0 s m和0 t n 从 t 中寻 找一个子列i t j使得S s i t j 最大 0 i j n 已有许多高效的算法可解决局部比对问题 而动 态规划算法也只要作一点小小的修改就可以用于局部比对 假设在下面讨论局部序列比对时依然使用公式 3 2 的打分函数 局部比对意味着不计删 除序列 t 前缀 或前缀0 t i与空位比对 的得分 这表明在对动态规划算法之得分矩阵进行初始 化时 按下述方式处理 3 8 局部比对也不计删除序列 t 后缀 或后缀j t n与空位比对 的得分 即 3 9 由式 3 8 可知 在得分矩阵初始化时 对第0行进行如下处理 0 j n 3 10 而其他行 除最后一行 的计算不变 由公式 3 9 可知 最后一行的计算应该是 3 11 同样 dm n依然是最优局部比对的得分 而匹配的子列i t j按如下方式寻找 j min k dm k dm n 然后由位置 m j 出发 反推比对路径 最终通过斜线 非空位 到达 0 i 3 2 33 2 3 寻找最大的相似子序列寻找最大的相似子序列 上面讨论的是两条序列比对的第一种变化形式 即从两条完整序列的比较演化为一条序列与 另一条序列的一部分进行比对 第二种变化形式是对两条序列都进行部分比较 例如 假设 s 和 t 是两条蛋白质序列 并且已知 s 和 t 具有功能上相关的子序列 而 s 和 t 的其他部分与该功 能无关 又如 假设一条很长的黑猩猩 DNA 序列 要求找出其中与人类基因组具有相似部分的 任何一条子序列 对于这种情况 采用全局序列比对方法 不可能找出高度相似的局部区域 需 要设计序列局部相似性的比较算法 下面假设得分函数只奖励匹配 即匹配奖励分值为 1 失 配罚分为 1 空位罚分为 1 这里使用的数据结构依然是一个 m 1 n 1 的矩阵 D 但是 对数组元素含义解释 与基本算法的有所不同 每个元素的值代表序列0 s i某个后缀和序列0 t j某个后缀的最佳比对 同子序列与完整序列比对一样 这种局部比对不计前缀的得分 所以新的边界条件是 另外 由于0 s i和0 t j总有一个得分为 0 的空后缀比对 见图3 10 因此 矩阵 D 中的 所有元素大于或等于 0 于是 新的递归计算公式为 3 14 阈值 0 意味着矩阵中的 0 元素分布区域对应于不相似的子序列 而正数区域则是局 部相似的区域 最后 在矩阵中找最大值 该值就是最优的局部比对得分 它所对应的点为序列 局部比对的末点 然后 反向推演前面的最优路径 直到局部比对的起点 3 2 43 2 4准全局序列比对准全局序列比对 所谓准全局比对就是在评价序列比对时不计 终端空位 end space 的得分或代价 在 下面的讨论中 以字符 代表 空位 表示插入或删除操作 我们仍然采用基本的比对算 法 但是 改进对序列两端空位的打分函数 终端空位是那些出现在序列第一个字符前或最后一个字符之后的空位 例如 图3 11 a 第二条序列中所有的空位都是终端空位 这两条序列长度相差较大 一条序列长度为8 而另一 条序列的长度为18 如果两条序列的长度相差比较大 那么 在最终的比对中一定存在很多空位 从而引起得分下降或代价增大 然而 如果忽略这些终端空位 则图3 11 a 所示的比对相当 好 6个匹配 1个失配 1个空位 但是 根据基本的序列比对算法 图3 11 a 的比对不是最优的 最优的比对似乎应该是 图3 11 b 所示的比对 因为其得分更高 同样多的空位 但除了空位以外全是匹配 尽管具 有更高的得分 而且第二条序列的字符全部匹配 但是 从寻找相似区域的观点来看 图3 11 b 所示的并不是我们所感兴趣的比对 为了匹配更多的字符 第二条序列被严重地割裂开 如果我 们在一段长序列中寻找一段与另一个短序列区域相似的区域 无疑图3 11 a 所示的排列更好 假设有两条序列 s 和 t 在这两条序列的比对中 s 后面的空位与 t 的后缀匹配 如果去掉这 部分比对 仅保留 s 与 t 前缀的比对 则被去掉的比对不再计分 因此 为了得到不计 s 右端空 位得分的最优比对 我们需要寻找 s 和 t 前缀的最优比对 由于在基本的序列比对算法中 矩阵 元素di j代表0 s i和0 t j的相似性 因而 可以取最后一行的最大值 取最后一行的最大值作为 比对的得分 就相当于不计序列 s 右端的空位得分 负分 注意 这里矩阵元素di j代表0 s i和0 t j 的最大相似性 并且计算公式就是 3 6 式 按照上述方法 我们可以同样取最后一列的最大值作为比对的得分 从而达到不计序列 t 后面空位的目的 进一步将两者结合起来 形成一种不计两条序列末尾空位的最优比对方法 现在 我们再讨论序列的前端空白 情况与末端空白相似 只要在设置初始条件时将矩阵 第0行和第0列元素的值设为 0 即达到不计序列前端空位得分的目的 综合上述分析结果 对两条序列进行准全局比较的算法如下 将矩阵的第一行和第一列元 素设置为 0 按行或列优先的次序计算矩阵 D 所有元素di j 0 i m 0 j n 取最后一行 和最后一列的最大值 而在搜索最优路径时 从矩阵最后一行和最后一列的最大值出发 按照基 本算法中的方法返回 准全局比较算法与基本算法在计算di j时的区别归纳为下列四个方面 1 第一行初始值为 0 表示不计第一条序列的前端空位 2 寻找最后一行的最大值 表示不计第一条序列的末端空位 3 第一列初始值为 0 表示不计第二条序列的前端空位 4 寻找最后一列的最大值 表示不计第二条序列的末端空位 对于最后一行和最后一列的另一种处理办法是 最后一行的横向移动不被空位罚分 最后 一列的纵向移动也不被罚分 这样 就可以允许在两条序列终端自由存在空位 当矩阵 D 所有元 素计算完以后 其右下角的值即为两条序列准全局比对的得分 反推路径即可得到序列的最佳比 对 3 2 53 2 5关于连续空位的问题关于连续空位的问题 我们定义 k 阶空位 是一个具有连续 空位 字符的区域 其空位字符的数目 k 1 对于 序列的突变 k 阶空位出现的可能性要大于 k 个不连续的空位 这是因为一个 k 阶空位对应于一 串字符的插入或者删除 对应于一个遗传突变事件 不连续的空位可能对应于多个不同的突变事 件 而一个突变的发生比多个突变同时发生的可能性要大 从这个意义上说 我们希望将 k 阶空 位与 k 个孤立空位在打分方面区别开来 k 阶空位的得分 或代价 应该比 k 个孤立空位的得分 高 或代价低 另外 在将 cDNA 与基因 DNA 比较时 通过改进的空位打分方式可以正确处理内含子 真 核生物的基因是非连续的 在编码区域的外显子之间插入了内含子 在基因转录过程中 内含子 被剔除 外显子连接起来形成完整的编码区域 并表现为 mRNA 通过反转录 根据 mRNA 形 成 cDNA 因此 同一个基因的 cDNA 序列与 DNA 序列的差别主要反映在内含子的删除 这样 在根据基因的 cDNA 寻找 DNA 时 我们也希望对大片段连续空位的处理与孤立空位的处理有所 区别 然而 直到现在 我们在算法中还没有区分单个空位及连续的多个空位 这意味着 k 阶空 位所带来的得分是一个线性的函数 设 p k 代表空位得分函数 其中 k 是连续空位的个数 则 p k kb 3 15 其中 b b 0 是单个空位得分的绝对值 对应于罚分 b 在本节中 将讨论解决序列比对连续空位问题的算法 该算法的时间复杂性是 O n3 该 算法与基本算法类似 主要差别在于得分函数不具有加和性 即不能将一个比对分成两个部分 而希望总的得分是两个部分的和 然而 当比对的分割点跨越 块 Block 边界时 得分的 加和性依然成立 任何一个比对都可以被唯一地分为若干个相继的块 有三类块 它们各自的组 成如下 1 两个字符的比对 2 与序列 s 的空位比对的 t 的最大连续字符序列 3 与序列 t 的空位比对的 s 的最大连续字符序列 最大 意味着不能再扩展 图3 12显示了一个序列比对及块的分割 上述第一类块的得分是 p a b a b 是两个比对的字符 其他两类的得分是 p k k 为连续的空位数 现在 对于一个比对的打分不再是按照单个比对的列进行 而是按块进行 只有当打破块 的边界时 得分的加和性才成立 对