檢驗雙重性質(zhì)特征的基因模糊聚類分析方法

祖 穎，朱 平，馬 沖

江南大學(xué) 理學(xué)院，江蘇 無錫 214122

1 引言

序列比較是計算生物學(xué)中最重要的任務(wù)之一，通過它可以判斷生物序列之間的差異或相似性。此外，準(zhǔn)確估計生物序列之間的相似性可以很好地應(yīng)用于相似性分析[1-2]、系統(tǒng)發(fā)育分析[3-4]等?；诒葘Φ姆椒ㄊ切蛄斜容^中使用的主要方法，如Smith-Waterman算法[5]和BLAST[6]，這些方法通常通過動態(tài)規(guī)劃原理得到最佳對齊結(jié)果[7]。然而，這些基于比對的方法在應(yīng)用中存在一些缺點[8]。首先，在利用比對的方法進(jìn)行序列對比時，序列必須在可靠對齊的前提下才能獲得良好的結(jié)果。由于一些序列發(fā)生了基因序列重排，如基因重組或水平基因轉(zhuǎn)移，所以難以獲得準(zhǔn)確地結(jié)果。其次，隨著低成本高通量測序技術(shù)的普及，更多完整的基因組和更多的數(shù)據(jù)可被利用。因此，利用比對的方法進(jìn)行序列對比是非常耗時的。

結(jié)合上述分析，在基于比對的方法之外，也產(chǎn)生了很多不基于比對的方法和應(yīng)用。其中K字詞頻率的方法大多數(shù)被應(yīng)用在快速對序列進(jìn)行物種分類的方面。最常見的基于K字詞頻率和Markov鏈模型的方法包括歐氏距離測量[9]、夾角余弦值測量[10]、標(biāo)準(zhǔn)化歐氏距離測量[9]。Kun等[11]在Torney等[12]、Reinert等[13]和Wan等[14]提出的和統(tǒng)計學(xué)方法基礎(chǔ)上重新加權(quán)定義，并取得了更準(zhǔn)確的聚類結(jié)果?；蛐蛄心承┢蔚闹匾栽谛蛄袑Ρ戎写嬖谥町?。當(dāng)一些片段的順序和頻率幾乎保持不變時，可以合理地假設(shè)這些片段對區(qū)分序列幾乎沒有作用。如果序列片段的順序和頻率有很大的差異，應(yīng)該假設(shè)這些片段對于區(qū)分序列有很大的作用?；谏鲜龇治?，通過利用重新加權(quán)定義的統(tǒng)計學(xué)方法，并將其應(yīng)用于基因序列聚類能夠獲得更準(zhǔn)確的結(jié)果。

核酸堿基對的相互作用是生物結(jié)構(gòu)中的重要作用[15]，DNA結(jié)構(gòu)也通過鏈內(nèi)堿基對相互作用來維持。它在DNA序列的穩(wěn)定性中起著至關(guān)重要的作用。GC堿基對的平均堆疊相互作用比AT堿基對強兩到三倍[16]，與將雙鏈DNA序列保持在一起的其他相互作用（氫鍵和疏水相互作用）相比，DNA序列中堿基對相互作用的能量明顯較大，因此在計算序列之間的距離時，堿基對的相互作用應(yīng)該被考慮進(jìn)去。

在本文中，將4個堿基（A，T，G，C）之間的二階馬爾可夫模型等價于具有16個二核苷酸狀態(tài)的一階馬爾可夫鏈模型，通過計算AA，AT，AG，AC，TA，TT，TG，TC，GA，GT，GG，GC，CA，CT，CG和CC之間的轉(zhuǎn)移概率，同時引入核酸堿基對的相互作用，將它作為一組權(quán)重值，通過模糊聚類分析檢驗雙重性質(zhì)特征的合理性，并利用歐氏距離和夾角余弦值距離這兩種統(tǒng)計學(xué)方法對比單一性質(zhì)和雙重性質(zhì)方法的準(zhǔn)確性。實驗結(jié)果表明雙重性質(zhì)方法更加準(zhǔn)確。

2 方法

2.1 Markov鏈模型

基于Markov鏈模型的轉(zhuǎn)移概率矩陣，一種直接用特征矩陣來刻畫基因序列的方法，具體描述如下。

對于一個給定的長為n的基因序列s=s1s2…sn，基因序列X(s)Markov鏈模型的轉(zhuǎn)移概率矩陣（Sij代表狀態(tài)ij的個數(shù)，Si代表狀態(tài)i的個數(shù)），即由狀態(tài)i轉(zhuǎn)移到狀態(tài) j的概率，因為堿基有4種，那么就有4個可能的狀態(tài)，所以a[i][j]共有4×4種可能的取值。把它們用矩陣表示成：

2.2 核酸堿基對的相互作用

多屬性決策是指在考慮多個屬性的情況下，選擇最優(yōu)備選方案或進(jìn)行方案排序的決策問題。如果第 j個屬性的值在所有決策中變化很小，認(rèn)為這個屬性對決策做出了小的貢獻(xiàn)，應(yīng)該給予一個小的權(quán)重。相反，如果第 j個屬性在所有決策中明顯不同，這個屬性被認(rèn)為在做出最佳決策中發(fā)揮重要作用，應(yīng)該給予很大的重視。在本文中，當(dāng)一個二核苷酸在序列對比作出的貢獻(xiàn)很小時，應(yīng)該分配較小的權(quán)重，當(dāng)一個二核苷酸在序列比對中貢獻(xiàn)很大時，應(yīng)該分配較大的權(quán)重。由于每個二核苷酸具有不同的重要性，考慮核酸堿基對的相互作用作為一組權(quán)重值，再構(gòu)造具有雙重性質(zhì)的特征矩陣，并用模糊聚類分析來檢驗該方法的合理性。

基于上述分析，本文考慮核苷酸兩個性質(zhì)：一個是基于堿基之間的轉(zhuǎn)移概率；另一個是核酸堿基對的相互作用。對于第二個性質(zhì)，在Markov鏈模型的轉(zhuǎn)移概率矩陣中引入核酸堿基對的相互作用，同時定義相互作用[17]：

相互作用=核酸堿基對的相互作用×堿基轉(zhuǎn)移概率

核苷酸之間相互作用的3個概率[17]如下：嘌呤（A，G）轉(zhuǎn)換嘧啶（C，T）概率是1/2，反之亦然；嘌呤轉(zhuǎn)換嘌呤或者嘧啶轉(zhuǎn)換嘧啶的概率是1/3；相同的核苷酸轉(zhuǎn)換的概率是1/6。這種概率的假設(shè)是根據(jù)在物種的進(jìn)化的過程中生物效應(yīng)的轉(zhuǎn)換。

例如從數(shù)據(jù)庫找出一個物種的基因序列片段：

S=AGGCGTTACTGCCTAA

通過Markov鏈模型計算出該序列對應(yīng)的轉(zhuǎn)移概率矩陣為MM，核苷酸相互概率矩陣為MI：

因此，表1中的每個基因序列都可以通過計算得到上述兩種矩陣，使用適當(dāng)?shù)姆椒ㄓ嬎銉尚蛄兄g的距離矩陣。本文使用歐氏距離d(A ,B)[18]，定義如下：

考慮兩個特征來計算序列之間的距離。那么序列A和B之間總的相似性dT(A ,B)[19]定義為：

其中，dM(A,B)表示是堿基的轉(zhuǎn)移概率的相似性，dI(A,B)是堿基的相互作用的相似性。

通過上述距離矩陣，得到16個物種的動態(tài)模糊聚類圖，從而檢驗雙重性質(zhì)方法的合理性，并根據(jù)不同的統(tǒng)計學(xué)方法比較單一性質(zhì)和雙重性質(zhì)方法得到的分層聚類結(jié)果。

3 16個物種的p53基因序列的歐氏距離矩陣

從NCBI數(shù)據(jù)庫里選取了鼢鼠、普通牛、人類等16個物種p53基因完整的CDS序列，來源及序列號見表1。

表1 16條p53基因完整的CDS序列的物種來源及序列號

基于上述的方法，通過Matlab計算表1中所有序列的4×4矩陣 MM和 MI，根據(jù)公式（1）和公式（2）計算序列之間的歐氏距離。利用計算得到的歐氏距離來分析16個物種的p53基因序列的關(guān)系，從而來檢驗雙重性質(zhì)特征方法的有效性。

4 模糊聚類

4.1 相似度計算

基因序列越相似，它們之間的特征值就越小?；诰嚯x方法兩樣本之間關(guān)系的度量時統(tǒng)計學(xué)中通常用的做法是，為了進(jìn)一步構(gòu)造模糊鄰近關(guān)系矩陣，刻畫各個序列之間的關(guān)系采用相似度方法。相似度r用來表示樣本之間相似程度的度量，r越接近0，說明兩個樣本之間的相似度就越低；r越接近1，說明兩個樣本之間的相似度就越高。定義兩基因序列A與B之間的相似度為：

其中，c和α為常數(shù)，這里為使rij∈[ ]0,1，本文選取c=1,α=1對距離矩陣進(jìn)行相似度計算。當(dāng)歐氏距離越大時，rij就越小，得到相應(yīng)的模糊鄰近矩陣 R=見表2。

4.2 模糊動態(tài)聚類圖

基于表2的模糊相似矩陣，用Matlab編程計算得到模糊等價矩陣，從而得到動態(tài)模糊聚類圖，如圖1所示。

從動態(tài)聚類圖1可以看出，當(dāng)選取閾值為λ=0.822 6時可分為2大類：編號1、6、2、4、16、13、11、5、8、3、7、9和14的物種分為一類，其均為陸生類動物；編號10、12和15的物種分為一類，其均為水生類動物。其中每一類物種又可以進(jìn)行分類，當(dāng)閾值為λ=0.901 9，可分為6類：在第一大類中哺乳類動物中人類（編號1）與大猩猩（編號6）關(guān)系最近，聚為一類；牛（編號2）、原始牛（編4）和羊（編號16）同屬于?？撇溉閯游?，聚為一類；家犬（編號11）和狗（編號13）聚為一類；鼠（編號5）、家鼠（編號8）和田鼠（編號14）同屬于倉鼠，聚為一類，彌猴（編號3）、食蟹猴（編號7）、日本獼猴（編號9）聚為一類；歐洲比目魚（編號10）、青鳉（編號12）和紫色球海膽（編號15）聚為一類。根據(jù)不同的需求，選取適當(dāng)?shù)拈撝担傻玫讲煌姆诸惤Y(jié)果。上述聚類結(jié)果符合實際，表明雙重性質(zhì)方法是合理的。

表2 由Markov轉(zhuǎn)移概率和核苷酸相互作用概率模糊相似矩陣

圖1 λ取不同的值時分割樹形圖

通過以上分析，發(fā)現(xiàn)利用核酸堿基對的相互作用，相當(dāng)于對堿基片段AA，AT，AG，AC，TA，TT，TG，TC，GA，GT，GG，GC，CA，CT，CG和CC加權(quán)，然后再構(gòu)造具有雙重性質(zhì)作用的距離矩陣，得到了符合實際的結(jié)果，從而可以合理地推測出堿基片段AA、TT、GG和CC在序列對比中作用較小，相反，堿基片段AC、AT、GC和GT的作用在序列對比中比其他的堿基片段作用都大。

5 層次聚類比較分析

這里用兩種統(tǒng)計學(xué)方法的層次聚類結(jié)果來比較單一性質(zhì)方法和雙重性質(zhì)方法，為了方便計算，先把兩種性質(zhì)矩陣相加，再計算序列之間的距離。利用Matlab軟件編程計算得到16個物種單一性質(zhì)和雙重性質(zhì)的歐氏距離矩陣[9]和余弦角距離矩陣[10]，根據(jù)4個距離矩陣得到聚類結(jié)果如圖2～5所示。

由圖2和圖3可知，它們的區(qū)別在紫色海膽（編號15），紫色海膽屬于水生物種，和歐洲比目魚（編號10）、青鳉（編號12）更親近，在單一性質(zhì)Markov鏈轉(zhuǎn)移概率矩陣下，紫色海膽（編號15）和人類（編號1）、大猩猩（編號6）歐氏距離分別為0.191 7、0.202 3，和歐洲比目魚（編號10）、青鳉（編號12）歐氏距離分別為0.166 7、0.209 0。由于它們歐氏距離差別不明顯，所以圖2中錯誤的將紫色海膽（編號15）和人類（編號1）、大猩猩（編號6）聚為一類。在雙重性質(zhì)距離矩陣下，紫色海膽（編號15）和人類（編號1）、大猩猩（編號6）歐氏距離分別為0.263 1、0.279 5，和歐洲比目魚（編號10）、青鳉（編號12）歐氏距離分別為0.215 7、0.275 9。在雙重性質(zhì)情況下，紫色海膽（編號15）和歐洲比目魚（編號10）、青鳉（編號12）歐氏距離都比人類（編號1）、大猩猩（編號6）歐氏距離要小，所以圖3將它們聚為一類，這表明雙重性質(zhì)方法比單一性質(zhì)方法的歐氏距離聚類結(jié)果效果更好。

圖2 單一性質(zhì)歐氏距離16條基因序列聚類譜系圖

圖3 雙重性質(zhì)歐氏距離16條基因序列聚類譜系圖

圖4 單一性質(zhì)余弦角距離16條基因序列聚類譜系圖

圖5 雙重性質(zhì)余弦角距離16條基因序列聚類譜系圖

由圖4和圖5可以看出，它們的區(qū)別在田鼠（編號14）屬于鼠科動物，和鼠（編號5）、家鼠（編號8）是一類，紫色海膽（編號15）的聚類，紫色海膽（編號15）屬于水生物種，和歐洲比目魚（編號10）、青鳉（編號12）更親近，所以具有雙重性質(zhì)的余弦角距離聚類結(jié)果更符合實際。

綜上所述，結(jié)合兩種距離方法的聚類結(jié)果對比發(fā)現(xiàn)，在提取基因特征中引入核酸堿基對的相互相用的方法比單一方法更準(zhǔn)確，是一種能夠較準(zhǔn)確反映生物學(xué)中一些物種相似性信息的方法。

6 結(jié)論

本文在Markov鏈模型的基礎(chǔ)上引入核酸堿基對相互作用，構(gòu)成了具有雙重性質(zhì)矩陣的方法，采用模糊聚類方法，做出了動態(tài)聚類譜系圖進(jìn)行聚類分析，結(jié)果符合實際；此外，還與單一性質(zhì)的特征矩陣方法作聚類結(jié)果對比，發(fā)現(xiàn)具有雙重特征性質(zhì)矩陣的方法更加準(zhǔn)確。

在本文方法中，考慮了不同二核苷酸的潛在重要性，通過核酸堿基對的相互作用反映這些差異。從實驗結(jié)果可以看出，雙重性質(zhì)方法比單一性質(zhì)方法更加準(zhǔn)確，這告訴人們，從多方面分析事物，才能更加準(zhǔn)確地做出判斷。

如何從二元組擴展到多元組，使核酸堿基對的相互作用得到充分利用，這將是今后的研究工作之一。