一般間隙近似無重疊模式匹配

武優西,陳 彤,閆文杰,高雪冬

(河北工業大學 人工智能與數據科學學院,天津 300401) (河北省大數據重點實驗室,天津 300401)

1 引 言

模式匹配(又稱字符串匹配)是計算機科學的基礎問題之一[1],其在生物信息學[2]、網絡入侵[3]以及數據挖掘[4]等諸多領域具有重要應用,其問題實質是在一個相對長的序列串S上查找相對較短的模式P所有出現的位置及個數,這里序列串S以及模式P必須采用相同的字母表中[5].

近年來,為了滿足實際需要,在傳統通配符的基礎上,研究者們致力于間隙約束的模式匹配,其模式串P可以表示為:P=p1[a1,b1]p2…[aj,bj]pj+1…[am-1,bm-1]pm(1≤j≤m),其中aj和bj分別為pj和pj+1之間通配符之間最小和最大通配符的個數.這種間隙約束的模式在序列模式挖掘中亦有深入探索,如Min等人[6]在間隙約束下探索了序列中字符作用不均等的模式挖掘方法;Wu等人[7]采用網樹結構對周期性一般間隙的序列模式挖掘問題進行研究;Ding等人[8]在無重疊條件下探索了間隙約束的序列模式挖掘問題;文獻[9]中提出免預設間隔約束的對比序列模式高效挖掘算法,解決了序列數據挖掘問題;Dong等人[10]在挖掘重復模式問題上,提出了e-RNSP算法;Tan等人[11]根據位置的頻繁模式檢測,提出具有弱通配符間隙算法.在上述研究中,均需采用間隙約束模式匹配技術實現模式支持度的計算,從而判定模式的頻繁性,因此這種間隙約束模式匹配是序列模式挖掘的核心與基礎[12].

當前研究主要在非負間隙下進行研究,即aj≥0,但一般間隙約束能得到更多有價值的信息.因此,Fredriksson和Grabowski[13]對一般間隙約束下模式匹配問題進行研究,柴等[14]也提出一般間隙及一次性條件的模式匹配問題,并提出了具有完備性的有效算法DNCP;文獻[15]中提到一種基于靈活間隙模式匹配的服務器,并提出RNAPattMatch算法.與其他多種間隙約束模式匹配相比,無重疊模式匹配具有獨特的計算性能[16],在序列模式挖掘中更容易發現有價值的頻繁模式[17].

然而精確匹配有時不能完全提供給我們更多有價值的信息.在實際應用中大多數情況屬于近似模式匹配[18].近似模式匹配問題是指在模式串中其中一個字符與序列串中某個字符進行匹配時,兩個字符位置相同,但字符可以不完全相同.按照距離計算公式,可將近似條件約束模式匹配分為Hamming距離[19]、δ距離[20,21]及編輯距離[22]等.Hu等人[23]運用字符串相似搜索的通用Gram過濾器,提出了GFilter算法.

有鑒于此,本文提出了一般間隙近似無重疊模式匹配問題(Nonoverlapping Approximation Pattern Matching with General Gaps,簡稱NOAPMGG).該問題具有以下4個特點:它是一種嚴格的近似模式匹配問題,在匹配出現的過程中,不要求模式串上的字符與模式串上的字符完全相同;滿足了無重疊條件約束,允許同一字符在多個位置使用多次,但不允許在同一位置上同一字符使用多次;在模式串的設定中,允許間隙設置上出現多個負間隙.

2 NOAPMGG問題

定義1(序列串).序列串記為S,S=s1…si…sn(1≤i≤n),由字符組成的長度為n的信息源,其中,si∈∑,∑表示為一種符號的集合.在不同應用當中,∑的意義也會隨之不同.例如,在DNA基因序列中,∑由{A,C,G,T}組合而成.例如,S=s1s2s3s4s5=agcag,序列串S的長度為5.

定義2(模式串).模式串記為P,P=p1[a1,b1]p2…pj[aj,bj]pj+1…[am-1,bm-1]pm(1≤j≤m),其中,m指模式串P的長度,pj∈∑;aj、bj分別指間隙約束的下界與上界且取值必須是整數,并滿足aj≤bj.若aj≥0且bj≥0,則稱此區間為非負間隙條件約束;若aj或bj至少有一個小于0,則稱此區間為一般間隙條件約束.

定義3(出現).出現記為L,如果一個位置索引L=,其中,0≤lj≤n,0≤j≤m,lj≠lj-1,且滿足si=pj且:

(1)

則稱L為模式串P在序列串S上的一個出現.例如,模式串P=p1[a1,b1]p2[a2,b2]p3= A[0,3]T[0,2]C,則說明模式串P的長度為3;在字符A和字符T間根據間隙約束可知,兩個字符間可通配0到3個字符;而字符T和字符C間根據間隙約束可知,兩個字符間可通配0到2個字符.

定義5(Hamming距離).Hamming距離記為T,即表示兩個(相同長度)字對應位不同的數量.例如:序列串S1=s1s2s3s4s5=abcab與序列串S2=s1s2s3s4s5=acbac,在兩個序列串中,第2、3、5位置字符不同,因此,Hamming距離為3.

3 網樹與NOPARLA算法

3.1 網樹的概念

網樹的數據結構是一種樹型結構的拓展,目前被應用于解決多種模式匹配問題[24]、序列模式挖掘問題[25]和圖問題[26],其定義如下:

定義6(網樹).網樹是一種非線性數據結構,是由根結點與子網樹組成,所有相連的路徑之間用來表示雙親與孩子關系.

性質1.網樹具有樹結構的一些特點,如雙親、孩子、葉子結點、根結點、層等.此外,其具有如下獨特特點:

1)一棵網樹中可以存在一個或多個根結點;

2)網樹中非樹根結點也會存在多個雙親結點;

3)任何一個結點到它的祖先結點路徑可能不止一條;

4)給定序列中的結點可以在網樹的不同層上出現;

5)一個結點的孩子結點一定在網樹的同一層上出現.

定義7(樹根葉子路徑).在網樹中,能從樹根結點到葉子結點的一條完整的路徑.

定義9(最右雙親策略).從最后一層的最大結點開始查找,選擇最右雙親結點向樹根方向尋找樹根的過程.首先判斷最右雙親結點是否滿足間隙約束條件,若滿足則進行匹配,若不滿足則判斷次右雙親結點,直至找到相匹配的雙親結點.

定義10(最左孩子策略).從第一層的最小樹根結點出發,選擇最左孩子結點向葉子結點方向尋找路徑的過程.首先判斷最左孩子結點是否滿足間隙約束,若滿足則進行匹配,若不滿足則判斷次左孩子結點,直至找到與之相匹配的葉子結點.

3.2 NOPARLA算法的設計及分析

3.2.1 創建一般間隙近似網樹

在文獻[16]中,提出利用網樹的結構特點可以有效地解決無重疊條件約束下的模式匹配問題.由于本文研究的是近似距離下的模式匹配問題,因此需要在文獻[16]的創建網樹問題上進行一些修改.在修改算法的過程中會存在兩個難點:在處理結點關系時,會出現Hamming距離的實時更新;在創建父子關系前,需要進行預判斷處理.

在創建一般間隙近似網樹前,我們需提到一些與近似距離相關的定義如下.

1)計算該葉子結點到根結點之間的Hamming距離d時樹根到葉子的路徑數;

(2)

1)計算該根結點到葉子結點之間的Hamming距離d的路徑數.

(3)

算法1.CreateNetTree

輸入:長度為m的模式串P,長度為n的序列串S,Hamming距離閾值T

輸出:網樹NetTree

1.forj=1 tomstep 1do

2.fori=1 tonstep 1do

3. 依據定義11創建并計算結點的NHDRL值;

4.endfor

5.endfor

6.returnNetTree;

例1.給定序列串S=tgttagt,模式串P= t[-1,1]t[-1,2]a[-1,1]t,Hamming距離為1,根據CreateNetTree算法可以建立一棵近似網樹,如圖1所示.

該網樹的創建步驟如下:

1)當j=1時,將序列串中所有字符一次掃描插入第一層nettree[1],根據定義13,對所有結點進行初始值定義;

圖1 模式串P在序列串S對應的網樹Fig.1 Nettree for pattern P in sequence S

3.2.2 最左孩子策略

為解決近似距離閾值問題,在創建一般間隙近似網樹為每個結點都進行了在不同Hamming距離下路徑數量的定義與計算.若滿足近似度閾值條件,則根據樹根葉子路徑數匹配該結點的孩子結點.該方法可以有效提高匹配效率,具體過程如下:從最小根結點出發,根據間隙約束條件與Hamming距離的設定,當上一層結點與下一層結點建立關系時,根據公式(2)選擇與之匹配的最小孩子結點.其算法如算法2.

算法2.RootToLeaf

輸入:NetTree,模式串P長度為m,序列串S長度為n,閾值T,當前結點F

輸出:occin,各層結點F[m]集合

1.occin.position[1]=t;

2.F[1]=nettree[1][t];

3.ifF[1].used=true orF[1].toleave=falsethen

4.ifp[1].start≠s[F[1]]thenapprox=1;

5.forj=1 tomstep 1do

6.fort=1 tonettree[j][i].children.size()step 1do

7.maxh=T-approx;

8.endfor

9.ifnettree[j][child].toleave=false thenapprox=

approx+1;

10.endfor

11.endif

12.returnoccinandF[m];

3.2.3 近似閾值監測機制

根據文獻[17]可知,基于從最小根結點向下采用最左孩子結點找到較多的樹根葉子路徑,若模式串P的間隙約束略大或間隙條件過多時,結點的孩子結點也會成倍增加,由于近似度閾值的約束會造成該結點選擇的孩子結點并不是最小孩子結點,因此僅依靠最左孩子策略并不能找到更多的出現,由此可見需要引入近似閾值檢測機制.具體過程如下:從網樹的第二層開始,若第一次檢測到si=pj,且該結點在Hamming距離0～T-1時有路徑,則重新向上根據最左雙親策略尋找與之匹配的結點.其算法如算法3.

算法3.RootToLeaf_Approx

輸入:NetTree,模式串P(長度:m),序列串S(長度:n),閾值T,當前結點Q

輸出:occin,各層結點Q[t]集合

1.forj=1 tomstep 1do

2.if(j≤m)then

3. 采用算法2,從結點Q[t]出發,采用最左孩子策略,獲得子路徑Q及對應的Hamming距離為approx的子模式;

4.forx=jdownto 2 step-1do

5.fory=1 toNetTree[x][y].parents.size()step 1do

6.maxh=T-approx;

7.endfor

8.ifNetTree[x][child].toleave=falsethenapprox=approx+1;

9.endfor

10.endif

11.endfor

12.returnoccinandQ[t];

3.2.4 剪枝策略

定理1.若同一網樹中有互不相交的兩條路徑分別為A、B,則A、B經過的路徑上所有結點可構成兩個無重疊出現.

由定理1可知,同一棵網樹中有兩條互不相交的兩條路徑,可構成兩個無重疊出現.因此,根據最左孩子策略,在網樹上剪掉一條樹根葉子路徑后,并不會影響尋找其他的無重疊出現.但在進行剪枝后會出現一些不能到達葉子結點的結點,只需從葉子到樹根逐層檢測,將不可抵達葉子結點的結點進行剪枝,這樣可避免回溯過程.其算法如算法4.

算法4.Pruning

輸入:NetTree,模式串P的長度m,occin

輸出:NetTree

1.fori=mdownto 1 step-1do

2.position=occin.position[i];

3.ifi

4.pos=occin.position[i+1];

5.position_b=NetTree[m+1][pos].parent[0];

6.ifposition_b

7.endif

8.forposition=1 toNetTree[i].size()step 1do

9.current=NetTree[i][position];

10.ifcurrent.used=false ¤t.child=false

thenbreak;

11. updatecurrent.used;

12.ifcurrent.used=truethenupdate NHDRL;

13.endif

14.endfor

15.endfor

16.returnNetTree;

3.2.5 求解算法NOPARLA算法

根據本文提出的問題,NOPARLA算法將作為求解算法.該算法原理:

1)根據算法1將問題轉化成一棵網樹并根據性質11和公式(2)對結點進行NHDRL值定義;

2)從下往上根據性質12和公式(3)對結點進行NHDLR值定義;

3)采用算法3進行最左孩子策略和近似監測機制;

4)當找到一個無重疊出現時,采用算法4將網樹進行剪枝;

5)重復第3、4步,直至找到所有無重疊出現.

其算法如算法5:

算法5.NOPARLA

輸入:長度為m的模式串P,長度為n的序列串S,閾值T

輸出:符合條件的無重疊出現集合C

1. 采用算法1創建一棵近似網樹NetTree;

2.forj=m-1 downto 1 step-1do

3.fori=1 tonstep 1do

4. update NHDLR;

5.endfor

6.endfor

7.fort=1 toNetTree[1].size()step 1do

8. 采用算法3 RootToLeaf_Approx(N,Q[t]),從結點Q出發,采用最左孩子策略;

9. 采用算法4 Pruning更新NetTree;

10.C=Q∪occin;

11.endfor

12.returnC;

例2.給定序列串S=tgttagt和模式串P= t[-1,1]t[-1,1]a[-1,2]t,Hamming距離為1.

計算過程如下:

1)模式串P在序列串S中所建立的對應網樹,如圖1所示.

2)根據公式(3)對網樹重新定義NHDLR值.

根據NOPARLA算法,可以找到三個無重疊出現,即<1,3,2,1>、<3,4,3,4>、<4,6,5,7>.

3.3 算法分析

由于網樹最多只有m層,每一層結點數最多只有n個結點,由此可見,每一個孩子結點最多會有W個雙親結點,即W=max(bj-aj+1),且aj為間隙約束的下界,bj為間隙約束的上界,近似閾值為T.因此,NOPARLA算法的空間復雜度為O(m*n*W*T).其中,m為模式串的長度,n為序列串的長度,W為模式串的最大間隙長度W=max(bj-aj+1),T為近似閾值.

在分析本文算法時間復雜度之前,先分析CreateNetTree算法、RootToLeaf算法、RootToLeaf_Approx算法以及Pruning算法的時間復雜度.易知CreateNetTree算法創建一棵一般間隙近似網樹的時間復雜度為O(m*n*W*T);在RootToLeaf算法中,每個結點最多有W個孩子結點,網樹共有m層,近似閾值為T,所以RootToLeaf算法的時間復雜度為O(m*W*T);RootToLeaf_Approx算法的時間復雜度為O(m2*W*T);根據Pruning算法可知,該算法的時間復雜度為O(m3*W*T).由于NOPARLA算法最多運行n-m次剪枝算法,由此本文算法的時間復雜度為O(m3*n*W*T).

4 實驗結果及分析

4.1 實驗環境及實驗數據

本文中實驗運行的軟、硬件環境為:Inter(R)Core(TM)i5-3210處理器、2.50GHz主頻、4.00GB內存、Windows8.1 專業版操作系統的計算機.完成實驗工具為Microsoft Visual C++6.0.本文所使用的測試序列串均為真實的生物數據,DNA序列均來自美國國家生物計算信息中心網站1.為體現模式串的間隙約束具有一般性規律,在本文中設計了九個具有一般性的模式串. 表1和表2分別給出了實驗中使用的序列串和模式串的特征.

表1 真實生物序列
Table 1 Sequences of real biological data

序號片段名稱位點片段長度S1Segment 1CY0585632286S2Segment 2CY0585622299S3Segment 3CY0585612169S4Segment 4CY0585561720S5Segment 5CY0585591516S6Segment 6CY0585581418S7Segment 7CY058557982S8Segment 8CY058560844

4.2 實驗結果

為了測試本文NOPARLA算法的性能,我們又提出了三種對比算法,即NOPALR算法、NOPARL算法和NOPALRA算法.這三種對比算法均采用網樹結構進行求解,其中NOPALR算法采用最右雙親策略,當找到一個無重疊出現時將網樹進行剪枝;NOPARL算法采用最左孩子策略,當找到一個無重疊出現時將網樹進行剪枝;NOPALRA算法采用最右雙親策略,再采用近似閾值檢測,當找到一個無重疊出現時將網樹進行剪枝.

本文選取表2中P1～P9中的9個模式串與表1中S1～S8中的8個序列串作為本文的對比實驗數據,在選取近似閾值T=1與T=2進行對比實驗,并對這四種算法的求解質量與求解速度進行對比.表3和表4分別給出了近似閾值T=1與T=2情況下,模式串P1～P9在序列串S1～S8上的出現總數.

表2 模式串
Table 2 Patterns

表3T=1時,序列S1～S8在模式P1～P9上的出現總數
Table 3 Total results ofS1～S8 inP1～P9 whenT=1

出現數P1P2P3P4P5P6P7P8P9NOPALR341045303193255514581725123027493199NOPARL353244143452268818821905161629303218NOPARLA351843573382273016791882140030943450NOPARLA360345563445265819722016169331473508

表4T=2時,序列S1～S8在模式P1～P9上的出現總數
Table 4 Total results ofS1～S8 inP1～P9 whenT=2

出現數P1P2P3P4P5P6P7P8P9NOPALR6468100866216499126192500211843604469NOPARL6575101236389533826012524227143604445NOPARLA6438101126220520427662713226846684828NOPARLA6664101976442508727882753240548044971

為了能更清晰的反映出四種算法分別在近似閾值T=1與T=2情況下,模式串P1～P9在序列串S1～S8上的時間性能,四種算法的平均運行時間見表5和表6.

4.3 實驗結果分析

1)在求解性能方面,NOPARLA算法整體比NOPARL算法、NOPALR算法與NOPALRA算法較好.

通過表3～表4,即T=1和T=2時,序列串S1～S8在模式P1～P9上的出現總數.從兩張表中,清晰看出NOPARLA算法從總體來說,求解性能優于其他三種算法.在模式串P3上,NOPARL算法求解性能優于NOPARLA算法;在模式串P4上,NOPALRA算法與NOPARLA算法均沒有NOPARL算法的求解性能好;產生這樣的原因一是由于序列串的長短會造成這樣的原因,由于在NOPARL算法與NOPALR算法并沒有判斷近似監測機制,因此,對網樹造成較小的影響;二是因為當序列串足夠長的時候,為了能找到能多有價值的出現,算法對近似距離做檢測,當滿足一定條件時,會觸發近似監測機制,使得產生的新無重疊出現比原有無重疊出現對網樹的影響減小,因而提高了算法的求解性能.

表5T=1時,序列S1～S8在模式P1～P9上的平均運行時間
Table 5 Average running time ofS1～S8 inP1～P9 whenT=1

平均時間P1P2P3P4P5P6P7P8P9NOPALR4392844394141164171714579301424NOPARL52235052548814332145176711331660NOPARLA50634052548613522029168410781625NOPARLA53937354351014812271183411491735

表6T=2時,序列S1～S8在模式P1～P9上的平均運行時間
Table 6 Average running time ofS1～S8 inP1～P9 whenT=2

平均時間P1P2P3P4P5P6P7P8P9NOPALR74041277955128116184369120244986NOPARL83846790877731316840409422405250NOPARLA77143880361130686908410222425521NOPARLA86449493477232747387429523855855

2)在求解性能方面,通常情況下NOPARLA算法比NOPALRA算法可以發現更多的出現.

從表3～表4可以看出,在模式串P1～P9上,當T=1時,NOPARLA算法有七次獲得了最多出現數,當T=2時,NOPARLA算法有八次獲得了最多出現數.

3)在時間運行方面,改進后的NOPARLA算法與NOPALRA算法運行時間總體比未改進的NOPARL算法與NOPALR算法較長.

從表5～表6可以看出,即T=1和T=2時,序列串S1～S8在模式P1～P9上的平均運行時間下,NOPARLA算法在運行時間上總體來說消耗最多,NOPALR算法消耗的時間最少,并且可以看出NOPARLA算法與NOPALRA算法消耗的時間,從整體來說均比NOPARL算法與NOPALR算法時間長,造成這種現象別的主要原因是在處理無重疊出現時,對近似距離做檢測,通過檢測可以找到更多有價值的出現,因此會增加時間的運行.

4)隨著近似閾值T的增加,找到的出現數與運行時間也會隨之增加.

通過表3～表4,即T=1和T=2時,序列串S1～S8在模式P1～P9上的出現總數.換言之,當T從0到2逐漸增大時,出現數也會隨之增加,產生這樣現象的原因是由于近似閾值T可轉換為近似距離為T-approx及近似距離為approx問題.從表5～表6可以看出,即T=1、T=2時,序列串S1～S8在模式串P1～P9上的平均運行時間也會隨之增加,這與上一章中提高的算法時間復雜度與空間復雜度分析一致.

通過本章做的大量實驗進行分析可知,與其他3種算法相比,NOPARLA算法從整體來說性能最好.

5 結 論

本文提出了一般間隙近似無重疊模式匹配問題,即NOAPMGG問題.該問題需要滿足無重疊模式的條件約束,為增加模式串的一般性而提出近似模式匹配,為進一步增加模式串的靈活性而提出一般間隙模式匹配問題,來解決NOAPMGG問題.本文提出NOPARLA算法,該算法根據網樹的結構特點,采用最左孩子策略,使用近似閾值機制找到一個無重疊出現.該算法的空間復雜度與時間復雜度分別為O(m*n*W*T)和O(m3*n*W*T),其中m為模式串的長度,n為序列串的長度,W為模式串的最大間隙長度,T為近似閾值.通過大量實驗證明了,NOPARLA算法具有較高的可行性與高效性.