艾 亮,馮 杰
(中央民族大學(xué) 理學(xué)院,北京 100081)
生物序列的相似性分析是生物信息學(xué)的重要研究方向之一。在早期研究中,通常采用多序列比對(duì)的方法對(duì)序列進(jìn)行比較分析,許多算法現(xiàn)在已經(jīng)非常成熟[1-3],例如使用較多的ClustalW算法。但多序列比對(duì)是基于同源序列片段間是鄰接保守的假設(shè),這與遺傳重組相沖突,而且當(dāng)樣本量較大或序列長(zhǎng)度較長(zhǎng)時(shí),比對(duì)算法的時(shí)間成本很高。因此,非比對(duì)方法[4]一經(jīng)推出,立即受到研究人員的廣泛關(guān)注。非比對(duì)方法不是具體比較基對(duì),而是將序列看成是一個(gè)整體并將其轉(zhuǎn)化為數(shù)值向量再進(jìn)行分析比較,其優(yōu)點(diǎn)是在計(jì)算機(jī)上計(jì)算迅速,且結(jié)果較準(zhǔn)確。
蛋白質(zhì)序列的比較分析方法大致分為兩大類:圖形表示方法和數(shù)值向量刻畫方法。圖形表示方法也稱為可視化方法,其基本思想是建立一組映射,將氨基酸映射成平面或空間的點(diǎn),然后將點(diǎn)連接起來(lái)得到空間曲線。進(jìn)一步地,我們還可以從這些圖形表示中提取生物序列的數(shù)值特征,利用這些數(shù)值特征進(jìn)行序列分析[5-12]。數(shù)值向量刻畫方法主要是將蛋白質(zhì)序列轉(zhuǎn)換為多維的數(shù)值向量,例如Chou K.C.[13]和Chen W.等[14]將氨基酸的20維頻率向量與理化性質(zhì)或者氨基酸之間的相互作用結(jié)合起來(lái)構(gòu)建(20+λ)維向量來(lái)表示蛋白質(zhì)序列,其中λ指的是理化性質(zhì)個(gè)數(shù)或者氨基酸相互之間作用的指標(biāo)數(shù)。賈美多等[15]結(jié)合氨基酸的5-字母分類模型和序列的k-字節(jié)模型,提取信息將序列轉(zhuǎn)化為一個(gè)30維向量,之后利用歐氏距離求得蛋白質(zhì)序列兩兩間的相對(duì)距離進(jìn)而構(gòu)建系統(tǒng)進(jìn)化樹。Xian-HuaXie等[16]使用氨基酸隨機(jī)和獨(dú)立放置的序列分布圖之間的相對(duì)偏差來(lái)定義序列間的差異。Li Y等[17]結(jié)合氨基酸的概率、平均出現(xiàn)位置概率和兩個(gè)相鄰氨基酸的馬爾科夫轉(zhuǎn)移概率分布來(lái)構(gòu)建蛋白質(zhì)數(shù)值向量表示。Yongkun Li等[18]將蛋白質(zhì)序列中20種氨基酸的數(shù)目、平均位置和位置的正則化中心二階矩結(jié)合起來(lái)構(gòu)成60維數(shù)值向量來(lái)衡量病毒之間的相似性。Lily He等[19]基于氨基酸的親水性指數(shù)、極性需求和側(cè)鏈的化學(xué)成分將氨基酸分成8類,之后融合蛋白質(zhì)序列中這8類氨基酸的數(shù)量、平均位置和位置的二階矩信息構(gòu)建24維特征向量進(jìn)行進(jìn)化分析。朱臣臣等[20]選擇3種氨基酸的理化性質(zhì)繪制蛋白質(zhì)序列的3D圖形,再基于氨基酸的位置信息構(gòu)建20個(gè)點(diǎn)集,分別求其轉(zhuǎn)動(dòng)慣量構(gòu)建23維特征向量。Stephen S.-T. Yau等[21]、Qi Dai等[22]和Yufeng Liu等[23]統(tǒng)計(jì)序列中所有的長(zhǎng)度為k的子串的頻率,將這些數(shù)字組成向量,使用該向量刻畫生物序列的特征。
蛋白質(zhì)是由氨基酸組成的,已有研究表明氨基酸的物理化學(xué)性質(zhì)對(duì)蛋白質(zhì)序列分類和進(jìn)化具有重要意義[24-25]。本文將氨基酸的10種理化性質(zhì)通過(guò)主成分分析濃縮為6組主成分,對(duì)每條蛋白質(zhì)序列,計(jì)算反映氨基酸理化性質(zhì)的6組主成分得分均值,再融合20個(gè)氨基酸的位置信息構(gòu)成一個(gè)26維的蛋白質(zhì)序列特征向量,最后利用歐式距離度量蛋白質(zhì)序列間的相似性并構(gòu)造系統(tǒng)進(jìn)化樹。通過(guò)對(duì)三組蛋白質(zhì)序列數(shù)據(jù)集的測(cè)試表明,本文的方法能將每條蛋白質(zhì)序列準(zhǔn)確聚類,結(jié)果與現(xiàn)有進(jìn)化關(guān)系一致,說(shuō)明了該方法的有效性。


表1 氨基酸的10種理化性質(zhì)Table 1 10 physicochemical properties of amino acids
為消除量綱不一的影響,先對(duì)10組氨基酸理化性質(zhì)進(jìn)行標(biāo)準(zhǔn)化處理,將其化為均值為0,標(biāo)準(zhǔn)差為1的數(shù)據(jù)框。然后對(duì)該20×10的氨基酸理化性質(zhì)矩陣進(jìn)行主成分分析,將10組變量的信息壓縮為幾個(gè)綜合變量,提取有效的主成分來(lái)表示20種天然氨基酸的理化性質(zhì)。主成分分析結(jié)果見表2。

表2 重要主成分的貢獻(xiàn)率Table 2 Contribution of significant principal components
由表2可以看到,前6個(gè)主成分的累計(jì)貢獻(xiàn)率為95.91%,遠(yuǎn)大于85%,可以認(rèn)為這6個(gè)主成分能代表原先10組理化性質(zhì)的絕大部分信息。計(jì)算這6個(gè)主成分的得分,即把原來(lái)的20×10的氨基酸理化性質(zhì)矩陣轉(zhuǎn)化為20×6的主成分得分矩陣(見表3)。

表3 主成分得分矩陣Table 3 Principal component score matrix

(1)


對(duì)于一條長(zhǎng)度為n的蛋白質(zhì)序列S=(s1,s2,…,sn),其中sj∈Ω,j=1,2,…,n,還可以基于每種氨基酸Ai(i=1,2,…,20)的平均位置[19]構(gòu)造一個(gè)20維的特征向量,如下所示:
(2)


為驗(yàn)證本文所提方法的有效性,用三組蛋白質(zhì)序列數(shù)據(jù)集[19- 20]進(jìn)行實(shí)驗(yàn),利用歐氏距離計(jì)算兩兩蛋白質(zhì)序列所對(duì)應(yīng)的26維特征向量之間的距離,然后利用UPGMA算法(該算法已嵌入到MEGA11軟件)構(gòu)建生物系統(tǒng)進(jìn)化樹。
9個(gè)物種的ND5蛋白質(zhì)序列信息在表4中給出。使用本文的方法,可以得到9物種ND5蛋白質(zhì)序列的一個(gè)9×26特征矩陣,然后計(jì)算兩兩間的歐氏距離可以得到相似性距離矩陣,結(jié)果見表5。

表4 9物種ND5蛋白質(zhì)序列信息Table 4 Information on 9 ND5 protein sequences

表5 9物種ND5蛋白質(zhì)序列相似性距離矩陣Table 5 The similarity/dissimilarity matrix of 9 ND5 protein sequences
觀察表5可以看出,普通黑猩猩和侏儒黑猩猩的相似性距離最小,為2.679,表示普通黑猩猩和侏儒黑猩猩間的親緣關(guān)系最近;大鼠和負(fù)鼠的相似性距離最大,為9.976,表示大鼠和負(fù)鼠間親緣關(guān)系最遠(yuǎn)。同時(shí),可以看到,人類、普通黑猩猩、侏儒黑猩猩和大猩猩這四個(gè)物種間的相似性距離比較小,說(shuō)明它們的蛋白質(zhì)序列相似性程度高,進(jìn)化關(guān)系上較為接近;長(zhǎng)須鯨和藍(lán)鯨間相似性距離也很小,說(shuō)明它們的進(jìn)化關(guān)系接近;負(fù)鼠和其他八個(gè)物種的相似性距離都很大,表明在進(jìn)化關(guān)系上與其它物種相比負(fù)鼠相對(duì)比較獨(dú)立。
進(jìn)一步利用相似性距離矩陣構(gòu)建物種進(jìn)化樹,結(jié)果如圖1所示。通過(guò)觀察發(fā)現(xiàn)9個(gè)物種被分成4個(gè)分支:第1個(gè)分支是侏儒黑猩猩、普通黑猩猩、人類和大猩猩,在這一分支中,侏儒黑猩猩和普通黑猩猩進(jìn)化關(guān)系更近,其次是人類,而后是大猩猩,這與進(jìn)化事實(shí)相符合;第2個(gè)分支是藍(lán)鯨和長(zhǎng)須鯨;第3分支為大鼠和小鼠;第4個(gè)分支為負(fù)鼠,與其他物種進(jìn)化關(guān)系較遠(yuǎn),單獨(dú)成一個(gè)分支。從進(jìn)化關(guān)系上看,侏儒黑猩猩、普通黑猩猩、人類和大猩猩都屬于靈長(zhǎng)目人科,藍(lán)鯨和長(zhǎng)須鯨都屬于鯨目須鯨科,大鼠和小鼠都屬于嚙齒目鼠科,負(fù)鼠屬于負(fù)鼠目負(fù)鼠科,本文的分析結(jié)果與實(shí)際進(jìn)化關(guān)系相一致。

圖1 9物種ND5蛋白質(zhì)序列的進(jìn)化樹Fig.1 The phylogenetic tree of 9 ND5 protein sequences
12個(gè)桿狀病毒蛋白質(zhì)序列信息見表6,使用本文所提方法對(duì)其構(gòu)建進(jìn)化樹,結(jié)果見圖2。由圖2可以看到,Alphabaculovirus病毒和Betabaculovirus病毒被分為兩大分支,并且Alphabaculovirus中的Group I和Group II也都形成各自的分支,與實(shí)際的病毒進(jìn)化關(guān)系一致。而文獻(xiàn)[7]沒(méi)有將Alphabaculovirus病毒和Betabaculovirus病毒形成兩個(gè)大分支,并且Group II中的6個(gè)病毒不在一個(gè)分支,HzSNPV、HaSNPV、HearNPV與Betabaculovirus病毒的進(jìn)化距離要比與Group II中的其他病毒的距離要小,這與實(shí)際的進(jìn)化關(guān)系不一致。文獻(xiàn)[8]和[20]雖然將Alphabaculovirus病毒和Betabaculovirus病毒形成了兩個(gè)大分支,但是Group II中的6個(gè)病毒并不在一個(gè)分支,Group I中的AcMNPV、BmNPV、RoMNPV各自與Group II中的三個(gè)病毒形成分支。

圖2 12個(gè)桿狀病毒蛋白質(zhì)序列的進(jìn)化樹Fig.2 The phylogenetic tree of 12 baculovirus protein sequences

表6 12個(gè)桿狀病毒蛋白質(zhì)序列信息Table 6 Information on 12 Baculovirus protein sequences
甲型流感病毒的一些亞型是根據(jù)H(血凝素類型)的編號(hào)(H1到H18)和N(神經(jīng)氨酸酶類型)的編號(hào)(N1到N11)來(lái)標(biāo)記的,最致命的甲流亞型是H1N1、H2N2、H5N1、H7N3和H7N9,本文選取了35個(gè)與這些重要亞型相關(guān)的蛋白質(zhì)序列。
使用我們的方法對(duì)該蛋白質(zhì)序列數(shù)據(jù)集構(gòu)建進(jìn)化樹,結(jié)果見圖3。由圖3可知,五種最致命的甲型流感病毒亞型H1N1、H2N2、H5N1、H7N3和H7N9各自形成5個(gè)分支,35個(gè)病毒都被正確聚類。相比之下,用ClustalW方法構(gòu)建的進(jìn)化樹則有3個(gè)甲型流感病毒亞型聚類錯(cuò)誤,如圖4所示,其中A/turkey/VA/505477-18/2007(H5N1),A/turkey/Ontario/FAV110-4/2009(H1N1)和A/turkey/Virginia/4135/2014(H1N1)沒(méi)能正確被聚類。并且在同一臺(tái)筆記本電腦下,ClustalW方法完成多序列比對(duì)需要花費(fèi)約7 s,而我們的方法將序列轉(zhuǎn)化為特征向量只需0.17 s。

圖3 本文方法構(gòu)建的35個(gè)甲型流感病毒蛋白質(zhì)序列的進(jìn)化樹Fig.3 The phylogenetic tree of 35 influenza A virus protein sequences constructed using our method

圖4 ClustalW方法構(gòu)建的35個(gè)甲型流感病毒蛋白質(zhì)序列的進(jìn)化樹Fig.4 The phylogenetic tree of 35 influenza A virus protein sequences constructed using ClustalW
新的非比對(duì)的蛋白質(zhì)序列相似性分析的方法,將蛋白質(zhì)序列轉(zhuǎn)化為數(shù)值向量時(shí),同時(shí)考慮了蛋白質(zhì)序列中20種天然氨基酸的數(shù)量、理化性質(zhì)和平均位置信息,最終將每條蛋白質(zhì)序列都轉(zhuǎn)化為唯一與之對(duì)應(yīng)的26維特征向量。該新方法在3個(gè)數(shù)據(jù)集上均獲得了準(zhǔn)確的聚類結(jié)果,這說(shuō)明該新方法在分析蛋白質(zhì)序列的相似性方面是有效的。此外,該方法不需要復(fù)雜的計(jì)算,而且簡(jiǎn)便快捷。