貝葉斯支端定年法推斷分異時間和演化速率

張 馳

(1 中國科學院古脊椎動物與古人類研究所，中國科學院脊椎動物演化與人類起源重點實驗室 北京 100044)

(2 中國科學院生物演化與環境卓越創新中心 北京 100044)

1 簡介

推斷類群的系統關系和分異時間是分支系統學分析的基礎。如何合理利用化石的形態和地質年代的數據來完成此類推斷一直是一個棘手的問題。傳統的推斷過程往往采用分步的策略。首先單獨從形態數據出發，利用簡約法構建類群的系統關系。這個關系只有拓撲結構而沒有時間的信息，而且只代表給定時間和搜索條件下最優的結果(即最簡約的樹或其合意樹)。然后固定這個拓撲結構，從化石年代出發，利用最小枝長法(Laurin,2004)或等距枝長法(Brusatte et al., 2008)確定樹中內部節點的分異時間。最小枝長法是把當前節點的時間往前推一個百萬年作為祖先節點的時間；而等距枝長法則是把祖先節點到后代節點的時間的中點最為當前節點的時間。至于演化速率，可以通過祖先狀態重建和每個樹枝上特征變化的次數，結合上一步推斷的時間來估計。這種分步計算的方式比較直觀，因此被用于實際數據分析(Wang and Lloyd, 2016)。然而，該策略存在諸多缺陷。首先，它每一步都忽略了推斷中的不確定性，包括樹的拓撲結構、分異時間和祖先特征狀態；其次，每一步都只利用了一部分數據信息，如建樹時只用到形態特征，定年時只用到化石年代；再次，定年的方式很主觀，對化石數量的增減很敏感，且不適用于現生類群；最后，整個過程缺乏一個嚴謹的統計學框架，無法對不同的模型假設進行檢驗。

近些年開發的貝葉斯支端定年法(Bayesian tip dating) (Ronquist et al., 2012;Gavryushkina et al., 2014; Zhang et al., 2016)很好地克服了上述問題。貝葉斯支端定年法把化石形態和年代數據整合在一次完整的計算過程中，能夠盡可能地利用數據信息，同時考慮了樹的拓撲結構、分異時間、演化速率以及化石年代的不確定性。該方法通過統計模型來描述特征的演化、類群的生滅以及化石的采樣等過程，并借助相對成熟的貝葉斯統計框架和計算方法來進行參數估計和模型選擇。但該方法相對來說比較復雜，需要較多的統計學知識，對古生物學家來說往往難于理解和上手，而且鮮有系統闡述支端定年法計算過程和參數意義的文獻(Gavryushkina and Zhang, 2020)。本文逐層剖析支端定年法的計算過程，解釋其中用到的重要模型和參數意義，旨在一定程度上為古生物學家分析實際數據提供參考。

本文首先介紹描述時間樹的石化生滅過程(fossilized birth-death process)模型(Stadler,2010), 然后介紹描述特征演化速率的寬松形態鐘(relaxed clock)模型，接著介紹描述特征狀態變化的Mk模型(Lewis, 2001), 再通過貝葉斯公式把上述模型聯系起來，最后介紹估計參數后驗分布的馬氏鏈蒙特卡羅(Markov chain Monte Carlo, MCMC)算法。附錄提供了計算中生代鳥類數據(Zhang and Wang, 2019)的MrBayes命令。

2 分異時間

時間樹(timetree)代表類群的系統關系和分異時間，它的概率分布可以通過石化生滅過程(Stadler, 2010)來給出。該過程描述了從這些類群的最近共同祖先(樹根)開始，分異、滅絕、采集化石和采樣現生類群這一系列事件的發生，并對應于一棵完整樹(圖1A)。但是實際數據分析中無法推斷這個完整樹，只能推斷和樣本相關的部分，即樣本樹(圖1B)。記每個樹枝的分異速率(或叫成種率)為λ, 滅絕速率為μ, 沿每個樹枝的化石采樣速率為ψ, 現生類群的采樣概率(或采樣比例)為ρ, 通過建立并求解一系列常微分方程，可以得到給定λ,μ,ψ,ρ時，樣本時間樹T= {τ, t}的概率分布，記為P(T|λ,μ,ψ,ρ), 其中τ代表拓撲結構，t代表以百萬年為單位的分異時間。

圖1 石化生滅過程可能得到的時間樹Fig. 1 Example timetree generated from the fossilized birth-death process

在MrBayes軟件中，該生滅過程以樹根時間t1為起始條件(圖2), 計算時要指定t1的先驗分布。通常這個先驗比較寬泛(最大范圍從0到無窮), 不過一般能從研究的類群預估一個更精確的范圍，例如，其下界不會早于最古老化石的年代。化石的時間可以固定為具體的數值(百萬年前), 也可以用一個均勻分布給出時間的上下界。在計算時還需要提供現生類群大致的采樣比例(ρ)。現生類群可以有兩種采樣策略，一種是均勻隨機采樣(random), 另一種是多樣化采樣(diversity) (Zhang et al., 2016), 可以根據實際數據的情況自行選擇，后者可能更符合高階元類群的采樣模式(如每個科只取一個代表的屬或者每個屬取一個代表物種)。對于分異、滅絕和化石采樣速率，程序為了設置先驗方便，重新參數化為d=λ-μ,v=μ/λ,s=ψ/ (μ+ψ)。d的默認先驗為指數分布(范圍從0到無窮),v和s的默認先驗為均勻分布(范圍從0到1,更一般地為貝塔分布)。這樣，時間樹包括分異時間等參數的先驗分布基本就確定了。

圖2 用于各個概率分布公式中參數的示例樹上化石F1和F2的特征狀態為0,現生類群S1和S2的特征狀態為1,內部節點的特征狀態用x0, x1, x2表示化石的年代分別為t3 = 100 Ma, t5 = 50 Ma,樹根的時間為t1, 其他分異時間為t2和t4記t = {t1, t2, t3, t4, t5}。各個枝上的特征演化速率為r = {r1, r2, r3, r4, r5}Fig. 2 Example parameters and symbols used in the probability distributions The character states for fossils F1 and F2 are 0, for extant taxa S1 and S2 are 1, for internal nodes are x0, x1, x2. The ages of fossils F1 and F2 are t3 = 100 Ma and t5 = 50 Ma The root age is t1 and the remaining divergence times are t2 and t4. Denote t = {t1,t2, t3, t4, t5}. The evolutionary rates on the branches are r = {r1, r2, r3, r4, r5}

需特別提到的是，有些數據只包含了化石而沒有現生類群。這時一般假設生滅過程在未到達現今時間點所有類群就都滅絕了。因此現生類群不論采用何種采樣策略和比例(程序默認為1.0), 都沒有樣本被采集到。不過由于MrBayes軟件的限制，最年輕的樣本總被顯示為現生類群(時間為0), 因此需要把整棵樹的時間軸進行相應的平移，或者說分異時間要加上年代最近的化石的年代。這只是結果顯示的問題而不是程序計算的錯誤(后續版本將修復這個問題)。

除了石化生滅過程之外，MrBayes還為時間樹提供了一個均勻分布的先驗(Ronquist et al., 2012)。它沒有生滅和采樣這些參數，只依賴于樹根時間t1, 因此只需設置t1的先驗分布。均勻分布時常被認為是沒有信息的先驗，但它實際上往往帶有很強的信息，在定年這個問題上，會影響對分異時間的估計。而石化生滅過程看似參數很多，但實際上其設置可以很靈活。例如，分異、滅絕和化石采樣速率都可以隨時間變化，在不同的時間段內各自獨立(Gavryushkina et al., 2014; Zhang et al., 2016)。這可能更符合實際生物學過程，同時還能推測凈成種速率和化石采樣速率隨時間的變化。

3 演化速率

形態特征的演化速率是指每百萬年每個特征期望的變化次數。對于給定的一段時間，演化速率越快，則特征最終期望的變化次數越多。一般給每個樹枝一個演化速率參數，記為r (圖2)。時鐘模型應用于形態數據時被稱為形態鐘模型，類比于用在分子數據時的分子鐘模型。嚴格鐘(strict clock)模型假設演化速率在各個樹枝上都相同，這通常不適用于形態數據，實際分析時往往需使用寬松鐘(relaxed clock)模型。寬松鐘模型可以分為兩類，一類為獨立速率，另一類為自相關(autocorrelated)速率，區別在于r的概率分布P(r)不同。

獨立速率模型假設枝上的演化速率彼此獨立，它們都服從均值相同的某個概率分布。常用的概率分布包括伽馬分布(Lepage et al., 2007)和對數正態分布(Drummond et al.,2006)。分布的均值也被稱為基準速率(base rate), 反映平均的演化速率大小。分布的方差則反映演化速率在樹枝之間變化的劇烈程度：方差較小時各個速率相差不大，這意味著演化速率在整棵樹上沒有明顯的差異；而方差越大，不同樹枝上速率的差異越明顯。

自相關速率模型假設后代樹枝上的演化速率依賴于臨近祖先那枝上的速率(例如r2和r5都依賴r1,r3和r4都依賴r2)。當前樹枝上的速率一般假設服從對數正態分布(Kishino et al., 2001; Thorne and Kishino, 2002), 其均值為臨近祖先節點的速率。同理，分布的方差也反映演化速率在樹枝之間變化的劇烈程度。

這兩類速率模型往往對分異時間的估計也有影響，這主要是因為自相關速率模型會傾向于速率的變化是漸進的，而獨立速率模型沒有這種限制，會更適應臨近樹枝間速率變化比較劇烈的情況。對化石形態數據來說，獨立速率模型可能更適用。

默認的情況下，形態數據矩陣中所有特征都共享每個枝上的演化速率，因此，這個速率代表的是所有特征的平均情況。如果需要考慮不同特征演化速率的異質性，就需要對特征進行分區。一般可以按照不同特征類型或不同身體部位或功能來分，每個分區內的特征共享一組演化速率，而分區之間特征演化速率的模式是獨立的，這樣就可以推斷不同部位或功能相關特征隨時間會發生怎樣的變化(Lee, 2016; Zhang and Wang,2019)。需注意的是，分區越多，每個分區內的特征數量就越少，因此能夠估計演化速率參數的信息就越少，會造成方差很大甚至參數個數超過特征數量導致無法進行參數估計。因此，要在考慮演化速率異質性和分區數量之間做一個權衡。

4 特征狀態變化

有了分異時間和演化速率，就可以計算在給定時間段t和演化速率r的情況下，形態特征從一個狀態變為另一個狀態的概率(稱為轉移概率)。這個概率由Mk模型(Lewis,2001)給出。Mk模型是描述特征狀態變化最簡單的模型，它假設狀態之間轉換的速率是相等的。這里只以兩個狀態的特征為例，用P00(r,t)表示狀態0保持不變的概率，P01(r,t)表示從0變到1的概率，P10(r,t)表示從1變到0的概率，P11(r,t)表示狀態1保持不變的概率，則P00(r,t) =P11(r,t) = 1/2 + 1/2 × exp(-2rt),P01(r,t) =P10(r,t) = 1/2 - 1/2 × exp(-2rt)。從公式中可以發現，時間t和速率r總是以乘積的形式出現，因此，在沒有化石年代信息時，兩者是不可識別的。換句話說，單單依靠形態數據來建樹，樹的枝長為時間和速率的乘積，即距離，以每個特征期望的變化次數為單位。只有同時利用化石形態和年代數據才能將分異時間和演化速率單獨估計出來。

以圖2為例，化石F1和F2的特征狀態為0, 現生類群S1和S2的特征狀態為1, 內部節點的特征狀態未知用x0,x1,x2表示。F1和F2的時間分別為100和50 Ma。那么根據Mk模型，給定時間樹T= {τ, t}和速率r時，特征狀態列0011的概率為P(0011|T, r) = ∑x0∑x1∑x2Px0x1(t1-t2,r1)Px1x2(t2-t4,r2)Px21(t4,r3)Px20(t4-t5,r4)Px11(t2,r5)Px00(t1-t3,r6)。

其中西格瑪符號代表對特征在內部節點所有可能狀態的求和。由于形態特征矩陣往往只包含可變的特征，因此這個概率還要除以所有可變狀態的概率，即P(0011 |T, r)/[1 -P(0000 |T, r) -P(1111 |T, r)]。帶有這一校正的Mk模型稱為Mkv模型(Lewis, 2001)。

假設形態矩陣中的特征都彼此獨立，那么就可以計算每一列特征在樹上的概率，再把這些概率乘起來。這個概率被稱為似然函數，表示為P(D|T, r), 其中D代表形態特征矩陣數據。

5 貝葉斯公式

在統計推斷時，參數都是未知的隨機變量，需要根據數據來估計它們的分布，即計算P(T, r,θ|D), 該分布稱為后驗分布，其中T= {τ, t}為時間樹，r為演化速率，θ代表其它參數(包括λ,μ,ψ等)。根據分層貝葉斯公式，可得P(T, r,θ|D) =P(D|T, r)P(r)P(T|θ)P(θ)/P(D)。等號右側分子中，第一項似然函數在第四節給出，第二項演化速率的先驗分布在第三節給出，第三項和第四項為時間樹及其參數的先驗分布在第二節給出。這樣公式分子中各項都可以計算了。分母P(D)是特征數據的概率，這需要計算對所有參數的多重積分，實際上基本無法給出解析表達式，只能通過數值算法進行近似。所以貝葉斯計算在絕大多數情況下會使用馬氏鏈蒙特卡羅算法。

6 馬氏鏈蒙特卡羅算法

馬氏鏈蒙特卡羅(MCMC)算法通過構造馬爾可夫鏈，使其平穩分布為要估計的后驗分布。這里為了簡化，只以一維參數的情形為示例(圖3)。實際分析中，參數一般是多維的(如τ, t, r,θ), 不過算法的思想是類似的。

MCMC采用的Metropolis-Hastings算法(Metropolis et al., 1953; Hastings, 1970)可分為如下幾步：

(1) 為θ設定任意初始值；

(2) 基于θ當前的值，建議一個新的值θ’, 例如θ’ ～ uniform(θ- w/2,θ+ w/2);

(3) 如果π(θ’) ＞π(θ), 就接受θ’; 否則，以概率α=π(θ’) /π(θ) 接受θ’;

(4) 如果θ’被接受，就更新θ=θ’; 否則，保持θ不變；

(5) 記錄θ的值，回到步驟(2)。

注意到，在計算概率α的時候，會計算后驗分布的比，這樣后驗分布分母的部分就約掉了，只剩分子部分的比。也就是說，只要能夠把分子部分寫出解析表達式，MCMC算法就可以使用來估計參數的后驗分布了。

計算結束后，就收集到一些θ的樣本。由于參數的初始狀態往往比較差，MCMC需要經過很多代才收斂到后驗概率密度比較高的地方，因此在估計后驗分布的時候會舍棄初始的一些樣本(burn-in), 只用MCMC鏈收斂后記錄的樣本來估計后驗分布。MrBayes默認舍棄前25%的樣本。同時，MCMC鏈還要迭代足夠多次，以保證有足夠多的有效樣本來估計參數的后驗分布。一般需要有效樣本大小(ESS)大于100。

實際運算中，最好獨立運行至少兩次MCMC, 以確保兩次的結果是一致的。有時鏈長不夠，或者不同的運算卡在不同的后驗分布區域，都會導致估計的結果不一致。這時調整MCMC的設置或者改善模型都可能幫助MCMC算法發揮更好的效能。使用Metropolis-coupled MCMC也是有效跨越多峰分布的手段(Lakner et al., 2008)。該算法同時運行多條MCMC鏈，一條為冷鏈(cold chain), 其余為熱鏈(hot chains), 熱鏈和冷鏈之間可以相互交換。當然MCMC樣本只從冷鏈中采集，熱鏈只是用來幫助跨越多峰的。MrBayes默認同時獨立運行兩次，每次運算使用四條鏈，其中一條為冷鏈，其余三條為熱鏈。

圖3 馬氏鏈蒙特卡羅算法在一維情形的示例參數θ的后驗分布π(θ)為估計的目標(曲線)Fig. 3 Illustration of the MCMC algorithm for one dimensional parameter The posterior distribution π(θ) of parameter θ is the target distribution to be estimated (curve)

7 討論

本文從貝葉斯統計計算的角度分層剖析了支端定年法的原理和計算過程。貝葉斯后驗分布包含先驗分布和似然函數，其中先驗分布的兩個重要組成部分：分異時間和演化速率模型，在定年分析中尤為關鍵，是影響定年準確性的主要因素。

石化生滅過程作為描述類群分異、滅絕和采樣的隨機過程，具有較大的靈活性。不過該模型還有待完善之處。現生類群的采樣方式可以是隨機的或多樣化的，這兩種情形可能都是極端，真實的采樣模式可能介于兩者之間，或者有的支系是隨機采樣，有的支系是多樣化采樣，但目前還沒有模型能夠支持這種情況。在更好的模型被開發出來之前，可能只能調整數據來盡量符合其中一種采樣策略。這種處理方式在現生類群比較少或根本沒有時一般問題不大。另外，分異速率、滅絕速率和化石采樣速率可以按分段的方式隨時間變化，不過在同一時間段內，所有樹枝都共享同一速率。對于不同支系明顯受到的選擇壓力不同或化石保存的完整性明顯不同等情況，按支系分段而非時間分段的模型(Barido-Sottani et al., 2020)可能更合適，雖然這一類模型本身也有其他限制。不論怎樣，石化生滅過程只是作為時間的先驗分布，當數據量比較大時(包括化石在樹上的分布程度和數量以及形態特征的數量和完整度), 數據在推斷中會起主導作用而先驗的影響減少。但是實際情況往往比較復雜，數據量也不盡如人意，這時考察不同先驗的影響就尤為重要。

演化速率的先驗，即形態鐘模型，也會和時間相互作用，從而影響對最終分異時間的估計。這種影響在化石較少或化石在樹上分布很不均時尤為明顯。這主要是因為化石形態數據只提供了距離的信息(每個特征期望的變化次數), 其為時間和速率的乘積(見第四節)。當缺少化石時，就沒辦法準確提供時間的信息，那么對于同樣的距離，可以是很長時間速率很慢，也可以是很短時間但速率很快，具體是怎樣就只能取決于時間和演化速率的先驗了。對于某些支系時間估計得明顯偏大或偏小，但又沒有化石來校正的情況，可以通過添加內部節點的校正分布來得到更合理的時間估計(O’Reilly and Donoghue, 2016)。在完全沒有化石只有現生類群時，節點定年法(另一類型定年方法)(Yang and Rannala, 2006)就是通過使用內部節點的校正分布來完成的。

描述形態特征狀態變化的模型也有很大的改進空間，其中涉及更多的建模和隨機模擬等工作，不是本研究的重點，這里只簡單討論一下Mk模型對定年可能的影響。Mk模型假設特征各個狀態間轉變的速率都是相等的。這種轉變有無序和有序之分(只對三個及以上狀態的特征)。無序是指特征可以直接從一個狀態變為任意其他狀態(如從0直接變為3), 而有序是指特征只能在臨近狀態間直接變化(如從0到1, 從1到2, 再從2到3)。顯然，有序需要更多次變化(也就是更長的距離)才能從當前狀態變為不相鄰的狀態，因此對有序的特征假設無序的變化會低估距離。更復雜的情況是，各個狀態間轉變的速率未必相等甚至相差很多，極端情況像Dollo特征甚至是不可逆的。這時使用Mk模型也會造成距離估計的偏差。在計算時還假設不同特征之間都是獨立的，如果有些特征有較強的相關性，則會導致演化距離的高估。前面提到，距離是時間和速率的乘積，因此在化石(時間信息)很豐富的理想情況下，距離的偏差會主要體現到演化速率上而對分異時間影響較小。但是分析實際數據時會更復雜一些，要具體問題具體分析。相關的工作還較少(Klopfstein et al., 2019), 需要更多后續研究來更詳細地考察。

最后提到貝葉斯計算使用的MCMC算法。該算法的策略與簡約法和似然法有明顯不同。簡約法尋找的是簡約樹長最小的樹，似然法尋求的是似然值最大時參數的估計(即最大似然樹)。因此在設計樹的搜索方法時，只要盡可能快速地找到最優的樹就可以了。MCMC算法是為了估計參數的后驗分布，這不僅僅是一個點，而是參數的空間。因此MCMC算法的效率涉及收斂(convergence)和混合(mixing)兩部分。收斂是指MCMC達到分布概率密度高的區域，混合是指MCMC能夠按概率分布進行取樣。提高收斂速度相對容易，可以通過如簡約樹長向導的方式來快速找到概率大的樹(Zhang et al.,2020)。提高混合則更困難，需要設計更好的建議(proposal)方法，這是貝葉斯計算的重點也是難點。

總之，貝葉斯支端定年法作為整合的分析方法，能夠結合化石形態和年代數據以及現生類群的形態和分子數據來推斷類群的系統關系，分異時間和演化速率，同時考慮了樹的拓撲結構、分異時間、演化速率以及化石年代的不確定性。但該方法仍處于發展初期，模型和算法的諸多方面還亟待完善，因此還有很多工作需要做。

Supplementary material can be found on the website of Vertebrata PalAsiatica (http://www.vertpala.ac.cn/EN/2096-9899/home.shtml) in Vol. 59, Issue 4.