序列比对——DNA双序列比对1. 基本算法2. BestScore 空间改进算法3. 分治法-线性空间改进算法4. k-band 时间改进算法5. 仿射罚分多序列比对1. SP measures2. 星比对算法3. 多核并行测试数据集mitochondrial genomes16s rRNACOVID-19RNA模拟数据模拟noncoding数据

序列比对——DNA

将两个或多个序列排列在一起,标明其相似之处。序列中可以插入间隔(通常用短横线“-”表示)。对应的相同或相似的符号(在核酸中是A, T(或U), C, G,在蛋白质中是氨基酸残基的单字母表示)排列在同一列上。这一方法常用于研究由共同祖先进化而来的序列,特别是如蛋白质序列或DNA序列等生物序列。在比对中,错配与突变相应而空位与插入或缺失对应。序列比对还可用于语言进化或文本间相似性之类的研究。

双序列比对

1. 基本算法

(1) 动态规划

篇幅有限不做展开,网络上已经有很多详细的资料,也可参考一下我写的一篇笔记

(2) 打分机制

  1. match 匹配 +1
  2. mismatch 不匹配 -1
  3. gap 缺失匹配 -2

根据打分机制写出动态规划函数,有序列s和t

具体代码可参考链接,如有错误,还望指正,感谢。

2. BestScore 空间改进算法

前文使用的算法,不难看出空间的复杂度是O(mn),m和n是序列的长度。

但是如果我们不需要回溯出具体比对方案,我们可以对其空间复杂度进行改进,改进后的空间复杂度为O(n),n为m和n中的最小值。因为不需要回溯具体路径,就只需要保存一行数据进行更新和迭代,即可算出最后的比分。

具体代码参考链接,如有错误,还望指正,感谢。

3. 分治法-线性空间改进算法

典型的时间换空间算法,对于长序列而言,电脑内存不够,可以采用此算法

(1) 算法思路

选取序列 的中间位点 , 那么在双序列比对 中 的对应只会有两种可能:

  1. ~

  2. ~ '-'

    注:j 值是遍历1..n取得分最高的值,通过不断分治,最后得到比对的结果

    具体代码参考链接,如有错误,还望指正,感谢。

4. k-band 时间改进算法

如果我们比对的是两个相似的序列,那我们可以采用k-band算法对其时间复杂度进行改进。

(1) 算法思路

可知两个相似的序列得到的最佳比对,其回溯路径一般在主对角线附近。那么就不需要对整个矩阵进行填充和计算,只需要对主对角线附近的区域进行计算,这个区域便称为k-band。

其中k值,代表从主对角线区域往外扩展的偏移量,我们将k值的初始值设置为两个序列长度的差值。k值的取值是不断迭代的,每次翻倍,直至矩阵获取的最大分值不再改变。

(2) 分析

采用k-band法对算法进行改进,算法的时间复杂度降为O(kn)。证明如下:

具体代码参考链接,如有错误,还望指正,感谢。

5. 仿射罚分

仿射罚分是对之间的打分机制进行了改进。因为在实际情况中,更倾向于连续的缺失对齐,而不是非常离散的缺失。所以将缺失对齐的首个gap分配更高的罚分,对后续连续的gap予以较低的罚分。

下图给出了一个例子,左边和右边的匹配按照之前的打分机制来看,所得到的score应该是一致的。但是实际上左边更符合实际的情况,因此引入了仿射罚分的机制,来尽可能的避免出现像右边的情况。

(1) 算法思路

仿射罚分还是利用动态规划的思想来解决,但是需要用到多状态转换,因为需要根据对齐过程中不同的状态进行打分。不难分析,在对齐的过程中一共有三个状态。那么便建立三个打分矩阵,对矩阵同时进行遍历求解,最后得到三个终值,其中最大的值,便为最佳对齐的score。

(2) 新打分机制

  1. match 匹配 +1
  2. mismatch 不匹配 -1
  3. gap 缺失匹配首个 -h
  4. gap 缺失匹配非首个 -g

(3) 三个状态

  1. a:s[i]~ t[j] 即序列s字符对齐序列t字符
  2. b:'—'~ t[j] 即序列s字符对齐gap,表示序列t中缺失
  3. c:s[i]~'—' 即gap对齐序列t字符,表示序列s中缺失

(4) 初始化三个矩阵的值

(5) 动态规划函数

具体代码参考链接,如有错误,还望指正,感谢。

多序列比对

1. SP measures

Sum of pairs (SP) measures : Scoring a alignment based on pairwise alignments.

SP measures 是一个基于双序列比对的多序列对齐打分机制

SP打分机制,是将多序列比对中的所有的双序列比对(PSA)都拿出来算一遍Score,最后再汇总加起来得到最终的得分。

打分函数

计算方式分为两种(注:p(-,-)=0):

  1. 对多序列对齐中的每一列得分进行计算,最后相加
  2. 对多序列对齐中的每一对PSA得分进行计算,最后相加

多序列比对仍然可以采用动态规划的方法来求解最佳对齐,但是其算法复杂度为,其中k为序列的个数,n为序列的长度。可以看到,多序列比对采用SP measures打分机制是一个NPC问题,只能采用启发算法求解一个近似最佳的对齐。

2. 星比对算法

算法实现思路

  1. 在序列中选择一个中心序列 ,中心序列的选择算法,在下文详细阐述
  2. 将中心序列和其余序列进行双序列比对,得到k-1个双序列比对
  3. 将k-1个双序列比对进行整合,整合原则“once a gap, always a gpa”

例子

  1. 找到中心序列

    遍历全部c值,得到每条序列的scores值,取得最大分值的序列做为中心序列

  2. 一一进行双序列比对

  3. 整合PSA,得到MSA

    具体代码参考链接,如有错误,还望指正,感谢。

3. 多核并行

由上文可知,建立得分矩阵选取中心序列是非常耗时的。因此在建立得分矩阵的这个环节,我们可以采取多核并行,加快程序的运行速度。由实际的运行速度来看,并行计算要比单核运行快3.6倍左右(根据电脑的核心数有关,测试电脑的核心数是4核)。

具体代码参考链接,如有错误,还望指正,感谢。

 

测试数据集

mitochondrial genomes

16s rRNA

COVID-19

RNA模拟数据

模拟noncoding数据