張文影，昌 攀，鐘 誠

(廣西大學 計算機與電子信息學院，廣西 南寧 530004)

0 引 言

目前主流的DNA序列分類方法[1,2]有基于序列比對的分類方法[3-5]和無需比對的分類方法[6-8]。盡管基于比對的分類方法可以獲得較高的分類精度，但仍存在對功能相似結構差異較大的DNA序列分類時精度較低、比對大量序列增加了計算復雜度[9]等弊端。

為了克服基于比對的序列分類方法的局限性，學者們提出了無需比對的DNA序列分類方法。常用的無需比對的分類使用了K元組[10]方法，文獻[11]通過實驗表明K元組較適用于核小體DNA序列的分類，在對核小體的分類中取得了較高的精度和敏感度。然而，基于K元組方法的計算結果很大程度上取決于參數(shù)K的設定，此外，當K值較大時，這些方法也面臨著特征向量維數(shù)較高和計算復雜等問題。另一些DNA序列分類研究則使用有效的壓縮算法對序列進行壓縮，進而分類DNA序列。文獻[12]介紹了3種DNA壓縮算法，并通過實驗闡明DNA壓縮算法在序列分類領域的應用。文獻[13]利用Voss映射將DNA序列映射成二進制指示序列，然后利用離散傅里葉變換(discrete fourier transformation，DFT)處理該指示序列，進而分析DNA序列之間的進化關系。這種方法只考慮了單個堿基，而忽略了DNA序列堿基化學性質、出現(xiàn)頻率和位置等重要生物學信息。本文結合K元組方法和離散傅里葉變換DFT提取序列堿基化學性質、出現(xiàn)頻率和位置的特征值，將任意長度的DNA序列映射成等長的特征向量，從而更準確地分類DNA序列。

1 分類算法

1.1 DFT與功率譜

DFT變換可以將數(shù)值序列從時域序列轉化為頻域序列，其中序列的信息也就相應地轉化為不同頻率下的幅度和相位，這樣就可以挖掘出序列中隱藏的重要信息。DFT已被應用于DNA序列分析研究中[14]。

利用DFT變換，長度為N的時間序列f(n)在特定的k頻率上的映射可表示為F(k)

(1)

F(k)在頻率k上的功率譜PS(k)[15]

(2)

當k=0時，頻率為0時的特定功率譜PS(0)為

(3)

1.2 提取DNA序列的特征向量

DNA序列是由4種堿基：腺嘌呤A、鳥嘌呤G、胞嘧啶C和胸腺嘧啶T組成的一維線性序列。任一條DNA序列可以被唯一描述為3種獨立的嘌呤(R)和嘧啶(Y)的分布、氨基(M)和羰基(K)的分布、強氫鍵(S)和弱氫鍵(W)的分布[16]。依據(jù)4種堿基的不同化學性質可以將堿基分為以下3類：

(1)嘌呤R=A、G，嘧啶Y=C、T；

(2)氨基M=A、C，羰基K=G、T；

(3)強氫鍵S=G、C，弱氫鍵W=A、T。

依據(jù)這3類堿基分類，可將任一條DNA序列映射成3條特定堿基組成的符號序列[17]。例如，DNA序列x=GTCCACTGGCATGGT可映射為3條符號序列φRY(x)=RYYYRYYRRYRYRRY，φMK(x)=KKMMMMKKKMMKKKK，φWS(x)=SWSSWSWSSSWWSSW。

參照文獻[17]的方法，將φβ(x)(β∈{RY,MK,WS})轉換成如下12條二進制指示序列H(x)={hRR(x),hRY(x),hYY(x),hYR(x),hMM(x),hMK(x),hKK(x),hKM(x),hWW(x),hWS(x),hSS(x),hSW(x)}

(4)

其中，α∈{RR,RY,YY,YR,MM,MK,KK,KM,WW,WS,SS,SW}，φβ(n)為序列φβ(x)的第n個符號，N為DNA序列的長度。

例如，對于序列φRY(x)、φMK(x)和φWS(x)可生成如下12條二進制指示序列：hRR(x)=00000001000010，hRY(x)=10001000101001，hYY(x)=01100100000000，hYR(x)=00010010010100，hMM(x)=00111000010000，hMK(x)=00000100001000，hKK(x)=10000011000111，hKM(x)=01000000100000，hWW(x)=00000000001000，hWS(x)=01001010000100，hSS(x)=00100001100010，hSW(x)=10010100010001。

設Sα為hα(x)中子序列α的數(shù)量，則可計算子序列α在二進制指示序列hα(x)中的出現(xiàn)頻率T(α)

T(α)=Sα/(N-1)

(5)

依據(jù)式(1)，將二進制指示序列hα(x)進行DFT變換后得到序列Fα(k)

(6)

其中，L=N-1為二進制指示序列hα(x)的長度。根據(jù)式(2)，F(xiàn)α(k)在特定頻率k上的DFT功率譜PSα(k)為

(7)

將帕賽瓦爾定理運用在DFT變換中得到

(8)

(9)

(10)

其中，E(k)是每個功率譜PSα(k)的加權值

(11)

于是，可提取得到的二進制指示序列hα(x)的離散傅里葉功率譜的特征值Mα。

式(10)充分體現(xiàn)了二進制指示序列hα(x)子序列種類、出現(xiàn)頻率和位置的序列特征信息。這樣，就可以將DNA序列x轉換為一個12維的特征向量(MRR,MRY,…,MSW)。

1.3 分類算法描述

為了完成DNA序列的分類，將提取的DNA序列特征值利用式(12)分別計算出待分類DNA序列對應特征向量與每個樣本DNA序列特征向量的距離，給定近鄰數(shù)目K，通過K-NN分類器[18]，選取K個距離最小的樣本DNA序列組成K-最近鄰樣本集合S，統(tǒng)計S中每種DNA序列類別出現(xiàn)的次數(shù)，然后選擇出現(xiàn)頻率最大的DNA序列的類別作為待分類DNA序列的類別。

(12)

算法流程如圖1所示。

圖1 IAF-DNASC算法流程

算法1形式描述了改進的無需比對的分類算法(簡記為IAF-DNASC算法)。

算法1：IAF-DNASC

Input：待分類DNA序列集X={x1,x2,…,xm}，樣本DNA序列集Y={y1,y2,…,yn}，近鄰數(shù)目K

Output：輸出DNA序列xi的類標簽

Begin

(1)按式(4)將DNA序列xi映射為12條二進制指示序列H(xi)，i=1,2,…,m；

(2)按式(4)將DNA序列yj映射為12條二進制指示序列H(yj)，j=1,2,…,n；

(5)按式(10)對H(xi)的功率譜xi_PSα(k)進行特征提取，得到一個12維特征向量，記為xi_Mα，i=1,2,…,m；

(6)按式(10)對H(yj)的功率譜yj_PSα(k)進行特征提取，得到一個12維特征向量，記為yj_Mα，j=1,2,…,n；

(7)對每個待分類序列xi，i=1,2,…,m，執(zhí)行以下步驟：

1)用式(12)計算xi與yj之間的距離：

3)輸出DNA序列xi的類標簽Label(xi)

其中，c表示所有可能的類標簽。

End

IAF-DNASC算法依據(jù)DNA序列不同堿基化學性質，結合K元組方法將DNA序列映射成12條二進制指示序列，二進制指示序列中保存了DNA序列的堿基化學性質，能充分表征DNA序列；對映射生成的二進制指示序列進行DFT變換得到離散傅里葉功率譜，綜合利用二進制指示序列的不同子序列出現(xiàn)頻率和位置信息，對功率譜采取匹配加權值的方法提取特征值，解決不等長離散傅里葉功率譜的比較問題，使得可采用歐氏距離準確計算出不同類別DNA序列各特征向量之間的距離，進而利用計算得到的距離值通過K-NN分類器準確分類DNA序列。

2 實 驗

2.1 實驗數(shù)據(jù)和環(huán)境

實驗選取了3個核小體DNA序列的數(shù)據(jù)庫：Caenorhabditis elegans(CE)、Homo sapiens(HM)和Drosophila melanogaster(DM)。這些數(shù)據(jù)可以在網(wǎng)站http://lin.uestc.edu.cn/server/iNucPseKNC/dataset下載。HM數(shù)據(jù)庫包含2273條核小體形成DNA序列和2300條核小體抑制DNA序列；CE數(shù)據(jù)庫包含2567條核小體形成DNA序列和2608條核小體抑制DNA序列；DM數(shù)據(jù)庫包含2900條核小體形成DNA序列和2850條核小體抑制DNA序列。核小體是染色體的基本組成單位，也是染色體中的最重要的重復單元。有關這3個核小體DNA序列的實驗數(shù)據(jù)提取和篩選的具體方法可參見文獻[19]。每個核小體DNA序列數(shù)據(jù)庫中包含了兩類核小體DNA序列：一類是核小體形成DNA序列樣本(positive data)，另一類是核小體抑制DNA序列樣本(negative data)。實驗利用精度(A)，敏感度(Se)和特異性(Sp)3個指標評價DNA序列分類性能

其中，TP表示核小體形成DNA序列正確分類的序列數(shù)，F(xiàn)P為核小體形成DNA序列錯誤分類的序列數(shù)，TN表示核小體抑制DNA序列正確分類的序列數(shù)，F(xiàn)N為核小體抑制DNA序列錯誤分類的序列數(shù)。

實驗使用的計算機為8核Intel(R) Core(TM) i7-6700K CPU @ 4.00 GHz 4.00 GHz處理器、內(nèi)存容量為32.0 GB。運行操作系統(tǒng)為Windows 10。軟件運行環(huán)境為Visual studio 2012，采用C++編程實現(xiàn)算法。實驗將本文的IAF-DNASC算法與已有的同類代表性算法correlation，UCD-gzip和UCD-deflate的分類性能進行比對分析。

2.2 實驗結果

在實驗過程中，我們利用十折交叉驗證法，通過隨機抽樣將每個數(shù)據(jù)集劃分為10個子集，每次選擇其中的9個子集作為訓練數(shù)據(jù)，剩余的1個子集作為測試數(shù)據(jù)，重復10次實驗，取實驗結果的平均值。K-NN分類算法中參數(shù)K的最優(yōu)取值取決于測試數(shù)據(jù)集，參照文獻[11]，實驗選取K=1,3,5,7,…，21，L=2。

圖2、圖3和圖4分別給出了4種算法在3個數(shù)據(jù)集CE、DM和HM上取不同K值時獲得的DNA序列分類精度、敏感度和特異性。

圖2 4種算法在CE數(shù)據(jù)集上的分類精度、敏感度和特異性

圖2的結果表明：在數(shù)據(jù)集CE上，對于不同的K值，IAF-DNASC的平均分類精度和敏感度均高于算法UCD-deflate、UCD-gzip和correlation的平均分類精度和敏感度，尤其是遠高于UCD-deflate和UCD-gzip的分類精度和敏感度；IAF-DNASC的特異性略低于UCD-gzip和UCD-deflate，但高于correlation的值。

圖3 4種算法在DM數(shù)據(jù)集上的分類精度、敏感度和特異性

由圖3可看出：在數(shù)據(jù)集DM上，對于不同的K值，IAF-DNASC算法的分類精度和特異性均高于算法UCD-deflate、UCD-gzip和correlation，IAF-DNASC算法的分類敏感度與UCD-gzip幾乎一樣，高于UCD-deflate和correlation。

圖4 4種算法在HM數(shù)據(jù)集上的分類精度、敏感度和特異性

由圖4可知，在數(shù)據(jù)集HM上，IAF-DNASC算法的分類精度和特異性在4種算法中都是最高的，IAF-DNASC算法分類時取得敏感度略低于correlation，但高于UCD-gzip和UCD-deflate。

上述這些結果表明，在數(shù)據(jù)集CE、DM和HM上，整體上，IAF-DNASC算法的分類性能優(yōu)于correlation、UCD-gzip和UCD-deflate，這是因為IAF-DNASC算法綜合利用堿基化學性質、出現(xiàn)頻率和位置信息，采用匹配加權值方法提取DNA序列DFT功率譜特征值，能夠充分挖掘和表征DNA序列的重要特征，利用這些特征信息可以更準確地計算出不等長DNA序列間的距離，從而可以更為準確地完成DNA序列的分類。

3 結束語

本文結合K元組方法和離散傅里葉變換，提出了改進的無需比對的DNA序列分類算法，該算法有效地提取了DNA序列堿基化學性質、出現(xiàn)頻率和位置特征信息，通過這些提取的特征信息可以更加準確地計算出不等長DNA序列之間的距離，利用計算得到的距離值通過K-NN分類器實現(xiàn)更準確地分類DNA序列，從而整體上獲得更高的分類精度、敏感度和特異性。RNA序列、蛋白質序列舉有不同的生化性質，下一步工作將研究將本文方法拓廣到RNA序列、蛋白質序列的序列挖掘中。