序列比对（待续）_双序列比对程序

一、双序列比对

动态规划（Dynamic programming）

动态规划 (Dynamic programming, DP) 是通过把原问题分解为相对简单的子问题的方式求解复杂问题的方法。动态规划常常适用于有重叠子问题和最优子结构性质的问题。

• 重叠子问题：自底向上
  下面三个函数分别使用了递归、递归+列表、迭代方法求解斐波那契数列。后两种方法仅对重叠子问题求解一次，因而时间复杂度急剧下降。其中迭代呈现的自底向上策略就是动态规划对重叠子问题的处理思路。

#f(n)=f(n-1)+f(n-2); f(0)=0; f(1)=1

#directly, inefficient, O(2^n)
def f(n):
    if n==0:
        return 0
    elif n==1:
        return 1
    else:
        return f(n-1)+f(n-2)

#top-down, O(n)
def f(n):
    memo=[-1]*(n+1)
    if n == 0:
        return 0
    if n == 1:
        return 1
    if memo[n] == -1:
        memo[n]=f(n-1)+f(n-2)
    return memo[n]

#down-top, O(n)
def f(n):
    memo=[-1]*(n+1)
    memo[0]=0
    memo[1]=1
    for i in range(2,n+1):
        memo[i]=memo[i-1]+memo[i-2]
    return memo[n]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

• 最优子结构
  简单理解就是问题的最优解是子问题最优解的组合。以比对为例，碱基之间的比对只有三种情况：相同、替换和空位。在使用替换矩阵和空位罚分对三种情况打分后，得分最高者即是该位置的最优比对（公式1,2）。在对所有可能比对情形打分后，即可获得最优比对。由此即可构建下述的双序列比对算法：N-W和S-W。更为详细的推演过程可参考coursera课程：《生物信息学: 导论与方法》。

全局比对（Needleman-Wunch）

全局是指对两条序列所有残基进行比对。分别求三种比对状态的得分，最大值即为该位置最优比对。通过回溯，得到最佳比对结果。
F ( i , j ) = { F ( i − 1 , j − 1 ) + s ( x i , y j ) , x i : y j F ( i − 1 , j ) + d , x i : − F ( i , j − 1 ) + d , y j : − F(i,j)=

$$\begin{cases} F(i-1,j-1)+s(x_i,y_j),\quad &x_i : y_j \\ F(i-1,j)+d, &x_i : - \\ F(i,j-1)+d,&y_j : - \end{cases}$$ F(i,j)=⎩⎪⎨⎪⎧​F(i−1,j−1)+s(xi​,yj​),F(i−1,j)+d,F(i,j−1)+d,​xi​:yj​xi​:−yj​:−​
公式1. Needleman-Wunch 算法的状态转移方程

局部比对（Smith-Waterman）

相较全局比对，局部比对状态转移方程中引入了一个最小值0，这意味着由于回溯过程可以被最小值中断，及比对是局部的。这更符合只有部分片段相似的生物学现象·，如结构域、内含子等。
F ( i , j ) = { F ( i − 1 , j − 1 ) + s ( x i , y j ) , x i : y j F ( i − 1 , j ) + d , x i : − F ( i , j − 1 ) + d , y j : − 0 F(i,j)=

$$\begin{cases} F(i-1,j-1)+s(x_i,y_j),\quad &x_i : y_j \\ F(i-1,j)+d, &x_i : - \\ F(i,j-1)+d,&y_j : - \\ 0 \end{cases}$$ F(i,j)=⎩⎪⎪⎪⎨⎪⎪⎪⎧​F(i−1,j−1)+s(xi​,yj​),F(i−1,j)+d,F(i,j−1)+d,0​xi​:yj​xi​:−yj​:−​
公式2. Smith-Waterman算法的状态转移方程

---

二、一对多比对

无论是数据库中存储、还是NGS产生的序列都具有数量庞大的特点。因而对于数据库搜索和NGS短序列比对问题，单纯使用双序列比对方法的开销巨大。为了降低开销到可用水平，此类问题一般先根据索引将query或reads定位到大致区间，然后再使用基于动态规划的算法进行比对。

基于哈希表的seed方法

• 原理
  将query分为长度固定的小片段（seed），以seed为键在目标序列（已转换成哈希表）进行查询，得到匹配的锚点。向锚点两侧进行延伸打分，保留符合条件的比对。
• BLAST
  https://zhuanlan.zhihu.com/p/54294281
• FASTA
  … …
• MAPQ

基于后缀的FM-index算法

• 数据结构

表1. 基于后缀的数据结构

• BWT
  下面的代码展示了 Burrows-Wheeler 转换的基本概念，包括Burrows-Wheeler matrix, Burrows-Wheeler transform等。

#burrows-wheeler transform
#rotation
str+='$'
matrix=[]
for i in range(len(str)):
    str=str[1:len(str)]+str[0]
    matrix.append(str)
#sort according right context
bw_matrix=sorted(matrix)
#bwt,i.e. last column
last_column=[ i[-1] for i in bw_matrix]
lc=''.join(last_column)
1
2
3
4
5
6
7
8
9
10
11
12

由于第一列和最后一列的相邻关系（Last/First Mapping），可以将bwt逆转换，得到原始字符串。以下代码从概念上展示了这一过程。

#first column and count
first_column=[ i[0] for i in bw_matrix]
fc=''.join(first_column)
count={}
[count.setdefault(k,fc.count(k)) for k in set(fc)]
#last column subscript array
sub_D={'a':0,'c':0,'g':0,'t':0,'$':0}
sub=[]
for i in lc:
    sub_D[i]+=1
    sub.append(sub_D[i])
#mapping
def fc_index(c):
    global idx
    if c == 'a':
        idx=0
    elif c == 'c':
        idx=count['a']
    elif c == 'g':
        idx=count['a']+count['c']
    elif c == 't':
        idx=count['a']+count['c']+count['g']
    else:
        pass
    return idx

ch=lc[0];ori=ch
i=0
n=0
while n < len(lc)-2:
    print(ch)
    print(i,fc_index(ch),sub[i])
    i=fc_index(ch)+sub[i]
    print(i,fc_index(ch),sub[i])
    ch=lc[i]
    ori+=ch
    n+=1

ori=ori[::-1]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

下图总结了上述的BWT算法过程。在了解BWT的概念后，紧接着的一个问题就是如何将之用做索引呢？

图1. Burrows -Wheeler Transform

• FM-index
  下面以序列tgc为例说明索引机制。由于BWT的F列是经过排序的，方便查找，因而首先在F列查找字符c，可以定位到3,4行。然后在L列3,4行寻找字符g，由于只有1个g，g定位在F列第5行。以此类推，将t定位在第6行。（图2. A-C）
  那么第6行对应于字符串的哪个位置呢？显然无法从bwt数据结构直接得到，因而又引入了后缀数组作为辅助（图2. D），但从表1中可以看到后缀数组数据量较大。为了减少数据量，可以只保存部分数组值，比如只保留0,2,4。如果想知道t的位置，根据FL mapping关系可知t的前一个字符为a，而a对应的数值为0，可得t的位置为1。
  此外还有一个不容忽视的问题，即L列搜索效率问题。不同于已经排序的F列，如果仅仅使用迭代方法扫描L列，那么所用的时间复杂度可粗略估算为O(N^2/16)。FM-index使用一个记录字符出现次数的数组（图2.E; bwa中也称为occ数组）来解决该问题。根据预先计算好的occ数组，可以以O(1)复杂度得到某一范围内的特定字符出现的次数和次序。和SA类似，occ也只储存部分值（被储存值称为checkpoint），而未被储存的值可根据附近的checkpoint计算得到。
  以上F, L, SA, checkpoint 就是FM-index 的基本数据结构，对人类基因组来说，其总体大小在1G左右。
  

图2. FM-index

• 软件
  根据陈凤珍等[3]对四款常见软件的测评结果，①比对率：Bowtie2 ≈ bwa > SOAP2 > MAPQ；②用时：MAPQ >> bwa > Bowtie2 > SOAP2。

---

三、多序列比对

… …

---

参考
[1]【陈巍翻译】A42：基因组研究中用到的索引
[2]【陈巍翻译】A44：FM索引
[3] 陈凤珍,李玲,操利超,严志祥.四种常用的生物序列比对软件比较[J].生物信息学,2016,14(01):56-60.