昌 攀,鐘 誠

(廣西大學(xué) 計算機(jī)與電子信息學(xué)院 廣西高校并行分布式計算技術(shù)重點實驗室,南寧 530004)

1 引 言

隨著后基因時代(Post-genome era)的到來,越來越多的DNA序列信息被挖掘出來.使用傳統(tǒng)的生物測序?qū)嶒灧治鯠NA序列已經(jīng)變得越來越復(fù)雜和低效.聚類分析技術(shù)為分析快速增長的DNA序列提供了切實可行的方法,使預(yù)測同類別基因的功能變得方便[1].屬于同一類別的DNA序列彼此間的差異性較小,生物學(xué)上的功能彼此相似,有利于分析物種間的進(jìn)化關(guān)系[2,3].

物種聚類需要計算DNA序列之間的相似度.目前有兩種常用的測量序列之間相似度的方法.一種是基于序列比對的DNA序列相似度計算方法,它通過獲得序列間的全局比對信息來評判DNA序列間的相似度,常用算法有BLAST[4]、FASTA[5]、UCLUST[6]和CD-HIT[7].盡管序列比對可以獲得較好的結(jié)果,但當(dāng)處理海量DNA序列時,不僅消耗大量的內(nèi)存和計算時間,而且很難獲得一個讓人滿意的比對結(jié)果[8].另一種是無需比對的DNA序列相似度計算方法,其基本思想是將任一條DNA序列轉(zhuǎn)化成一個特定長度的特征向量,以快速計算出DNA序列之間的相似度.

已有的一些無需比對的DNA序列相似度計算方法結(jié)合統(tǒng)計學(xué)模型生成特征向量,典型的馬爾可夫模型(Markov Model)[9,10]被廣泛應(yīng)用于生物信息學(xué)的研究中.然而,馬爾可夫模型存在著一些不足.Lu等人在文獻(xiàn)[11]中指出:無需比對的馬爾可夫模型的假設(shè)削弱了模型的適用性.K元組方法[12,13]是一種常用的無需比對的DNA序列相似度計算方法,它通過一個長度為K的滑動窗口對DNA序列進(jìn)行分割,統(tǒng)計每個元組出現(xiàn)的頻數(shù),構(gòu)造一個基于元組頻數(shù)值的4K維特征向量,然后利用特征向量計算DNA序列之間的相似度.運(yùn)用K元組方法可以將任一條DNA序列轉(zhuǎn)化成一個定長的特征向量,有助于解決不同DNA序列長度不等長的問題.K元組方法的不足是:不能完全描述DNA序列的全部特征,對于子序列特征的提取也未能考慮到子序列的位置特征信息.文獻(xiàn)[14]提出的DMK算法融入K元組的位置特征,使用“信息熵”重新計算K元組的特征值,使計算得到的DNA序列相似度優(yōu)于K元組算法計算得到的DNA序列相似度.文獻(xiàn)[15]的CPF算法采用“信息熵”的方法,按照生物信息學(xué)的知識將DNA序列堿基進(jìn)行分類,通過發(fā)掘DNA序列的局部特征計算DNA序列的相似度對物種聚類.相對于K元組方法,CPF算法提高了聚類結(jié)果的準(zhǔn)確度.近年來,文獻(xiàn)[16,17] 將DNA序列的信息數(shù)值化,利用離散傅里葉變換(DFT)思想來處理DNA數(shù)值化序列,使測量DNA序列的突變率的效率更高,為度量物種的親緣性關(guān)系提供了一種新思路.文獻(xiàn)[16,17]的方法只考慮了組成DNA序列的4種堿基{A,T,C,G}各自的二進(jìn)制指示序列,忽略了DNA序列子序列所固有的種類、含量和位置等生物特征,生成的系統(tǒng)進(jìn)化樹未能準(zhǔn)確呈現(xiàn)物種的進(jìn)化關(guān)系.

本文通過離散傅里葉變換提取DNA序列中子序列種類、含量和位置3種重要生物特征,給出一種改進(jìn)的應(yīng)用于物種聚類的無需比對的DNA序列相似度計算算法AFCS_DFT.AFCS_DFT算法能計算出更準(zhǔn)確的DNA序列相似度和更準(zhǔn)確聚類物種.

2 DNA序列相似度計算

2.1 傅里葉變換與功率譜

在信號處理過程中,離散傅里葉變換(Discrete Fourier Transformation,DFT)將時間域上的序列轉(zhuǎn)換成頻率域,使時間域上的一些重要信息顯現(xiàn)出來.DFT變換被廣泛應(yīng)用于DNA研究(例如基因功能預(yù)測、編碼蛋白質(zhì)區(qū)域檢測、層次聚類等)中[18].離散傅里葉變換的功率譜可以揭示DNA核苷酸的位置分布和頻率的狀況.

在信號處理過程中,對于一段長度為N的信號序列χ(n)在特定的k頻率下進(jìn)行離散傅里葉變換轉(zhuǎn)換成序列X(k)[19]:

(1)

PS(k)=|Χ(k)|2,k=0,1,…,N-1

(2)

根據(jù)式(1)和(2),當(dāng)k=0時,計算出頻率為0時的特定功率譜PS(0):

(3)

當(dāng)χ(n)為二進(jìn)制指示序列時,PS(0)的值等于二進(jìn)制指示序列中1 的個數(shù)總和的平方.

2.2 AFCS_DFT算法

本文給出的改進(jìn)的應(yīng)用于物種聚類的DNA序列相似度計算方法AFCS_DFT(Alignment- Free Computing Similarity using Discrete Fourier Transformation)的思想是:首先,按照組成DNA序列的堿基生物特性[15],將任一條DNA序列生成3條等長的特征序列,從特征序列中辨別出子序列的種類、含量和位置3種重要生物特征,將這3條等長的特征序列映射成相應(yīng)的12條二進(jìn)制指示序列,利用離散傅里葉變換和后面給出的式(14),從功率譜(能量)中提取出這12條二進(jìn)制指示序列的生物特征向量值,從而將一條DNA序列轉(zhuǎn)換成一個24維的特征向量,最后使用歐式距離計算得到更準(zhǔn)確的DNA序列相似度.

通過分析組成DNA序列的4種堿基:胞嘧啶(C)、鳥嘌呤(G)、腺嘌呤(A)、胸腺嘧啶(T)的化學(xué)特性,將這4種堿基分成如下3類[14]:

1)嘌呤集R={A,G},嘧啶集Y={C,T};

2)胺基集M={A,C},酮基集K={G,T};

3)弱氫鍵集W={A,T},強(qiáng)氫鍵集S={C,G}.

對于某一個特定的堿基集,任一條DNA序列都可以映射成特定堿基集合組成的序列S(n)[14].例如,對于序列S:ACAGACACCATGGTGCACCTCCT,將S映射到嘌呤集R和嘧啶集Y中得到S1序列:

S1:RYRRRYRYYRYRRYRYRYYYYYY

將S映射到胺基集M與酮基集K中,得到S2序列:

S2:MMMKMMMMMMKKKKKMMMMKMMK

將S映射到弱氫鍵集W與強(qiáng)氫鍵集S中,得S3序列:

S3:WSWSWSWSSWWSSWSSWSSWSSW.

對于映射到特定集合上的序列S(n),采用K元組(這里K=2)算法[12,13],參照文獻(xiàn)[14]的方法生成相應(yīng)子序列的二進(jìn)制指示序列fα(n):

(4)

其中子序列α對應(yīng)于嘌呤集R和嘧啶集Y的{RR,RY,YR,YY}序列集合、對應(yīng)于胺基集M與酮基集K的{MM,MK,KM,KK}序列集合以及對應(yīng)于弱氫鍵集W與強(qiáng)氫鍵集S的{WW,WS,SW,SS}序列集合,亦即子序列α∈{RR,RY,YR,YY,MM,MK,KM,KK,WW,WS,SW,SS}集合.例如,對于S1序列可以分解成4條獨(dú)立的二進(jìn)制序列:

RR型的fRR序列,fRR=0011000000010000000000;

RY型的fRY序列,fRY=1000101001001010100000;

YR型的fYR序列,fYR=0100010010100101000000;

YY型的fYY序列,fYY=0000000100000000011111.

于是,可以將任一條長度為N的DNA序列分解成長度為L=(N-1)的12條二進(jìn)制指示序列fα(n),n=1,…,(N-1),fα(n)∈{fRR,fRY,fYR,fYY,fMM,fMK,fKM,fKK,fWW,fWS,fSW,fSS}.

通過進(jìn)一步分析序列fα(n)的組成,我們可以發(fā)現(xiàn)fα(n)可以保留了DNA序列子序列的種類、含量和位置3種重要生物特征.然后,利用離散傅里葉變換對fα(n)進(jìn)行子序列特征值的提取,特征值中將保留DNA序列子序列的種類、含量和位置3種重要生物特征.

為了提取DNA序列子序列的含量,我們給出式(5)以計算出每種特定的子序列α在二進(jìn)制指示序列fα中的含量C(α)計算如下:

(5)

其中Nα為二進(jìn)制指示序列fα(n)中相應(yīng)子序列α的數(shù)量,即二進(jìn)制指示序列fα(n)中“1”的數(shù)量.

利用離散傅里葉變換對fα(n)中的12條二進(jìn)制指示序列{fRR,fRY,fYR,fYY,fMM,fMK,fKM,fKK,fWW,fWS,fSW,fSS}提取得到的特征值保留了DNA序列子序列的種類、含量和位置3種重要生物特征,這將比采用K元組算法[11,12]提取出更多的DNA序列信息,克服了文獻(xiàn)[15,16]提取單一生物特征的不足,解決了CPF算法[14]計算復(fù)雜的問題.

以下是運(yùn)用離散傅里葉變換提取12條二進(jìn)制指示序列

1www.ncbi.nlm.nih.gov

中DNA序列子序列的種類、含量和位置3種特征值的具體方法.

根據(jù)式(1),將二進(jìn)制指示fα(k)進(jìn)行DFT變換后得到序列Fα(k):

(6)

其中L=N-1為二進(jìn)制指示序列fα(n)的長度.根據(jù)式(2),Fα(k)在特定頻率k下的DFT功率譜PSα(k)表示為:

PSα(k)=|Fα(k)|2,k=0,1,…,L-1

(7)

(8)

根據(jù)廣義帕賽瓦爾定理(Parseval's theorem):一個信號所含有的能量(功率)恒等于此信號在完備正交函數(shù)集中各分量能量(功率)之和[19].將帕賽瓦爾定理運(yùn)用在離散傅里葉變換中得到[19]:

(9)

(10)

(11)

(12)

離散傅里葉變換(DFT)存在對稱性[19].于是,式(12)可以調(diào)整為:

(13)

(14)

為了計算各序列特征向量之間的相似度,需要考慮以下約束條件:對于給定DNA序列集合S,序列x,y∈S,x與y之間的相似度d(x,y)需要滿足以下3個性質(zhì)[20]:(1)非負(fù)性:d(x,y)≥0,d(x,y)=0,當(dāng)且僅當(dāng)x=y;(2)對稱性:d(x,y)=d(y,x)恒成立;(3)滿足三角不等式:對于z∈S,d(x,y)≤d(x,z)+d(z,y).

采用歐氏距離計算任意2條DNA序列x和y特征向量之間的距離:

(15)

計算得到的距離d(x,y)滿足上述3個性質(zhì).

3 實 驗

實驗使用的計算機(jī)為4核AMD Athlon (tm)ⅡX4 645 3.10GHz處理器、內(nèi)存容量4GB.運(yùn)行操作系統(tǒng)Windows 7.軟件運(yùn)行環(huán)境是Visual studio 2012.采用C++語言編程實現(xiàn)算法.

2組實驗數(shù)據(jù)分別是:10個物種的全β-球蛋白基因序列[21]和美國國立生物技術(shù)信息中心NCBI1提供的30種甲型流感病毒的神經(jīng)氨酸苷酶基因序列.實驗首先運(yùn)用AFCS_DFT算法將任一條DNA序列映射成相應(yīng)的特征向量,然后根據(jù)式(15)計算任意個物種DNA序列之間特征向量的相似度,最后將DNA序列之間的相似度矩陣導(dǎo)入應(yīng)用軟件MEGA 6[17]中,使用構(gòu)造UPGMA樹[22]的方法構(gòu)建系統(tǒng)進(jìn)化樹進(jìn)行DNA序列的聚類.實驗將本文給出的AFCS_DFT算法與已有的同類算法CPF[15]和DFT[17]構(gòu)建的系統(tǒng)進(jìn)化樹進(jìn)行比對分析.

3.1 10個物種的全β-球蛋白基因序列構(gòu)建的系統(tǒng)進(jìn)化樹

文獻(xiàn)[21]中給出的10個物種的全β-球蛋白基因序列常用作測試物種之間的進(jìn)化關(guān)系.運(yùn)行AFCS_DFT、DFT和CPF算法分別計算得到10個物種全β-球蛋白基因序列之間的相似度.依據(jù)得到的全β-球蛋白質(zhì)基因序列之間的相似度,利用軟件MEGA 6[17],使用構(gòu)造UPGMA樹[22]方法構(gòu)建10個物種的全β-球蛋白基因序列的系統(tǒng)進(jìn)化樹如圖1-圖3所示,圖中橫坐標(biāo)代表物種相似度的度量的標(biāo)尺.

從生物信息學(xué)的角度分析,物種 Human與Gorilla、Chimpanzee的親緣關(guān)系遠(yuǎn)大于Human與Goat、Bovine的親緣關(guān)系[21].從圖1-圖3的結(jié)果可知,AFCS_DFT與CPF構(gòu)建的系統(tǒng)進(jìn)化樹符合這一特征,而DFT構(gòu)建的系統(tǒng)進(jìn)化樹結(jié)果卻相反.AFCS_DFT與CPF可以將哺乳動物與非哺乳動物(Gallus是非哺乳動物,其他物種是哺乳動物)區(qū)分開來同時與Bovine和Goat物種相比,Human、Gorilla、Chimpanzee和Lemur物種的相似度與Mouse和Rat物種更接近,這些結(jié)果與文獻(xiàn)[21,23]的研究結(jié)果一致.而DFT構(gòu)建的系統(tǒng)進(jìn)化樹未能將卵生動物Gallus與胎生動物Human、Lemur、Gorilla、Chimpanzee、Goat、Bovine、Mouse、Rat物種區(qū)分.

圖1 AFCS_DFT算法構(gòu)建的系統(tǒng)進(jìn)化樹Fig.1 Phylogenetic tree constructed by AFCS_DFT

圖2 DFT算法構(gòu)建的系統(tǒng)進(jìn)化樹Fig.2 Phylogenetic tree constructed by DFT

圖3 CPF算法構(gòu)建的系統(tǒng)進(jìn)化樹Fig.3 Phylogenetic tree constructed by CPF

另一方面,與CPF相比,AFCS_DFT符合生物信息學(xué)上具有高度親緣關(guān)系的物種DNA序列之間差異較小的特征.此外,AFCS_DFT增大了親緣關(guān)系稀薄的物種DNA序列之間的差異.

上述結(jié)果表明:AFCS_DFT、CPF在10個物種的全β-球蛋白基因序列數(shù)據(jù)實驗中,發(fā)掘物種之間進(jìn)化關(guān)系的性能優(yōu)于DFT;另一方面,AFCS_DFT在顯現(xiàn)物種之間的進(jìn)化關(guān)系上優(yōu)于CPF.

3.2 30種甲型流感病毒的神經(jīng)氨酸苷酶基因序列的系統(tǒng)進(jìn)化樹

為了測試AFCS_DFT、DFT和CPF算法聚類神經(jīng)氨酸苷酶全基因序列的準(zhǔn)確性,我們從美國國立生物技術(shù)信息中心NCBI中隨機(jī)選取了30種H1N1-H1N6的神經(jīng)氨酸苷酶全基因序列(NA基因)作為測試序列.根據(jù)AFCS_DFT、DFT和CPF計算得到的DNA序列之間的相似度,應(yīng)用MEGA 6軟件,使用構(gòu)造的0種H1N1-H1N6的神經(jīng)氨酸苷酶全基因序列的系統(tǒng)進(jìn)化樹分別如圖4-圖6所示,圖中的橫坐標(biāo)代表物種相似度的度量的標(biāo)尺.

圖4 AFCS_DFT構(gòu)建的系統(tǒng)進(jìn)化樹Fig.4 Phylogenetic tree constructed by AFCS_DFT

圖5 DFT構(gòu)建的系統(tǒng)進(jìn)化樹Fig.5 Phylogenetic tree constructed by DFT

從圖4-圖6的結(jié)果可知,AFCS_DFT、CPF可以正確地將30種甲型流感病毒H1N1-H1N6的神經(jīng)氨酸苷酶基因序列聚類,而DFT未能將30種甲型流感病毒H1N1-H1N6的神經(jīng)氨酸苷酶基因序列中第2、11、21、24、25號的DNA序列正確聚類,聚類錯誤率達(dá)16.7%.此外,從圖4、圖6可知,與CPF相比,AFCS_DFT更能反映了親緣關(guān)系密切的物種DNA序列之間的差異性較小、進(jìn)化關(guān)系較接近的特征.

這些結(jié)果表明,AFCS_DFTCPF對物種聚類的準(zhǔn)確性高于DFT;在對同一類別物種的聚類中,與CPF相比,AFCS_DFT縮小了親緣關(guān)系較為密切的物種之間DNA序列的差異性,能正確地聚類DNA序列,符合進(jìn)化水平越相近的物種的序列越接近的特征,顯現(xiàn)物種之間的進(jìn)化關(guān)系與類別差異.

圖6 CPF構(gòu)建的系統(tǒng)進(jìn)化樹Fig.6 Phylogenetic tree constructed by CPF

4 結(jié)束語

本文給出的改進(jìn)的AFCS_DFT算法提取了DNA子序列種類、含量和位置3種生物特征,利用提取的特征信息,可以更準(zhǔn)確地計算DNA序列之間的相似度、更準(zhǔn)確地生成系統(tǒng)進(jìn)化樹,更準(zhǔn)確地聚類物種,揭示了生物的進(jìn)化特性.在同一類別聚類中,AFCS_DFT算法縮小了親緣關(guān)系較為密切的物種之間DNA序列的差異性,能夠準(zhǔn)確聚類DNA序列.離散傅里葉變換和對功率譜提取特征向量的計算比較耗時,下一步工作將研究并行化AFCS_DFT算法,以提高效率;以及研究將本文給出的方法應(yīng)用于DNA序列模糊聚類.

[1] Demuth J P,De Bie T,Stajich J E,et al.The evolution of mammalian gene families [J].PLoS One,2006，1(1):e85.doi:10.1371/journal.pone.0000085

[2] Zhao B,Duan V,Yau S S.A novel clustering method via nucleotide-based Fourier power spectrum analysis [J].Journal of Theoretical Biology,2011,279(1):83-89.

[3] Eisen J A.Phylogenomics:improved functional predictions for uncharacterized genes by evolutionary analysis [J].Genome Research,1998,8(3):163-167.

[4] Altschul S F,Gish W,Miller W,et al.Basic local alignment search tool [J].Journal of Molecular Biology,1990,215(3):403-410.

[5] Lipman D J,Pearson W R.Rapid and sensitive protein similarity searches[J].Science,1985,227(4693):1435-1441.

[6] Edgar R C.Search and clustering orders of magnitude faster than BLAST [J].Bioinformatics,2010,26(19):2460-2461.

[7] Pham T D,Zuegg J.A probabilistic measure for alignment-free sequence comparison [J].Bioinformatics,2004,20(18):3455-3461.

[8] Yin Chang-chuan,Stephen S T Yau.An improved model for whole genome phylogenetic analysis by Fourier transform [J].Journal of Theoretical Biology,2015,382:99-110.

[9] Kantorovitz M R,Robinson G E,Sinha S.A statistical method for alignment-free comparison of regulatory sequences [J].Bioinformatics,2007,23(13):249-255.

[10] Freno A.Selecting features by learning markov blankets [C].Proceedings of 11th International Conference on Knowledge-Based Intelligent Informational and Engineering Systems/17th Italian Workshop on Neural Networks,Lecture Notes in Computer Science,2007,4692:69-76.

[11] Lu G,Zhang S,Fang X.An improved string composition method for sequence comparison [J].BMC Bioinformatics,2008,9(Sup.6):S15.

[12] Blaisdell B E.A measure of the similarity of sets of sequences not requiring sequence alignment [J].Proceedings of the National Academy of Sciences of the United States of America,1986,83(14):5155-5159.

[13] Vinga S,Almeida J.Alignment-free sequence comparison-a review[J].Bioinformatics,2003,19(4):513-523.

[14] Wei D,Jiang Q S,Wei Y J,et al.A novel hierarchical clustering algorithm for gene sequences [J].BMC Bioinformatics,2012,13:174.

[15] Bao Jun-peng,Yuan Rui-yu,Bao Zhe.An improved alignment-free model for DNA sequence similarity metric[J].BMC Bioinformatics,2014,15:321.

[16] Yin Chang-chuan,Chen Ying,Stephen S T Yau.A measure of DNA sequnece similarity by Fourier Transform with applications on hierarchical clustering[J].Journal of Theoretical Biology,2014,359(24):18-28.

[17] Hoang Tung,Yin Chang-chuan,Zheng Hui,et al.A new method to cluster DNA sequence using Fourier power spectrum[J].Journal of Theoretical Biology,2015,372:135-145.

[18] Tenreiro Machado J,Costa A C,Quelhas M D.Fractional dynamics in DNA[J].Communications in Nonlinear Science and Numerical Simulation,2011,16 (8):2963-2969.

[19] Oppenheim A V,Schafer R W,Buck.L R.,et al.Discrete-time signal processing [M].Englewood Cliffs:Prentice-hall,2010.

[20] Otu H H,Sayood K.A new sequence distance measure for phylogenetic tree construction[J].Bioinformatics,2003,19 (16):2122-2130.

[21] Feng J,Hu Y,Wan P,et al.New method for comparing DNA primary sequences based on a discrimination measure [J].Journal of Theoretical Biology,2010,266(4):703-707.

[22] Sokal R R.A statistical method for evaluating systematic relationships[J].The University of Kansas Science Bulletin,1958,38:1409-1438.

[23] Cao Y,Janke A,Waddell P J,et al.Conflict among individual mitochondrial proteins in resolving the phylogeny of eutherian orders [J].Journal of Molecular Evolution,1998,47(3):307-322.