祖 穎，朱 平，馬 沖

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

1 引言

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

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

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

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

2 方法

2.1 Markov鏈模型

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

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

2.2 核酸堿基對(duì)的相互作用

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

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

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

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

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

S=AGGCGTTACTGCCTAA

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

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

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

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

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

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

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

表1 16條p53基因完整的CDS序列的物種來(lái)源及序列號(hào)

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

4 模糊聚類

4.1 相似度計(jì)算

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

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

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

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

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

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

圖1 λ取不同的值時(shí)分割樹(shù)形圖

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

5 層次聚類比較分析

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

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

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

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

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

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

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

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

6 結(jié)論

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

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

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