DNA匹配的若干算法研究及分析

摘要：

对于DNA序列的查找、标记分析其基因特性是一个热门话题。在长度数以兆计的DNA序列中通过比较查找许多特征子序列，需要高效的匹配算法。本文总结了多种DNA匹配算法，并比较其优劣。

问题描述：

在DNA长序列间的比较能够发现隐藏的基因信息，而在DNA中查找特征序列能帮助发现病变，这些都是很有用方向。比如通过检测癌变基因来发现癌症，通过发现固定结构序列来查找功能与基因的对应关系。

DNA是一个一维的化学元素序列，可以把它看成是一个巨大的线性磁带，元素是{A,C,G,T}。不同的元素组合表征了不同的化学结构。一个有机体的序列特征能够表征它的特点。如果我们比较两个序列X,Y。已知X表示某种有害特征，如果Y与X像类似，我们就能猜测Y也是表示有害信息。其中一个研究目的就是检测某人的基因中是否包含不良基因片断。[1]

匹配问题：

问题一：精确匹配

在X中查找是否包含Y字串。

最简单而且保证能确保结果的做法是穷举法。

算法描述[2]如下：

Begin initialize A,x,text,n <- length[text], m <- length[x]

s <- 0

while s <= n – m

if x[1..m] = text[s+1..s+m]

then print “pattern occurs at shift” s

s <- s+1

return

end

但是这种方法的效率很低，时间复杂度是O((n-m+1)m)。

进而我们可以通过Boyer-Moore算法提高匹配效率。

Begin initialize A,x,text,n <- length[text], m <- length[x]

F(x) <- last-occurrence function

G(x) <- good-suffix function

s <- 0

while s <= n – m

do j <- m

while j > 0 and x[j] = text[s+j]

do j <- j-1

if j = 0

then print “pattern occurs at shift” s

s <- s+G(0)

else s <- s+ max[G(j), j-F(text[s+j])]

return

end

但是在现实应用中，单单精确匹配是不够的，更多的时候我们需要的是带有容错性的匹配算法。为了更好的了解模糊匹配算法，我们先要给出以下几个定义：

定义一：编辑距离[2]

字串X,Y的编辑距离是，X经过多少步基本操作能转变为Y。

基本操作包括：

替换：X的一个字符被替换成Y中的一个字符。

插入：Y的一个字符插入到X，因此X的长度增加一个字符。

删除：X的一个字符被删除，因此X的长度减少一个字符。

定义二：任意匹配字符

如果一个字串包含字符，表示这一位上的字符与任意字符匹配。

定义三：容错匹配

允许匹配的两个字串有一定数目的字符是不匹配的。错误量可以用编辑距离来表示。

其中，如何计算编辑距离的方法如下：

通过动态优化的方法：

OPT(i,j) = min[_xiyi + OTP(i-1, j-1),  + OPT(i-1, j),  + OPT(i, j-1)]

N	8	6	5
A	6	5		5	5
E	4		2	4	4
M	2		3	4	6
-		2	4	6	8
	-	N	A	M	E

其中中表示最佳匹配。

但是这种最佳算法的时间复杂度也只能是O(mn)。

还有其他的算法来寻找包含替换和间隔的字串对齐[3]，但是这些方法需要有特殊硬件的支持。

最近的研究，提出了若干高效的方法：

其中一个是SAX( Symbolic Aggregate approXimation)。主要思想是通过把时序序列转换成字符串，然后通过统计规律，快速的找到序列中的突变片断。

比如：

对于像心电图这样规律的时序信号，在突变的情况少，通过ＳＡＸ方法能够节省许多计算量，而且非常简单。

简化运算的推理过程如下：

通过一系列近似计算，把粒度变大，简化计算量，但同时能够有很好的精确度。

下面介绍如何把时序信号转变成字符串。

方法就是选择合适的区间，把每个区间段中的信号映射到一个字母上。如下图：

然后统计字符串序列的频率，转换成位图

比如：对于ＤＮＡ序列：

CCGTGCTAGGGCCACCTACCTTGGTCCGCCGCAAGCTCATCTGCGCGAACCAGAACGCCACCACCTTGGGTTGAAATTAAGGAGGCGGTTGGCAGCTTCCAGGCGCACGTACCTGCGAATAAATAACTGTCCGCACAAGGAGCCGACGATAAAGAAGAGAGTCGACCTCTCTAGTCACGACCTACACACAGAACCTGTGCTAGACGCCATGAGATAAGCTAACA

频率分别是

C: 0.20

A: 0.26

T: 0.24

G: 0.30

制作成位图:

同样的，对于两两字母组合，甚至更多子串的组合，也能做出位图

然后，根据频率的不同，染色：得到下列位图

这时候，比较两个时序序列的异同，只要比较两个位图就可以了。因为对于相似的序列，他的位图必然是相近的。

另外一个方法是LSH（locality-sensitive hashing）。LSH是由Indyk 和 Motwani在1998年提出，然后由Gionis 等人针对高维几何问题优化。目的是把字符串的子串匹配简化到一个更容易控制的精确匹配问题。

原理如下：

有两个字符串s1,s2，他们的长度都是d,字符集是Σ。对于固定的r <>如果s1,s2最多只有r个不同子字符，我们称s1,s2相似。

为了检测两个字串是否相似，我们构造以下随机过滤器。从集合{1...d}中选择k个位置i_1,i₂,...i_k；为了简化计算，位置可以重复。我们定义函数f，Σ ^d->Σ ^k

f(s) = <>1], s[i₂], ...s[i_k]>

f被称为LSH函数。我们认为两个字串相似当且只当f(s1)= f(s2)

对于真正相似的字串，他们满足f(s1)=f(s2)的概率 > (1-r/d)^k

所以对于相似的字串，很容易被识别出来，对于不相似的字串，也有很大的概率分辨出来。

但是也可能把不相似的字串识别成相似，因为刚好满足f(s1)=f(s2)条件。

要判断两个DNA是否相似，只要设定好k,然后通过f函数，就能很快的算出结果。

当然，除了SAX,LSH还有很多其他方法。其中这两个方法是我觉得最有特色的，和平常的思路不一样。也觉得能够在许多其他领域上应用。

参考文献：

[1] Algorithm design. Jon Kleinberg, Eva Tardos Pearson education 2006

[2] Pattern Classification. Richard O. Duda, Peter E. Hart, David G. Stork. 2004

[3] Huang,x. and Miller, W. (1991) A time-efficient, linear-space local similarity algorithm. Adv. Appl. Math., 12, 337-357

[4] E. Keogh, J. Lin and A. Fu (2005). HOT SAX: Efficiently Finding the Most Unusual Time Series Subsequence. In Proc. of the 5th IEEE International Conference on Data Mining (ICDM 2005), pp. 226 - 233., Houston, Texas, Nov 27-30, 2005.

[5] Li Wei, Eamonn Keogh and Xiaopeng Xi (2006) SAXually Explict Images: Finding Unusual Shapes. ICDM 2006.

[6] Li Wei and Eamonn Keogh (2006) Semi-Supervised Time Series Classification. SIGKDD 2006

[7] Xiaopeng Xi, Eamonn Keogh, Christian Shelton, Li Wei & Chotirat Ann Ratanamahatana (2006). Fast Time Series Classification Using Numerosity Reduction. ICML.

[8] Kumar, N., Lolla N., Keogh, E., Lonardi, S. , Ratanamahatana, C. A. and Wei, L. (2005). Time-series Bitmaps: A Practical Visualization Tool for working with Large Time Series Databases . In proceedings of SIAM International Conference on Data Mining (SDM '05), Newport Beach, CA, April 21-23. pp. 531-535

[9] Ratanamahatana, C. A. and Keogh. E. (2004). Everything you know about Dynamic Time Warping is Wrong. Third Workshop on Mining Temporal and Sequential Data, in conjunction with the Tenth ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD-2004), August 22-25, 2004 - Seattle, WA.

[10] Mayur Datar Nicole Immorlica Piotr Indyk Vahab S. Mirrokni Locality-sensitive hashing scheme based on p-stable distributions Proceedings of the twentieth annual symposium on Computational geometry 2004

[11] Similarity search in high dimensions via hashing A Gionis, P Indyk, R Motwani - Proc 25th VLDB Conference, 1999

[12] Jeremy Buhler Ecient Large-Scale Sequence Comparison by Locality-Sensitive Hashing Bioinformatics 17(5) 419{428, 2001

[13] C Yang MACS: Music Audio Characteristic Sequence Indexing for Similarity Retrieval
- IEEE Workshop , 2001

[14] Piotr Indyk Rajeev Motwani Approximate nearest neighbors: towards removing the curse of dimensionality ACM 1998

[15] Jeremy Buhler Efficient large-scale sequence comparison by locality-sensitive hashing BIOINFORMATICS Vol. 17 no. 5 2001 Pages 419–428