airbugs

一個特徵（trait）的不一樣的性狀（states）對一個生物系群分化及絕滅的影響有多大呢？ 例如：多形性（顏色變化多端）的鳥是不是比顏色單調的鳥有較高的分化速度？ 類似的問題一直是生物學家很有興趣研究的方向，因為它們都關係到生態及演化學最基本的大問題“為什麼世界上有這麼多種生物？為什麼有些生物的類型比其他類型的生物在地球上成功很多（例如族群量比較大或是分化成比較多的物種）？” 自從2007年的BiSSE model被發表之後，這個Trait dependent speciation & extinction model就一直被擴充。由本來只能fit兩個不一樣的states到多個states，甚至到continuous trait。最大的成就就是這個model可以用來比較不一樣的trait states對一個生物系群分化及滅絕速率的影響。 Ng & Smith (2014)這一篇文章有對這類型分析的完整回顧。這個model最大的假設就是你有完整的採樣（missing data很少並且是隨機的：不會有某個系群有特別完整個sampling然後其它的系群很多missing data），另一個前提假設就是整個系群的分化過程都是因為你所選取的Trait造成的！ 大家都知道這兩個假設（尤其是第二個）真的是很難達到又極度天真，但是不這樣假設的話就沒東西可以算了。所以在分析這些資料之前，先確定你有一個採樣非常完整的親緣關係樹，並且你相信你要討論的不同性狀對這個系群的演化過程影響甚巨，那就可以繼續走下去。 這個類型的model有三個主要的參數：lambda，mu，及q（Ng & Smith 2014的圖）。 0（白）與1（黑）代表一個trait的兩個不同states。lambda是擁有不一樣state的speciation rate（e.g., 多常看到0生0），mu是擁有不一樣state的extinction rate（e.g., 如果一個物種的trait state是0，它的extinction event多常發生？），然後q是不同states之間轉換transition rate（e.g., 多常看到原本是state 0的生物轉變成state 1）。其它的models都包含一樣的參數，只是變得比較複雜。 R有一個很方便的package可以做BiSSE類型的分析叫做diversitree。 以下就用我的研究類群Dynastes甲蟲的親緣關係樹（species tree output from *BEAST）分析speciation & extinction rates是不是與兩個很有趣的biological/ecological traits有關：（1）與高海拔地海拔居住的關係，（2）與居住在北美或南美洲的關係。 library('diversitree'); #先把需要的library load到R dynastes1 <- read.nexus('dynastes.tre'); #用read.nexus這個function把*BEAST的MCC species tree output讀到R裡面去，然後給他一個名字叫做dynastes1 下面這一步是這個dataset特別需要的，其它dataset不一定要： dynastes <- drop.tip(dynastes1, 'Dht'); #我自己的研究發現Dhb與Dht也許只是不同的地理族群關係，彼此之間的supports for evolutionary independence from iBPP並不高。所以需要把其中一個去掉，去掉之後的新tree叫做dynastes。 可以把兩個親緣樹都plot出來比較看看： par(mfrow = c(1,2)); #一次畫兩個並排的圖 plot(dynastes1); #左圖 plot(dynastes); #右圖 左邊比右邊多了一個taxon Dht： 右邊的tree就是接下來分析需要用的tree。 ecodata <- read.csv('ecodata.csv'); #把相對應每個taxon的ecodata叫進去R 我是存成有header的csv類型格式（下面是用wrangler開來看：0表示高海拔，1表示低海拔）： eco <- ecodata$state;#這是要分析的0跟1的data names(eco) <- ecodata$species;#給每個0跟1一個對應到tree上面的taxon的名字 接下來就是solve一堆likelihood functions： lik <- make.bisse(dynastes, eco);#使用BiSSE model p <- starting.point.bisse(dynastes); fit.full <- find.mle(lik, p); #full model，所有的參數都要估計 然後在設定一些constrained models（兩個例子） lik.1 <- constrain(lik, lambda1 ~ lambda0);#constrain高海拔跟低海拔對speciation rate的影響相同 lik.q <- constrain(lik, q10 ~ q01);#constrain兩個states有一樣的transition rate fit.1 <- find.mle(lik.1, p[argnames(lik.1)]);#找到maximal likelihood的點估計 fit.q <- find.mle(lik.q, p[argnames(lik.q)]); 然後做likelihood ratio tests看不同models（full model vs. constrained）有沒有統計上的差異： anova(fit.full, eq.l = fit.1, eq.q = fit.q); 當把假設高低海拔居住有相同的speciation rate時，這個model快要顯著的比full parameterized model差。但是如果是constrained高低海拔住之間的transition rate相同，對model fit的影響不大。 除了likelihood的分析之外，我們也可以用MCMC來估計model裡面的參數值。這邊我用Full model來做MCMC search： prior <- make.prior.exponential(1/(2*(p[1]-p[3])));#如果你不喜歡exponential，設其它的也可以 set.seed(1); #不一定要，但是可以看可重複性（如果你有興趣的話） tmp <- mcmc(lik, fit.full$par, nsteps = 100, prior=prior,lower=0,w=rep(1,6), print.every=0); #先跑一個short run來粗估一下參數的範圍 w <- diff(sapply(tmp[2:7], range));#把參數範圍算出來 mcmc2 <- mcmc(lik, fit.full$par, nsteps = 100000, prior = prior, w=w, print.every = 100);#真正的long run #總共跑10的5次方generations（Mac power book大概跑2個小時吧？[出去吃飯沒注意到]） 然後就是做圖看MCMC估出來的參數值： col1 <- c('gray','black'); par(mfrow = c(3,1), mar = c(3,4,1,1)); profiles.plot(mcmc2[2:3], col.line=col1, xlab='', ylab='');#這兩個是lambda profiles.plot(mcmc2[4:5], col.line=col1, xlab='', ylab='');#mu profiles.plot(mcmc2[6:7], col.line=col1, xlab='', ylab='');#q title(xlab = 'rate', ylab = 'posterior prabability density', outer = TRUE, line = -1); 上圖是speciation rates，中間extinction rates，下面transition rates。黑色是低海拔（1），灰色是高海拔（0）。 估計出來的結果居住在不同海拔高度有顯著不同的speciation rates（低海拔顯著高過高海拔）。除此之外，海拔高度對extinction及transition rates沒有顯著的影響！ 我們可以同樣的比較居住在南美洲及北美洲有沒有顯著不一樣的分化模式。因為地理的分佈可以是廣布的（可以同時被發現在北美洲及南美洲），所以這個分析需要用到改良的GeoSSE model。 大部分分析用的R code都與BiSSE model相同，只有以下兩點需要注意： （一）geodata的輸入需要注意到0代表”廣布“，然後1跟2分別代表不同地理區（在我的例子就是南美或北美）。 （二）把make.bisse變成make.geosse，因為是要用GeoSSE model就要設定geosse的likelihood function 一樣從頭到尾計算完畢就可以得到下面這個圖： 基本上估計的參數全部都彼此overlap。所以住在南美或是北美對不同evolutionary independent lineage的speciation，extinction，及transition rates都沒有明顯的影響！ 最後，當然這些model有很多限制也有很多天真的前提假設，估出來的結果也會受到很多不同東西的影響。比如Dan跟Emma就指出如果一個lineage曾經有significant shift in diversification rate（做個簡單的LTT plot或是用MEDUSA或BAMM之類的models都可以測試看看），那這些BiSSE family的Trait dependent speciation & extinction model就很有可能會估計錯誤參數。所以在分析data之前，還是要先了解一下自己的採樣過程，問的問題，data的組成，然後再看看適不適合這類型的分析。 如果一切都不是問題，你也相信你的dataset超棒，那就分析看看你研究的系群是不是因為某些geological events或是key innovations造成在不同類群之間有不同程度的分化！