潘克邁, 秦雪瑞, 劉雄恩

(福建農(nóng)林大學(xué)計(jì)算機(jī)與信息學(xué)院，福建 福州 350002)

分子系統(tǒng)發(fā)育推算中采用的DNA進(jìn)化模型如JC69、K80、F81、F84、HKY85、TN93、REV94等，均為描述4種核苷酸之間變換的馬爾可夫隨機(jī)過(guò)程[1,2]，用于分子進(jìn)化分析的同源DNA多序列比對(duì)數(shù)據(jù)中出現(xiàn)的比對(duì)缺失(gap)成為上述模型應(yīng)用的障礙.目前通用的數(shù)據(jù)預(yù)處理方法是將含gap位點(diǎn)成對(duì)或整列刪除.DNA序列比對(duì)數(shù)據(jù)代表了初始的同源性假設(shè)，忽略gap的數(shù)據(jù)處理方式會(huì)丟失分子進(jìn)化信息，造成系統(tǒng)發(fā)育樹(shù)枝長(zhǎng)估算的偏差和進(jìn)化距離的偏低估計(jì)[3-5].因此，在分子系統(tǒng)發(fā)育分析推算中融合gap信息是必要的.

Mcguire et al[6]將gap視作核苷酸的第5種狀態(tài)，首次提出了融合gap信息的JC69+gap、F81+gap、F84+gap等5-狀態(tài)模型，基于模型進(jìn)化距離的估算仍會(huì)出現(xiàn)偏差.林碧嬌等[7]在上述3種模型基礎(chǔ)上引入新的參數(shù)，提出JC69+gap′、F81+gap′、F84+gap′等改進(jìn)模型，區(qū)分了indel與substitution在性質(zhì)上的差異.基于改進(jìn)模型推導(dǎo)的核苷酸變換概率5階矩陣可以應(yīng)用于系統(tǒng)發(fā)育重建方法中的最大似然法.改進(jìn)模型F81+gap′減小了進(jìn)化距離估算的偏差，但其中F84+gap′模型的參數(shù)較多、計(jì)算復(fù)雜且未能推導(dǎo)出相應(yīng)的距離計(jì)算公式[8].目前流行的系統(tǒng)發(fā)育重建軟件，如PHYLIP[9]、PAUP[10]、MEGA[11]、PAML[12]、MrBayes[13]、PhyML[14]等，未將上述5-狀態(tài)模型作為計(jì)算模型的參考選項(xiàng).

秦雪瑞等[8]將序列比對(duì)后出現(xiàn)的gap視為統(tǒng)計(jì)抽樣過(guò)程中產(chǎn)生的隨機(jī)缺失數(shù)據(jù)，提出對(duì)同源DNA多序列比對(duì)數(shù)據(jù)中的gap進(jìn)行核苷酸插補(bǔ)等數(shù)據(jù)預(yù)處理方法；給出了核苷酸最近鄰插補(bǔ)算法，并比較分析插補(bǔ)前后序列基于4-狀態(tài)模型及比對(duì)序列基于5-狀態(tài)模型的序列間進(jìn)化距離；對(duì)核苷酸最近鄰插補(bǔ)的數(shù)據(jù)預(yù)處理方法進(jìn)行了評(píng)估.

本文基于最大似然法進(jìn)一步提出DNA序列g(shù)ap核苷酸最大似然插補(bǔ)算法，并比對(duì)缺失位點(diǎn)刪除、核苷酸最近鄰插補(bǔ)以及本文提出的最大似然插補(bǔ)等3種數(shù)據(jù)預(yù)處理方法在進(jìn)化距離估算、系統(tǒng)發(fā)育重建等方面的效果，旨在提高DNA分子系統(tǒng)發(fā)育重建樹(shù)枝長(zhǎng)估算的精確度.

1 方法和原理

1.1 最大似然估計(jì)的一般方法

最大似然法是一種參數(shù)估計(jì)的統(tǒng)計(jì)學(xué)方法，它利用總體的相關(guān)概率密度函數(shù)以及樣本所提供的信息來(lái)求解未知參數(shù)的估計(jì)量.其原理是：從總體中隨機(jī)抽取n個(gè)樣本單元，使樣本觀測(cè)值出現(xiàn)的可能(概率密度函數(shù)似然值)達(dá)到最大的參數(shù)估計(jì)量是最合理的[15].

若離散型總體X的分布函數(shù)為P{X=x}=p(x;θ),θ∈Θ，其中θ是一個(gè)未知參數(shù)或幾個(gè)未知參數(shù)組成的參數(shù)向量，Θ是θ的可能取值范圍.設(shè)X1、X2、 …、Xn是來(lái)自X的樣本，x1、x2、…、xn是對(duì)應(yīng)于樣本X1、X2、…、Xn的樣本觀測(cè)值.已知樣本X1、X2、 …、Xn取到觀察值x1、x2、…、xn的概率，則事件{X1=x1,X2=x2,…,Xn=xn}發(fā)生的概率為：

(1)

這一概率隨θ的取值而變化，L(θ)稱(chēng)為樣本的似然函數(shù).最大似然估計(jì)法固定樣本觀察值x1,x2,…,xn，在θ取值的可能范圍(Θ)篩選使似然函數(shù)L(x1,x2,…,xn;θ)達(dá)到最大參數(shù)值，作為參數(shù)θ的估計(jì)值.即

(2)

為了優(yōu)化計(jì)算，實(shí)際應(yīng)用中通常以對(duì)數(shù)似然函數(shù)的最大值來(lái)估算θ，即

(3)

1.2 比對(duì)缺失位點(diǎn)的核苷酸最大似然插補(bǔ)算法

基于最小進(jìn)化原理[16]，本文假設(shè)同源DNA比對(duì)序列在核苷酸插補(bǔ)前后采用4-狀態(tài)DNA進(jìn)化馬爾可夫模型，與運(yùn)用最大似然法推算的進(jìn)化樹(shù)拓?fù)湟恢?，即二者系統(tǒng)發(fā)育重建樹(shù)的差異僅在于枝長(zhǎng)的估計(jì)值.基于這個(gè)假設(shè)，首先在插補(bǔ)前將gap位點(diǎn)刪除后選擇4-狀態(tài)模型，運(yùn)用上述系統(tǒng)發(fā)育重建的最大似然法建立系統(tǒng)發(fā)育樹(shù)拓?fù)洌⒐浪阒﹂L(zhǎng)和模型的其它參數(shù)；再通過(guò)使單個(gè)位點(diǎn)似然函數(shù)值達(dá)到最大的方法估算gap位點(diǎn)最大可能的核苷酸狀態(tài).

以圖1所示的3個(gè)同源序列系統(tǒng)發(fā)育樹(shù)拓?fù)錇槔?個(gè)序列在一特定位點(diǎn)上的核苷酸狀態(tài)觀察數(shù)據(jù)分別為A、-、C，分枝結(jié)點(diǎn)為0和4.連接到第j個(gè)結(jié)點(diǎn)的枝長(zhǎng)tj定義為平均每個(gè)位點(diǎn)核苷酸置換的期望數(shù)，在重建系統(tǒng)發(fā)育樹(shù)拓?fù)鋾r(shí)已臨時(shí)估算.估算結(jié)點(diǎn)2序列在該缺失位點(diǎn)上最大似然核苷酸的似然函數(shù)：

(4)

式中，πx0是結(jié)點(diǎn)0序列核苷酸的頻率.核苷酸插補(bǔ)的最大似然法需要估計(jì)的參數(shù)是gap位點(diǎn)的核苷酸θ，θ=x,x∈{A,T,C,G}.

同源DNA序列中g(shù)ap位點(diǎn)核苷酸最大似然插補(bǔ)算法描述如下：

Algorithm Nucleotide Interpolation by MLI// maximum likelihood interpolation

Begin

Input multi-aligned DNA sequencesX

GetX′ by deleting all sites with gaps inX

Infer phylogenetic treeT by maximum likelihood withX′

Fors← 1stTo the last gap site inXDo

Begin

k← number of sequences with gap ats

Fori← 1stTo the 4k-th combination of nucleotide for k gaps atsDo

ComputeLiunderT, likelihood value for thei-th combination

IfLh= max (Li) ThenX1..k,s← theh-th combination of nucleotide,X1..k,s∈{T,C,A,G}

Forj← 2ndTo the last sequence Do

IfXj,s≠X1,sThen break Else continue loop

Ifj> count of sequences Then delete sitesElse remains

End

OutputXafter nucleotide interpolation at gap sites

End

與最近鄰插補(bǔ)法一樣，插補(bǔ)后，若各序列在該位點(diǎn)的核苷酸完全相同，則刪除該位點(diǎn)(整列)，否則保留插補(bǔ)后位點(diǎn).剔除插補(bǔ)后核苷酸相同的位點(diǎn)，因?yàn)椴逖a(bǔ)后這種位點(diǎn)在分子進(jìn)化分析中不提供進(jìn)化信息，而且會(huì)減低序列間進(jìn)化距離的估算精度[8].

對(duì)于長(zhǎng)度為n含有s個(gè)gap位點(diǎn)的m個(gè)多序列比對(duì)數(shù)據(jù)矩陣Xm×n，每個(gè)gap位點(diǎn)所在列最多存在m-1個(gè)gap，則核苷酸最大似然法中插補(bǔ)算法的時(shí)間復(fù)雜度為O(4m-1s).

2 結(jié)果與分析

2.1 測(cè)試序列

本文選取3組文獻(xiàn)采用的同源DNA序列[8]，用ClustalX進(jìn)行多序列比對(duì)后，分別采取比對(duì)缺失位點(diǎn)整列刪除(Del)、核苷酸最近鄰插補(bǔ)(MNI)和核苷酸最大似然插補(bǔ)(MLI)等3種數(shù)據(jù)預(yù)處理，估算序列間進(jìn)化距離；并用MEGA中的最大似然法在K80模型下進(jìn)行500次系統(tǒng)發(fā)育重建的測(cè)試試驗(yàn).

利用計(jì)算機(jī)重復(fù)60次模擬，結(jié)果如圖2所示.系統(tǒng)發(fā)育樹(shù)拓?fù)涞腄NA進(jìn)化過(guò)程，其中一組僅有substitution事件，另一組包含indel事件.DNA分子進(jìn)化的計(jì)算機(jī)模擬是基于分子鐘假說(shuō)[17,18]，采取固定時(shí)序隨機(jī)抽取位點(diǎn)和沿樹(shù)進(jìn)化序列的方法[19].模擬過(guò)程記錄真實(shí)樹(shù)各分枝枝長(zhǎng)以及8個(gè)同源DNA終端序列.30組僅有substitution事件的DNA多序列比對(duì)數(shù)據(jù)不含gap，在無(wú)需數(shù)據(jù)預(yù)處理的情形下重建系統(tǒng)發(fā)育樹(shù)，以重建樹(shù)枝長(zhǎng)與真實(shí)樹(shù)枝長(zhǎng)的誤差作為基準(zhǔn)誤差；另30組記錄包含indel事件的同源序列分別在3種預(yù)處理下重建樹(shù)的枝長(zhǎng)；并分析不同預(yù)處理方法在系統(tǒng)發(fā)育重建時(shí)的誤差.

圖1 系統(tǒng)發(fā)育樹(shù)拓?fù)淠澄稽c(diǎn)示意圖Fig.1 A phylogenetic tree of three homologous DNA sequences at a site

2.2 測(cè)試結(jié)果

2.2.1 不同預(yù)處理方法的進(jìn)化距離 表1、2、3分別顯示第1、2、3組文獻(xiàn)數(shù)據(jù)[8]在Del、MNI和MLI數(shù)據(jù)預(yù)處理下序列間進(jìn)化距離的估算結(jié)果.

表1 猿類(lèi)7個(gè)物種線(xiàn)粒體DNA序列在不同預(yù)處理和模型下的進(jìn)化距離Table 1 Evolutionary distances of mitochondrial DNA sequences of 7 apes under different preprocessing and models

表2 睡蓮科6種植物核糖體DNA中ITS在不同預(yù)處理和模型下的進(jìn)化距離Table 2 Evolutionary distances of rDNA ITS2 region of 6 Nymphaeaceae plants under different preprocessing and models

2.2.2 不同預(yù)處理方法下的系統(tǒng)發(fā)育重建 圖3、4、5分別顯示第1、2、3組文獻(xiàn)數(shù)據(jù)[8]在Del、NNI和MLI等3種數(shù)據(jù)預(yù)處理下系統(tǒng)發(fā)育重建的結(jié)果(分別位于圖的上、中、下部分).

從表1～3可知：比對(duì)缺失位點(diǎn)的核苷酸插補(bǔ)，無(wú)論是MNI算法還是MLI算法預(yù)處理后的序列間進(jìn)化距離都大于Del預(yù)處理下的距離估算.在同源多序列比對(duì)序列矩陣中直接刪除含有g(shù)ap位點(diǎn)的列，突變信息丟失是造成其序列間進(jìn)化距離偏低估計(jì)的原因.

從圖3～5可知：MNI和MLI核苷酸插補(bǔ)的預(yù)處理方法下系統(tǒng)發(fā)育重建的進(jìn)化樹(shù)拓?fù)渑cDel預(yù)處理下重建進(jìn)化樹(shù)的拓?fù)涫且恢碌模哉箼z驗(yàn)結(jié)果具有很高且?guī)缀跻恢碌闹眯哦?其次，證實(shí)了核苷酸插補(bǔ)前后基于4-狀態(tài)DNA進(jìn)化馬爾可夫模型的系統(tǒng)發(fā)育重建進(jìn)化樹(shù)具有相同的拓?fù)浣Y(jié)構(gòu)，差異僅在于枝長(zhǎng)的估計(jì)值.

表3 側(cè)耳屬8種真菌25S rRNA序列在不同預(yù)處理和模型下的進(jìn)化距離Table 3 Evolutionary distances of 25S rRNA sequences of 8 Pleurotus fungus under different preprocessing and models

圖3 第1組數(shù)據(jù)在3種預(yù)處理后重建的系統(tǒng)發(fā)育樹(shù)Fig.3 Reconstructed phylogenetic trees from processed data of group 1 by 3 preprocessing methods

圖4 第2組數(shù)據(jù)在3種預(yù)處理后重建的系統(tǒng)發(fā)育樹(shù)Fig.4 Reconstructed phylogenetic trees from processed data of group 2 by 3 preprocessing methods

圖5 第3組數(shù)據(jù)在3種預(yù)處理后重建的系統(tǒng)發(fā)育樹(shù)Fig.5 Reconstructed phylogenetic trees from processed data of group 3 by 3 preprocessing methods

2.2.3 模擬系統(tǒng)發(fā)育的真實(shí)樹(shù)與重建樹(shù)枝長(zhǎng)的誤差 3種預(yù)處理下模擬同源序列系統(tǒng)發(fā)育重建樹(shù)的枝長(zhǎng)與模擬真實(shí)樹(shù)枝長(zhǎng)的誤差統(tǒng)計(jì)結(jié)果(表4)表明：MNI和MLI插補(bǔ)的預(yù)處理方法下的枝長(zhǎng)誤差，無(wú)論是絕對(duì)誤差還是均方根誤差均小于Del預(yù)處理下的枝長(zhǎng)誤差，這是直接刪除含有g(shù)ap位點(diǎn)丟失進(jìn)化信息所導(dǎo)致的.而MNI和MLI兩者的枝長(zhǎng)誤差較小，且基本相同.經(jīng)統(tǒng)計(jì)，在30組模擬序列的MNI和MLI之間存在7.556%±2.366% gap位點(diǎn)插補(bǔ)核苷酸的差異,這種差異在不同分枝枝長(zhǎng)估算結(jié)果上體現(xiàn)出來(lái).

表4 重復(fù)模擬的系統(tǒng)發(fā)育樹(shù)枝長(zhǎng)在3種數(shù)據(jù)預(yù)處理下的誤差Table 4 Deviation of estimated branch length of simulated phylogenetic trees under 3 preprocessing methods

3 小結(jié)

真實(shí)序列和模擬序列試驗(yàn)均證實(shí)了MNI和MLI兩種核苷酸插補(bǔ)算法能夠有效減小序列間進(jìn)化距離的偏低估計(jì)；在系統(tǒng)發(fā)育重建中，兩種核苷酸插補(bǔ)的預(yù)處理方法能夠保持重建進(jìn)化樹(shù)的拓?fù)渑c流行的預(yù)處理方法一致的結(jié)構(gòu)，并能提高枝長(zhǎng)估算的精確度.而本文提出的核苷酸最大似然插補(bǔ)法與最近鄰插補(bǔ)法在枝長(zhǎng)估計(jì)誤差的減小程度上基本相同，雖然最近鄰插補(bǔ)的算法簡(jiǎn)捷，其時(shí)間復(fù)雜度[O(ms)]低于最大似然插補(bǔ)算法[O(4m-1s)]，但鑒于最大似然法良好的統(tǒng)計(jì)學(xué)性質(zhì)，推測(cè)其插補(bǔ)預(yù)處理后的數(shù)據(jù)更可靠.