计算DNA序列模式特征的匹配算法

戴胜冬，杨　昆

(杭州电子科技大学计算机学院，浙江 杭州 310018)

摘要：分析了DNA序列特征计算过程中的特殊性，提出了一种基于“空间换时间”的模式匹配算法，设计了以map数据结构来存储中间结果的方案，使得扫描DNA序列一次即可同时计算所有元组模式在该序列中出现的次数。实验结果及分析表明，算法提升了DNA序列模式特征计算的效率，较好地解决了计算DNA序列模式特征的问题。

关键词：算法；生物信息学；DNA序列；模式匹配

DOI: 10.13954/j.cnki.hdu.2015.01.018

计算DNA序列模式特征的匹配算法

戴胜冬，杨昆

(杭州电子科技大学计算机学院，浙江 杭州 310018)

摘要：分析了DNA序列特征计算过程中的特殊性，提出了一种基于“空间换时间”的模式匹配算法，设计了以map数据结构来存储中间结果的方案，使得扫描DNA序列一次即可同时计算所有元组模式在该序列中出现的次数。实验结果及分析表明，算法提升了DNA序列模式特征计算的效率，较好地解决了计算DNA序列模式特征的问题。

关键词：算法；生物信息学；DNA序列；模式匹配

DOI:10.13954/j.cnki.hdu.2015.01.018

收稿日期：2014-06-26

基金项目：国家自然科学基金资助项目(60903086)

通信作者：

作者简介：戴胜冬(1990-)，男，安徽巢湖人，在读研究生，表观遗传数据分析和挖掘.杨昆副教授，E-mail: yangkun@hdu.edu.cn.

中图分类号：TP311.1

文献标识码：：A

文章编号：：1001-9146(2015)01-0088-05

Abstract:This paper analyzed the specificity of DNA feature calculation， proposed an algorithm based on “space for time” idea， which designed a scheme of using map data structure to store the intermediate calculation results， so that each DNA sequence needs to be scanned just once to calculate all pattern features. Experiment results and analysis show that the algorithm can effectively improve the efficiency of calculating DNA feature， and better solve the problems common pattern calculating algorithms have in calculating DNA feature.

0引言

DNA序列特征的计算在生物信息学的研究中具有着重要的基础性地位，具体的序列特征有很多类型，而元组模式是其中一个重要且常见的类别[1-4]。对DNA甲基化状态预测就需要计算DNA序列中元组模式出现的次数[5-10]。例如，文献[9]通过计算人类淋巴细胞21号染色体的甲基化数据特征，发现序列模式是指示CpG岛甲基化状态的三类特征之一；文献[10]利用乳腺癌高甲基化基因的上游CpG岛的模式特征数据，分析得出序列特征可以作为标识来识别癌症高甲基化基因。然而如何高效率的量化DNA序列中元组模式特征已经成为生物信息学研究中的一个基本问题。在对DNA序列进行元组模式特征的计算中，DNA序列的数目往往非常庞大，单个DNA序列的长度也比较长，同时需要量化的元组模式数量较多、组合多变，这使得传统的模式匹配算法在DNA序列特征计算中时间效率不理想。

本文分析了上述DNA序列特征计算中存在的问题，针对其特殊性，提出了基于“空间换时间”思想的DNA序列模式计算的Map-based算法。算法设计了以元组模式为键元组模式在序列中出现次数为值的map数据结构，利用map数据结构来存储中间结果，使得扫描DNA序列一次即可同时计算所有元组模式在该序列中出现的次数。同时，设计了一个直观的单序列顺序遍历算法。作为对比参照，实现了常规的Boyer-Moore、KMP模式匹配算法，并在不同规模的真实数据上测试算法的性能。实验结果和分析表明，本文提出的Map-based算法在不同数据规模下都有着最好的性能。

1问题的描述和相关算法

1.1　问题的描述

DNA序列形式上是由A(腺嘌吟碱基)、T(胸腺嘧啶碱基)、C(胞嘧啶碱基)、G(鸟嘌吟碱基)4个碱基字符随意组合而成的一定长度的字符串。元组模式则是一定长度的DNA IUPAC code组成的片段，DNA IUPAC code中不仅包含A、T、C、G 4个碱基字符，还包含一些通配符，如R表示A or G，本文设计的算法不考虑其中的gap字符。将片段长度为n个碱基字符的元组模式称为n元组，例如元组模式GCCTRG长度为6个字符，称之为6元组。DNA序列的元组模式特征的量化，就是计算元组模式在DNA序列中出现的次数，例如DNA序列ATGCTGATGCC，元组模式ATGC在其中一共出现了2次，则DNA序列ATGCTGATGCC的元组模式ATGC的特征量化后即为2。实际计算过程中，要同时计算的DNA序列是由数量达到成千上万的DNA序列组成的序列集，元组模式特征则是由数量几百的元组模式序列组成的模式集。

1.2　相关的算法

KMP算法[11]是一个常用且高效率的模式匹配算法，该算法在常规顺序扫描匹配模式算法的基础上，预先构造模式本身的“失配推进表”，每当一趟匹配过程中出现字符比较不等时，不需要回溯扫描指针，而是利用已经得到的“部分匹配”的结果，实际上就是查询“失配推进表”得到在模式当前位置失配的情况下模式可以向右“滑动”的最远的一段距离，继续进行比较，此算法可以在序列长度为m，元组模式长度为n时在O(n+m)的时间数量级上完成串的模式匹配操作。

另一个久负盛名的模式匹配算法是BM(Boyer-Moore)算法，它是一种非常高效的字符串搜索算法。它由Bob Boyer和J Strother Moore设计于1977年[12]。此算法仅对搜索目标字符串(关键字)进行预处理，而非被搜索的字符串。虽然Boyer-Moore算法的执行时间同样线性依赖于被搜索字符串的大小，但是通常仅为其它算法的一小部分，不需要对被搜索的字符串中的字符进行逐一比较，而会跳过其中某些部分。通常搜索关键字越长，算法速度越快。

2Map-based算法

综合DNA元组模式特征计算中的特殊性，考虑采用在存储和查询上时间效率都较为优秀的Map来记录和查询DNA序列的元组模式特征及数量，在含有通配符的元组模式扩展阶段，Map被用来存储展开模式集。本文提出的算法如下描述。

算法1MAP-BASED-COMPUTE-PATTERN-COUNT(S，P，K)

输入：序列集合S，包含了所有的待计算的DNA序列。元组模式集合P，包含了所有的元组模式。元组模式长度的集合K，如需要计算的元组模式包含了2元组、3元组、4元组、6元组，则K={2，3，4，6}。

1letcount-mapbe a map

2letexpanded-patterns-mapbe a multimap

3letwc-patterns-mapbe a map

4letnowc-patterns-mapbe a map

5for each patternp∈P

6ifpcontains wildcards

7mapcount-map[p]to 0

8for each patternep∈expanded patterns ofp

9mapwc-patterns-map[ep]to 1

11mapexpanded-patterns-map[ep]top

12else then

13mapcount-map[p]to 0

14mapnowc-patterns-map[p]to 1

15for each sources∈S

16fori=0 tos.length

17for eachk∈K

18sub-string=SUB-STRING(s，i，k)

19ifwc-patterns-map[sub-string]==1

20for eachkeyinexpanded-patterns-map[sub-string]

21count-map[key]=count-map[key]+1

22ifnowc-patterns-map[sub-string]==1

23count-map[sub-string]=count-map[sub-string]+1

本算法借助map来存储预处理阶段的中间结果，一方面合理且高效地解决了包含通配符的模式计算问题，另一方面借助map较高的查找效率，可以快速地完成模式的定位和模式的计数操作。程序在预处理阶段构造的map虽然占用了较多的内存空间，但是在DNA序列的模式计算中，通过map的内存空间换取更快的查找效率来节省计算任务的时间是完全可取的。

3实验

本文实验采用真实的DNA序列数据，该数据来自于预测人类基因甲基化程度的原始数据集(http://rulai.cshl.edu/HDMFinder/M-set.fa)[13]，其中DNA序列的平均长度为15 370个碱基字符，原始数据未进行其他处理。元组模式来自文献[3]中的序列特征模式motif: 1个二元组，64个三元组，256个四元组，133个六元组，其中有30个包含通配符的六元组，共计454个元组模式。实验的软件条件为Java 7，Windows 7 32-bit，硬件条件为Intel CPU i5 3.2 GHz双核，内存2.0 GB。

对本文中实现的4种模式计算算法，分别在DNA序列个数为100、200、400、800、1 600、3 200的规模下测试算法的效率，在每个具体的数据上重复运行算法5次，统计每个算法运行时间的平均值和标准差，实验结果如表1所示。

对比表1中数据可以看出，运行速度最慢的是BM算法，最快的是Map-based算法，Map-based算法比第二快的单个序列顺序遍历算法平均快2.45倍。同时可以看出，实验数据中标准差都一致的较小，相对平均值最大的占2.84%，表明实现的4个算法性能稳定，实验结果可靠。同时分析各算法耗时的方差可以看出，普遍上Map-based算法方差都远远低于其他模式匹配算法，说明Map-based效率上的稳定性更好。

表1　4个算法在不同数据上的执行时间　ms

从实验结果可以看出来，常规的BM和KMP算法在本文实验条中表现非常不佳，这主要是因为BM和KMP算法比较适用于模式比较复杂、长度较长情况下单个序列的模式计算，不适合序列模式个数多、长度短、复杂度低的情况。在算法开始匹配之前，两者都要对模式进行预处理，得出一些辅助数据来加速算法运行，在序列模式计算中这个预处理过程要进行m×n次，m、n分别为序列集合和元组模式集合的元素个数，算法所需的辅助数据的重复计算会给性能带来很大的损失，模式匹配上的加速不足以弥补模式预处理消耗的时间。

为了比较4种算法对常规任务的计算效率，即对单个长度较长的复杂模式的计算效率，基于原始数据集，随机截取长度为200的5个DNA片段分别作为5个单一复杂模式，在长度为125 000、250 000、500 000、1 000 000、2 000 000、4 000 000的DNA序列上应用本文的4种算法来分别计算每个复杂模式，统计每个算法在5个单一模式的计算上耗费时间的平均值和标准差，实验结果如表2所示。

表2　4个算法在不同序列长度上单个模式特征的计算时间　ms

从表2数据可以看出，在上述实验条件下Map-based算法速度最慢，平均比BM、KMP、单个序列顺序遍历算法分别慢47、46、97倍。这个结果说明，模式较为复杂、模式长度较长的情况下，对于单个模式计算上常规算法具有优势。综合表1、表2中的结果可以得出，相较于常规算法，本文提出的Map-based算法能更好地适用于DNA序列模式的计算。

4结束语

本文提出的算法在应用到DNA的序列特征计算时，能很好地应对DNA序列数目庞大、单个DNA序列的长度较长，需要计算的元组模式数量多、组合多变的情况，实验结果表明，在实际应用中该算法不仅效率更高而且性能稳定性也较高，并且易于实现，具有一定的应用推广价值。

参考文献

[1]Fan S C， Zhang X G. Progress of bioinformatics study in DNA methylation[J]. Prog Bio-chem Biophys，2009，36(2):143-150.

[2]戴胜冬，杨昆，顾靓，等.DNA序列的甲基化特征提取软件[J].生物信息学，2012，10(1):5-7.

[3]杨昆，张彦斌，戴胜冬，等.DNA甲基化的重要特征[J].生物物理学报，2012，28(11):910-922.

[4]Edwards J R， O’Donnell A H， Rollins R A， et al. Chromatin and sequence features that define the fine and gross structure of genomic methylation patterns[J]. Genome Research，2010，20(7):972-980.

[5]Fang F， Fan S， Zhang X， et al. Predicting methylation status of CpG islands in the human brain[J]. Bioinformatics，2006，22(18):2 204-2 209.

[6]Bhasin M， Zhang H， Reinherz E L， et al. Prediction of methylated CpGs in DNA sequences using a support vector machine[J]. FEBS Letters，2005，579(20):4 302-4 308.

[7]Fan S， Zhang M Q， Zhang X. Histone methylation marks play important roles in predicting the methylation status of CpG islands[J]. Biochem. Biophys. Res. Commun.，2008，374(3):559-564.

[8]凡时财，邹见效，徐红兵，张学工.人类基因组CpG岛甲基化概况的预测[J].科学通报，2010，55:1 329-1 334.

[9]Bock C， Paulsen M， Tierling S， Mikeska T， Lengauer T， Walter J. CpG island methylation in human lymphocytes is highly correlated with DNA sequence， repeats， and predicted DNA structure[J]. PLoS Genetics，2006，2(3):e26.

[10]Lv J， Su J， Wang F， et al. Detecting novel hyper methylated genes in breast cancer benefiting from feature selection[J]. Computers in biology and medicine，2010，40(2):159-167.

[11]Knuth D E， Morris， Jr J H， Pratt V R. Fast pattern matching in strings[J]. SIAM journal on computing，1977，6(2):323-350.

[12]Boyer R S， Moore J S. A Fast String Searching Algorithm[J]. Communications of the ACM，1977，20(10):762-772.

[13]Das R， Dimitrova N， Xuan Z， et al. Computational prediction of methylation status in human genomic sequences[J]. Proceedings of the National Academy of Sciences，2006，103(28):10 713-10 716.

Pattern Matching Algorithm for Calculating DNA Feature

Dai Shengdong， Yang Kun

(SchoolofComputer，HangzhouDianziUniversity，HangzhouZhejiang310018，China)

Key words: algorithm; bioinformatics; DNA sequence; pattern matching