分子進(jìn)化研究中多重突變的校正方法

李可群

（同濟(jì)大學(xué)化學(xué)科學(xué)與工程學(xué)院，上海200092）

準(zhǔn)確地估計(jì)物種分歧時(shí)間和推斷其進(jìn)化歷史是分子進(jìn)化遺傳學(xué)和系統(tǒng)發(fā)育學(xué)的一個(gè)重要課題[1]。目前分子系統(tǒng)發(fā)育分析在計(jì)算物種分歧時(shí)間時(shí)大多基于分子進(jìn)化速率恒定的“分子鐘”假說(shuō)，但絕大多數(shù)生物分子在長(zhǎng)時(shí)間尺度和不同譜系的進(jìn)化速率并不恒定，從而計(jì)算得到的結(jié)果與化石年齡往往存在較大的偏差。如計(jì)算得到的哺乳動(dòng)物中真獸類(lèi)的物種分歧時(shí)間平均比化石早50%～90%。基于分子數(shù)據(jù)的系統(tǒng)發(fā)育研究方法還不能很好地闡明哺乳動(dòng)物基部問(wèn)題以及哺乳動(dòng)物進(jìn)化過(guò)程所經(jīng)歷的快速輻射事件[2]。為此，本文作者提出了一種不基于進(jìn)化速率恒定假說(shuō)的物種分歧時(shí)間計(jì)算方法[3]，并使用它發(fā)現(xiàn)了蛋白質(zhì)和核苷酸序列分子的絕對(duì)進(jìn)化速率與其進(jìn)化時(shí)間之間存在定量關(guān)系，即“活化能公式”[4]。我們還計(jì)算出了“寒武紀(jì)生物大爆發(fā)”時(shí)期一些主要物種類(lèi)群的物種分歧時(shí)間，所得結(jié)果與化石年齡很接近，明顯好于現(xiàn)有基于進(jìn)化速率恒定假說(shuō)的分子鐘方法得到的結(jié)果[5]。本文將繼續(xù)對(duì)分子進(jìn)化研究中兩種主要多重突變類(lèi)型，即回復(fù)突變和平行突變的校正方法進(jìn)行探討。

1 分子進(jìn)化與泊松分布

1.1 位點(diǎn)突變、序列分子突變與泊松分布

泊松分布一般可用來(lái)描述稀有隨機(jī)事件[6]，當(dāng)一個(gè)事件以固定的速率隨機(jī)且獨(dú)立地發(fā)生時(shí)，這個(gè)事件在某一時(shí)間段內(nèi)發(fā)生的次數(shù)可用泊松分布來(lái)描述，它的概率密度函數(shù)為。若兩個(gè)獨(dú)立的泊松分布并列進(jìn)行或先后發(fā)生，它們的平均發(fā)生次數(shù)分別為λ1和λ2，則總的隨機(jī)事件仍為泊松分布，總的隨機(jī)事件的平均發(fā)生次數(shù)為兩個(gè)獨(dú)立泊松分布平均發(fā)生次數(shù)之和，即λ=λ1+λ2，這一性質(zhì)被稱(chēng)為泊松分布的可加性[7]。

核苷酸和蛋白質(zhì)序列分子的位點(diǎn)突變率都很低，如蛋白質(zhì)序列分子，其位點(diǎn)突變速率的數(shù)量級(jí)大多約為10-9位點(diǎn)/年[8]。一般核苷酸和蛋白質(zhì)序列分子在較小的時(shí)間尺度內(nèi)可視為進(jìn)化速率恒定，因此在一個(gè)突變速率恒定的小時(shí)間段ta內(nèi)，核苷酸和蛋白質(zhì)序列分子單個(gè)位點(diǎn)的突變次數(shù)可認(rèn)為服從泊松分布：

式（1）為單個(gè)核苷酸或蛋白質(zhì)序列分子位點(diǎn)在該小時(shí)間段發(fā)生j次突變的概率，ka為其絕對(duì)進(jìn)化速率，kata為平均突變次數(shù)，分子進(jìn)化中被稱(chēng)之為遺傳距離。特別地，單個(gè)核苷酸或蛋白質(zhì)序列分子位點(diǎn)不發(fā)生突變的概率為p(x=0)=e-kata。對(duì)于一個(gè)有n0個(gè)被比較位點(diǎn)的核苷酸或蛋白質(zhì)序列分子，若不考慮回復(fù)突變，有[3]

式（2）中d為核苷酸或蛋白質(zhì)序列分子相對(duì)于其被比較祖先序列分子的差異位點(diǎn)數(shù)，k為它們的絕對(duì)進(jìn)化速率，t為自其祖先序列分子而來(lái)的進(jìn)化時(shí)間，p為這兩個(gè)核苷酸或蛋白質(zhì)序列分子的序列差異率。

1.2 分子進(jìn)化與泊松分布可加性

核苷酸和蛋白質(zhì)序列分子進(jìn)化速率在長(zhǎng)時(shí)間尺度內(nèi)并不恒定，序列分子各位點(diǎn)的突變速率也存在差異[9]。為了研究方便，我們把核苷酸序列分子的進(jìn)化路徑拆分成由很多個(gè)足夠小的時(shí)間段組成，以便在每個(gè)小時(shí)間段內(nèi)核苷酸序列分子上各位點(diǎn)的突變速率分別恒定。同時(shí)分子進(jìn)化和分子系統(tǒng)發(fā)育分析中一般認(rèn)為生物序列分子的位點(diǎn)突變之間是獨(dú)立的[10]，這樣便得到一個(gè)所謂的“泊松分布矩陣”。由于泊松分布具有可加性，“泊松分布矩陣”中同一行，即每個(gè)小時(shí)間段核苷酸序列分子各位點(diǎn)的泊松分布相加可得到一個(gè)新的泊松分布，新的泊松分布的遺傳距離等于組成它的各位點(diǎn)泊松分布的遺傳距離之和；而“泊松分布矩陣”中每一行相加而成的新泊松分布還可縱向逐個(gè)相加，它們的遺傳距離也一樣可以縱向逐個(gè)相加，后一相加關(guān)系即為文獻(xiàn)[3]中所提及的遺傳距離加和公式,或稱(chēng)遺傳距離可加性。“泊松分布矩陣”中的每一列，即核苷酸序列分子的每個(gè)位點(diǎn)，它們?cè)谶M(jìn)化路徑上各小時(shí)間段的泊松分布可以相加，后者的遺傳距離亦可相應(yīng)地相加。相加而成的各位點(diǎn)泊松分布還可相加成一個(gè)總的泊松分布，它們的遺傳距離也可相加。這樣，我們所研究的核苷酸序列分子在其進(jìn)化路徑上累積發(fā)生的突變可用該總的泊松分布來(lái)準(zhǔn)確描述。

蛋白質(zhì)序列分子的進(jìn)化情形與核苷酸序列分子相似。不過(guò)兩者也有不同的地方，編碼序列中三個(gè)連續(xù)的核苷酸分子組成的密碼子對(duì)應(yīng)一個(gè)氨基酸，其中引起氨基酸變化的叫非同義密碼子，不引起氨基酸變化的密碼子叫同義密碼子。蛋白質(zhì)序列分子的“泊松分布矩陣”與核苷酸序列分子的“泊松分布矩陣”相比，缺失由同義密碼子組成的那部分泊松分布（它們的遺傳距離為零）。因此，對(duì)于相同的進(jìn)化時(shí)間段，蛋白質(zhì)序列分子的進(jìn)化速率一般小于其對(duì)應(yīng)的核苷酸序列分子進(jìn)化速率。

2 物種分歧時(shí)間計(jì)算中回復(fù)突變和平行突變的校正方法

2.1 蛋白質(zhì)序列分子中回復(fù)突變和平行突變的校正方法

物種分歧時(shí)間計(jì)算一般多使用蛋白質(zhì)序列分子。我們首先討論蛋白質(zhì)序列分子中回復(fù)突變和平行突變的校正。

計(jì)算生物物種分歧時(shí)間時(shí)，因祖先蛋白質(zhì)序列分子難以得到，所以一般通過(guò)比較現(xiàn)生同源蛋白質(zhì)序列分子來(lái)計(jì)算分子絕對(duì)進(jìn)化速率和物種分歧時(shí)間。文獻(xiàn)[3]給出了圖1所示的兩個(gè)同源蛋白質(zhì)序列分子與它們共同祖先蛋白質(zhì)序列分子之間的遺傳三角形關(guān)系，一個(gè)蛋白質(zhì)序列分子所有已突變氨基酸位點(diǎn)與其同源蛋白質(zhì)序列分子相同位點(diǎn)相比等同于與它們祖先蛋白質(zhì)序列分子相同位點(diǎn)相比。

若遺傳三角形中兩個(gè)同源蛋白質(zhì)序列分子的絕對(duì)進(jìn)化速率分別為k1和k2,與它們共同祖先序列分子的序列差異率分別為p′1和p′2,即和。另設(shè)兩個(gè)同源蛋白質(zhì)序列分子自它們共同祖先序列分子的物種分歧時(shí)間為t，則根據(jù)式（2）有

圖1 [3] 回復(fù)突變和平行突變校正方法的示意圖

而兩個(gè)同源蛋白質(zhì)序列分子共同突變位點(diǎn)的突變概率為p′1p′2=(1-e-k1t)(1-e-k2t)。兩個(gè)同源蛋白質(zhì)序列分子比較時(shí)需考慮以下兩種突變類(lèi)型的校正，即需考慮它們對(duì)序列差異率的影響：

（1）回復(fù)突變：同源蛋白質(zhì)序列分子的氨基酸位點(diǎn)多次突變后形成與其祖先序列分子同一位點(diǎn)相同的氨基酸。若一個(gè)蛋白質(zhì)序列分子發(fā)生回復(fù)突變后，其同源蛋白質(zhì)序列分子對(duì)應(yīng)的相同位點(diǎn)未發(fā)生突變，則會(huì)使兩者的序列差異率減少，故需進(jìn)行校正。

由于存在20種基本氨基酸，若它們彼此突變的概率相同，兩同源蛋白質(zhì)序列分子的回復(fù)突變概率分別為式（3）中各自突變概率的。不難理解，兩個(gè)同源蛋白質(zhì)序列分子中一個(gè)蛋白質(zhì)序列分子發(fā)生回復(fù)突變，而另一個(gè)蛋白質(zhì)序列分子相同位點(diǎn)未發(fā)生突變的概率分別為p1=0.05(1-e-k1t)×e-k2t，p2=0.05(1-e-k2t)×e-k1t。兩個(gè)同源蛋白質(zhì)序列分子相同位點(diǎn)同時(shí)發(fā)生回復(fù)突變的情形放在平行突變校正中考慮。

（2）平行突變：兩個(gè)同源蛋白質(zhì)序列分子的相同位點(diǎn)同時(shí)突變成同一種氨基酸。平行突變也會(huì)使兩個(gè)同源蛋白質(zhì)序列分子的序列差異率減少。平行突變的概率p″為

式（4）中后兩個(gè)乘積項(xiàng)分別為兩個(gè)同源蛋白質(zhì)序列分子的相同位點(diǎn)同時(shí)突變成某一種相同氨基酸的概率，它們分別為式（3）中各自序列分子突變率的。而式（4）中第一個(gè)乘積項(xiàng)表示有20種可能形成的氨基酸。兩個(gè)同源蛋白質(zhì)序列分子的回復(fù)突變和平行突變導(dǎo)致兩者序列差異率減少之和為

由文獻(xiàn)[3]可知，未校正時(shí)，兩個(gè)同源蛋白質(zhì)序列分子的差異率為1-e-(k1+k2)t，故兩個(gè)同源蛋白質(zhì)序列分子在校正回復(fù)突變和平行突變后的序列差異率p為兩者之差，即

可以看出，需進(jìn)行校正的回復(fù)突變和平行突變概率之和為其序列差異率的。式（5）為校正回復(fù)突變和平行突變后物種分歧時(shí)間的計(jì)算公式，或多重突變校正后圖1中遺傳三角形的定量關(guān)系式。由文獻(xiàn)[3]可知，在該文獻(xiàn)表1中不引入較大相對(duì)誤差的情況下，使用替代式（5）中e-(k1+k2)t，可取得較好的計(jì)算效果，即式（5）可用0.475(e-k1t+e-k2t)=0.95-p替代。

蛋白質(zhì)序列決定蛋白質(zhì)的結(jié)構(gòu)，蛋白質(zhì)結(jié)構(gòu)又決定蛋白質(zhì)的性能。文獻(xiàn)[11]指出，氨基酸的親水性和疏水性對(duì)蛋白質(zhì)性能的影響很大，蛋白質(zhì)序列分子突變的一般性規(guī)律是：同為親水性或疏水性氨基酸之間突變較容易，而親水性和疏水性氨基酸之間的突變較難發(fā)生。組成蛋白質(zhì)的20種基本氨基酸中親水氨基酸和疏水氨基酸各占約10種[11-12]，因此在實(shí)際計(jì)算中發(fā)現(xiàn)，同源蛋白質(zhì)序列分子比較時(shí)需校正的回復(fù)突變和平行突變概率之和取為其序列差異率的可取得較好效果，即可使用下列計(jì)算公式：

用式（6）和本文作者給出的計(jì)算框圖成功地計(jì)算出哺乳動(dòng)物[13]和原生動(dòng)物[14]的物種分歧時(shí)間，結(jié)果與化石年齡相符很好，明顯好于現(xiàn)有基于分子鐘假說(shuō)的計(jì)算結(jié)果。文獻(xiàn)[15]也指出，典型的蛋白質(zhì)序列分子回復(fù)突變概率為其正向突變率的。需說(shuō)明的是，式（6）得到的絕對(duì)進(jìn)化速率是真實(shí)值的2倍[3]。

2.2 核苷酸序列分子中回復(fù)突變和平行突變的校正方法

核苷酸序列分子中回復(fù)突變和平行突變的校正方法基本上與蛋白質(zhì)序列分子相似。但核苷酸只有4種，即嘌呤核苷酸和嘧啶核苷酸各2種。核苷酸序列分子的位點(diǎn)突變可分為嘌呤核苷酸和嘌呤核苷酸之間，嘧啶核苷酸和嘧啶核苷酸之間的轉(zhuǎn)換突變，以及嘌呤核苷酸和嘧啶核苷酸之間的顛換突變。核苷酸序列分子中轉(zhuǎn)換概率一般大于顛換概率，如線粒體分子中轉(zhuǎn)換與顛換概率比為15∶1或更高[16]，其核苷酸轉(zhuǎn)換突變中出現(xiàn)回復(fù)突變和平行突變的概率與顛換突變中出現(xiàn)回復(fù)突變和平行突變的概率比值為225∶1或更高，前者遠(yuǎn)大于后者,后者的概率可忽略。此時(shí)同源核苷酸序列分子比較時(shí)需校正回復(fù)突變和平行突變的概率之和為其序列差異率的。按照類(lèi)似于蛋白質(zhì)序列分子回復(fù)突變和平行突變的校正方法，可以得到核苷酸序列分子的相應(yīng)校正公式為0.25(e-k1t+e-k2t)=0.5-p。由此得到的核苷酸序列分子絕對(duì)進(jìn)化速率也是真實(shí)值的2倍。

由于蛋白質(zhì)序列分子位點(diǎn)實(shí)際可發(fā)生突變的氨基酸種類(lèi)約為10種，而上述研究體系中核苷酸序列分子位點(diǎn)實(shí)際可發(fā)生突變的核苷酸種類(lèi)為2種。若10種氨基酸和2種核苷酸各自相互之間突變概率相同，根據(jù)文獻(xiàn)[3]，蛋白質(zhì)和核苷酸序列分子自其祖先蛋白質(zhì)或核苷酸序列分子的最大遺傳距離分別約為2.3和0.7，蛋白質(zhì)序列分子的最大遺傳距離大于核苷酸序列分子。

3 結(jié)束語(yǔ)

根據(jù)泊松分布的相加性，我們提出了核苷酸和蛋白質(zhì)序列分子進(jìn)化的“泊松分布矩陣”模型，它能幫助我們進(jìn)一步了解這些序列分子進(jìn)化的機(jī)理，并可用來(lái)推導(dǎo)遺傳距離可加性等計(jì)算公式。我們還提出了分子進(jìn)化研究中回復(fù)突變和平行突變的校正方法，推導(dǎo)出了校正兩種突變后的物種分歧時(shí)間計(jì)算公式，實(shí)際計(jì)算表明它們能得到滿(mǎn)意的結(jié)果。在后續(xù)工作中，我們將進(jìn)一步探討文獻(xiàn)[3]中替代公式能提高計(jì)算效果的原因和計(jì)算時(shí)物種選擇規(guī)則，以期為解決目前依然爭(zhēng)議較多的生物物種系統(tǒng)發(fā)育問(wèn)題提供新的思路和方法。