（七）R语言生物序列比对——Smith-Waterman局部比对算法

实现原理与代码

像全局序列匹配一样，局部序列比对的目的也是找到两个序列之间的相似度。 Smith-Waterman这一局部比对算法的原理和Needleman-Wunsch全局比对算法一样，也是运用了动态规划(DP)的思想。具体关于Needleman-Wunsch全局比对算法可参考我的这篇博客。Smith-Waterman算法经常两序列长度相差较大时使用，或需要匹配两序列中相似的片段。在的工作原理和Needleman-Wunsch也非常相似，只是在创建初始矩阵时第一行和第一列的初始值均设为0，进行DP时多了0这个选项，且最后回溯的时候是从整个矩阵中的最大值开始回溯到0。
该算法的第一步也是构造得分函数，具体如下：

match <- 2
mismatch <- -1
gap <- -2

scoreMat <- matrix(misMatch, 5, 5)
dimnames(scoreMat) <- c(c("A", "T", "C", "G", ""), c("A", "T", "C", "G", ""))
diag(scoreMat) <- match
scoreMat[5,] <- gap
scoreMat[,5] <- gap
#scoreMat如下
#   A  T  C  G  
#A  2 -1 -1 -1 -2
#T -1  2 -1 -1 -2
#C -1 -1  2 -1 -2
#G -1 -1 -1  2 -2
#  -2 -2 -2 -2 -2

在决定好得分后，程序会把两个需要比对的序列分别作为一个新矩阵的列和行标题。接着，矩阵的第一行和第一列所有元素初始值设置为0。接下来程序会从[2,2]开始遍历元素。在使用DP算法填满整个表后，会进行回溯，获得具体的匹配序列。
假设$a_{i,j}$为含需要对比序列的矩阵，其中元素的值的计算如下：
$$a[i,j]=max\{\begin{array}{l}a[i-1,j-1]+scoreMat[seq1[i-1],seq2[j-1]]\\ a[i,j-1]+scoreMat[5,5]\\ a[i-1,j]+scoreMat[5,5]\\ 0\end{array}$$
其中第一行$scoreMat[seq1[i-1],seq2[j-1]]$的值为该元素对应碱基的得分，第二行与第三行中的$scoreMat[5,5]$代表gap。
该部分代码实现如下所示：

alignmat <- matrix(0, length(seq1)+1, length(seq2)+1)
alignmat[,1] <- seq(0, gap*length(seq1), gap)
alignmat[1,] <- seq(0, gap*length(seq2), gap)
dimnames(alignmat) <- list(c("",seq1), c("",seq2)) #在每个序列前加入一个空字符
  
for(i in 2:nrow(alignmat)){ #因为[1,1]被填了，从[2,2]开始
  for(j in 2:ncol(alignmat)){
    alignmat[i,j] <- max(alignmat[i-1,j-1]+scoreMat[seq1[i-1],seq2[j-1]],
    					alignmat[i-1,j]+scoreMat[5,5],
    					alignmat[i,j-1]+scoreMat[5,5],
    					0) #和全局比对比起来多了个0的选项
  }
}

我们用序列AAAGAATTCAAA和GATT，我们就获得了一个如下所示的alignmat：

回溯的时候，要从表里最大的元素开始，回溯到值为0的元素。具体回溯过程和我的这篇博客中描述的全局比对算法过程相似。但不同的是，局部比对算法尽量保持了较短序列的连续性，将其作为一个整体与较长序列的部分进行比对。
排列好后的序列如下：

AAGAATTGTA
--GAAAT---

全部代码

由于该算法经过多次优化，笔者将不自行设计。请从该链接下载使用的代码，提取码：ngw8，文件原名为pairAlign.R，因为敏感词被屏蔽了，使用时请改回原名。以下的代码介绍的主要是使用方法：

install.packages("BiocManager") #调用包
BiocManager::install() #[a/s/n]输入a下载全部附属包
BiocManager::install("Biostrings") #install biostrings
suppressMessages(library(Biostrings)) #我们要用的得分矩阵要调用这个包

source("pairAlign.R") #无法载入的话使用完整文件路径
pairAlign("GATTACCTAGATTTCTAACGA", "GAATTCTAG", 
		nucleotideSubstitutionMatrix(2,-1), 
		gapOpening = -2, gapExtension = -2, type = "local")
		#前两个参数为需要比对的序列，蛋白质和DNA序列都可以
		#第三个参数为用到的得分矩阵。DNA可以用nucleotideSubstitutionMatrix，蛋白质可以用"BLOSUM50"
			#nucleotideSubstitutionMatrix(2,-1)定义了match得分为2，mismatch得分为-1
		#第四和第五个参数为gap的惩罚。默认的是开gap和延长gap的惩罚相同，但实际上开gap的惩罚应大于延伸gap的惩罚，因为现实中序列出现断裂的可能远小于两个断裂序列间距边远的可能
		#最后的type是序列比对的方法。"local"就是Smith Waterman局部比对

运行结果如下：

注意，这里匹配完成展示的只有长序列中相应的位置，而不是显示两个完整的序列。

结束语

在生物信息学中，获得两个的全局序列比对结果能够帮助我们推测两序列的相似度/同源性。下篇博客将会介绍多序列比对和BLAST的原理及使用方式，敬请期待！