airbugs

現在生態演化學的大家越來越多使用基因體階級的資料類型了，也越來越多人用這些大量的資料來研究族群分化及物種分化層級的問題。這些都是非常棒的事及發展方向！

族群基因體學研究有一個非常常被使用的軟體Structure（或是ADMIXTURE），它可以根據你輸入的遺傳資料幫你計算出最可能能夠解釋你收集到的資料中的遺傳變異需要有幾個不同遺傳來源的族群（k）。然後再幫你計算每一個你輸入的個體（individual）的遺傳資料有多少百分比是來自各個不一樣的遺傳來源。

這是一個很方便很好用的軟體，它的結果圖像化（visualization）也非常直觀一目瞭然：比如說下面這一個圖很明顯的說明你輸入的遺傳資料最合適的解釋目前觀察到的遺傳變異的遺傳群（Wright-Fisher population; k）是2: 紅色及藍色。蠻多個體（individual）的遺傳組成可以純粹由其中一個遺傳群解釋（要不是藍色就是紅色），但是中間有一些個體的遺傳組成同時包含來自兩個遺傳群的來源。

像上圖這樣的結果，有時候會讓人想到中間黃框框內的個體是不是藍色及紅色兩個族群的雜交後代呢？這樣的推論當然是可能的，而且也是非常發生的推論，各大研討會及期刊發表上面都可以看得到。但是我這裡想說的是雜交並不是唯一一個可以解釋目前看到上圖這樣的結果（pattern）的假說，有些特定的族群分化歷史也可以導致一樣的Structure分析結果。

 

在Lawson 2018年寫的文章介紹如何看待STRUCTURE分析結果的文章，其Fig. 2就有對這種結果不要太“過多解釋(overinterpret)”的介紹。

如上圖所示，b panel 的三個STRUCTURE結果是極度相似的，看到這樣的STRUCTURE很容易讓研究人員就會做a panel最左邊的歷史推論（recent admixture：P2族群是雜交後代），近年來也很多研究會討論panel a中間那個歷史模型（過去的遺傳交流，只是其中一個種源滅絕了：Ghost admixture：P2族群是雜交後代）出現的可能性。以上這兩個解釋或多或少都有用歷史上的遺傳交流（或雜交）來解釋目前觀察到的遺傳多樣性。

然而，同樣用電腦模擬，如果中間這一個P2族群在很近代有一小群人出走建立了遺傳多樣性很小的P1族群，那STRUCTURE分析一樣會算出一樣的結果！雖然Recent Bottleneck這一個歷史模型完全沒有加入任何遺傳交流的歷史參數（parameter），STRUCTURE因為估出來最佳的k依然是3，而因為P1是近代出現遺傳多樣性很小的族群（佔了一個k），唯一的解方就是把P2當成是有很多遺傳來源的族群！

要如何分辨你的資料到底是哪一種歷史造成的結果當然還是有很多辦法可以解決，如果有興趣可以好好的看看Lawson 2018這一篇文章（雖然已經有很多citations，但我還是覺得它得到的關注不夠多）。這裡主要想分享的就是大家data量越來越大，分析也可以做很多，看到patterns後會有很多很可能有趣大家也都這樣推測的歷史原因（都很好），但是不要忘了其他也可能的alternative hypotheses！

Lawson 2018 的文章連結：

https://www.nature.com/articles/s41467-018-05257-7

 

好久沒更新了，來寫一個大家好像不太在意的小地方Spurious Corrlation：

有時候在做生態學研究的時候，因為每個生物個體會有很多本身的差異，所以在做分組比較的時候很多測量值會要標準化。比如鍬形蟲的大顎長度會與體長有關，越大隻的雄蟲越長。鍬形蟲的翅膀長度也會與體長有關，越大隻翅鞘就越長。那如果研究大顎長度與翅鞘長度是不是有相關，能不能把兩個測量值（大顎長度及翅鞘長度）對體長做標準化（變成ratio）在來看兩個比例的關係呢？

要做當然是可以，但是還沒分析就可以知道這兩個標準化後的ratio一定會變成有直線相關性，因為我們把兩個不相關的測量值用了一樣的denominator：Spurious correlation就出現了！

 

簡單的舉個例子，假設有一個生物身高是如下隨機分佈：

sim_height <- rnorm(100, mean = 170, sd = 10)

此生物的手長也是隨機分佈如下：

sim_armL <- rnorm(100, mean = 60, sd = 5)

這兩個測量值是不會有直線相關性的！如下：

> fit1<-lm(sim_armL~sim_height)
> summary(fit1)

P value約0.6

 

如果假設這個生物的體重也是隨機分佈如下：

sim_weight <- rnorm(100, mean = 90, sd = 3)

這時候我們這個生物共有三個隨機分佈的測量值（如下）：

 

接著我們把體長跟手長先對於體重做標準化（純粹亂做）：

> new_height <- sim_height/sim_weight
> new_armL <- sim_armL/sim_weight

然後你就會發現標準化後的ratio們彼此就有顯著相關了！：

> fit3 <- lm(new_armL~new_height)
> summary(fit3)

P < 0.05！是不是很妙！

 

所以大家做研究標準化很重要是沒錯，但是要謹慎選用方法，不然就會找到假的相關性，做出錯誤的詮釋。研究有時候會聽到建議說結果沒有顯著也許是因為沒有好好標準化...但是大家也要小心只有標準化後才看到的顯著結果！

 

寫這一篇文章說的內容其實在生物學研究相關的統計方法已經被討論很久了（但是有興趣或是會特別注意的人好像不多），有興趣可以看按下面這一篇文章：

https://www.jstor.org/stable/2983064

因為我做的研究與物種分化及如何客觀的判別什麼樣的東西可以被稱為不同“種”，不同種與不同“族群”的生物學差異在哪，很多時候都會遇到有興趣的朋友或是同學問說為什麼不試試看雜交？或者是看到不同物種被定義的時候（例如我研究的不同長戟大兜蟲的“種”），常常能得到的反應就是『可是這些不同長戟可以交配生下沒有問題的下一代，為什麼是不同種呢？』。

但是生物種與雜交不孕（或下一代生育率下降）有直接的關係嗎？為什麼很多人（包含生物老師，生物學家及大眾）想到不同生物的種類總是會聯想到雜交不孕呢？

Mayr在定義生物種的觀念時候僅說了生物種是在自然的狀況下包含一群可以自由繁殖下一代的個體或族群們。不同物種之間因為“生殖隔離”在自然的情況下不會相互繁殖。但就因為這個“生殖隔離”一詞，一些當時不同意Mayr觀點的學者就嘲弄地說出『Mayr的生物種根本上就是雜交不孕的隔離觀念』，也說出了經典的『老虎跟獅子可以生下一代，難道他們是同一種嗎？』的句子。這兩個原本是來嘲弄生物種觀念的句子卻成為後來多數人對生物種觀念的認知，真的是始料未及的成果。

生物多樣性及螞蟻研究很有名的學者E. O. Wilson雖然在1953年就寫過文章抨擊這些攻擊Mayr生物種觀念的人根本就是不暸解族群的概念，也審慎的解釋了為什麼老虎跟獅子是有效的生物種：因為它們在印度等地有共域的機會，彼此有不同的生物學，“自然狀況”下有完整的“生殖隔離”不會產生遺傳混雜的下一代。然而大家就好像沒聽懂，繼續把生物種的觀念跟雜交不孕連想再一起。

這應該是近代生物學上被錯誤訊息影響最大的一個觀念了吧，而且最後假訊息竟然成為大多數人認知中的事實。對Wilson & Brown Jr. 1953文章有興趣的可以在下面連結中找到全文：

https://academic.oup.com/sysbio/article/2/3/97/1674750

至於”生殖隔離“是什麼意思，到現在很多生物學家還是有不一樣的想法及定義（但可以確定的就是生殖隔離並不是雜交不孕的同義詞）。什麼樣的狀況下才能稱作“完整”的生殖隔離就更是很難有共識了。最近Phil Trans Royal Soc有一期專刊討論生殖隔離被完成的期間可能會發生的事目前的暸解，一些相關研究例子及未來研究方向及展望，有興趣的可以從下面這個連結中的專刊導讀開始看起：

https://royalsocietypublishing.org/doi/10.1098/rstb.2019.0528

 

＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃

最後想說的是能把東西養一起嘗試不同種能不能產生下一代是不錯的研究精神，但是這種研究討論的問題是檢測兩種是否“雜交不孕”，而並非兩種是不是不同“生物種”，雜交不孕與是否為不同生物種是沒有一定的相關性的。雜交不孕能探討的問題是兩個個體的基因型產生的下一代在細胞或發育層次上互不相容，此問題本身是一個生物學上很有趣的課題。

很多不同生物種都可以雜交產生健康有繁殖力的下一代：像各種鸚鵡及水族產業中各式各樣雜交品系的魚。

同樣的生物種的不同個體也可以產生不相容的下一代：像人類鐮刀型血球異型和子父母可能會產生有病變的存活率低同型和子下一代（相容的上一代產生不相容的下一代：難道你要說這兩種會有雜交不孕下一代的人類基因型是不同“生物種”嗎？）。

＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃＃

最近常常逛我們博物館的館藏標本，
突然發現Field Museum有好多大王花金龜的標本！！！

最近抽空照了相，但是只挑有"TYPUS"的箱子拍照，
以後也許會有人再重新整理Goliathus的分類學，
或甚至是研究這個屬的分化和體色適應！！

留下這一些照片，希望大家都查得到。



其他更多的照片在相簿裡面。

當兩種親緣關係很近的物種分佈在同一個地方的時候，
這兩種物種是如何得以共存呢？

如果他們的食性都差不多，居住偏好及繁殖習慣都一樣，
這兩個物種要不是彼此競爭排除（有一種會絕種），
不然就會是混回去變成一種（不再有兩種）。

除非有生物學上的分化，這兩個物種才有辦法共存。
其中一個有利於物種共存的機制叫做性狀替換：
Character Displacement.
因為競爭的關係，兩物種在共域時演化出不同的性狀。

實際上的例子也不少，例如在南洋大兜的類群中：
高卡薩斯跟阿特拉斯在東南亞某些地區分布是重疊的。
重疊分布的地區這兩種有很不一樣的生殖器特徵，
雄蟲的大小差異也很大。
但是在不重疊的地區（單一種自己生活的地區），
生殖器變化可以很大，也可能跟其他種的樣式重疊。
雄蟲的大小也可以很類似。

-------------------------------------------------------------------------------


在我研究的姬兜蟲中也可以發現類似的現象：
兩種不同種姬兜共域的時候，雄蟲的頭角及胸角發育型態會非常不同。
例如在西馬來半島：
黑金鋼姬兜（X. s. tanahmalayu）與龍牙姬兜（X. beckeri）共域，
但是黑金鋼姬兜的雄蟲可以長長角型，龍牙姬兜卻是粗角形！

在婆羅洲的例子：
X. wiltrudi與婆羅洲龍牙姬兜（X. pachycera）共域，
X. wiltrudi是長角型，婆羅洲龍牙姬兜是粗角形。

有趣的是親緣關係上馬來黑金剛與婆羅洲龍牙關係比較接近，
龍牙姬兜則是X. wiltrudi的近源種（如下圖顯示）！


*姬兜的手繪圖是Silvestre (2004)的圖。

-------------------------------------------------------------------------------

那這裡的發現可以支持不同姬兜共域時有發生性狀替換的現象嗎？
還是說只是剛好這兩個例子中的蟲子剛好在共域的時候演化出不同型，
跟性狀替換沒有關係，只是演化上的巧合（contingency）？

這個親緣關係樹上總共有10個共域分布的例子，
其中共域分布的兩個物種剛好一個長角型另一種是粗角形的例子有8次。
在10次中發現8次共域不同型，會可能只是演化上的巧合嗎？

*註：兩次例外是菲律賓的菲律賓姬兜（X. phillipinensis）跟鬃毛姬兜（X. pubescens）。
這兩種有族群在Luzon還有Leyte共域。我把它們都算是長角型....

-------------------------------------------------------------------------------

為了回答上面的問題：
在10次中發現8次共域不同型，會可能只是演化上的巧合嗎？
我用電腦模擬來計算根據我的姬兜親緣關係樹：
純粹只考慮演化進程（evolutionary contingency）這狀況下，
10次共域中發生8次不同型的機率多大？

首先我比較不同的性狀（character states）彼此之間轉換的速率，
看看是不是不一樣。我的結果說明轉換速率差不多：
（這裡是用R ape package裡面的ER跟ARD比較）

計算出來Equal rate的rate matrix看起來像這樣：


有了rate matrix（Q）就可以來模擬性狀的演化！

我用了真實的姬兜親緣關係樹做了一千次的模擬，
選了其中隨機6次在下圖中顯示：
因為transition rate並不快，近源種很常是一樣的character state（紅色或黑色）。



結果呢？
1000模擬中，有56次顯示10個共域例子可以產生8個以上的不同型。
..........................
尷尬的結果：P (>=8) = 56/1000 (0.056)。




所以也許有些時候純粹演化上的機緣巧合也能解釋部分的看似性狀替換的現象！
當然我還需要回去看更多文獻，還有再去檢查一下另一個博物館的標本，
因為有些種類到底是粗角形還是長角型我並不清楚，目前結果是用猜的...
也有可能我的tree太小了，沒有足夠的共域例子，所以統計力不夠...

但是至少不是看圖說故事，我們可以簡單的檢測看看性狀替換的假說是否可以用簡單的evolutionary contingency來解釋。

距離我前老闆發表Statistical phylogeography這一篇文章已經經過15年了！親緣地理學的研究也慢慢的從之前根據gene tree來講故事變成用電腦模擬來檢測不一樣的歷史模型對所得到分子證據的解釋度。但是既使現在還是可以看到很多文章就是sequence一些loci，重建幾棵gene tree，然後就根據這些基因樹來推測故事...
這樣的方法用在分化很久的不同種之間的親緣關係分析基本上不會有太大的問題，因為基本上lineage都已經sorting的差不多，不同的lineages也已經有累積自己獨特的變異。
但是用於種內的不同個體，或是探討近緣種姊妹種或複合種群內不同個體的關係時，這些gene trees與不同個體之間的geneological history是沒有一定的相關性的。在近期分化的個體之間如果看到不一樣的基因型基本上都是由祖先多態性（ancestral polymorphism）留下來的結果，所以所重建的基因樹歷史並不會與個體間或族群間的分化事件相關（除非你是研究病毒）。用這樣的資料來討論分化的歷史，只是亂看錯誤的圖然後自己編故事。

現在可以做族群分化+親緣地理model test的電腦程式越來越多了，有些也越來越平易近人。DIYABC就是其中一個很容易上手的程式。他幾乎所有的分析都可以用“按”的！這讓沒有programing 經驗的研究員方便很多。以下我就用一個簡單的例子來玩玩怎麼檢測不同歷史模型對所得到真實分子資料的解釋度--哪個model對data有最好的fit？


在姬兜蟲（genus Xylotrupes）裡面有一個種叫做siamensis。這一個種在過去被分成很多不一樣的種或是亞種（例如socrates，tonkinensis等等），其最主要的原因就是雄性成蟲的頭角和胸角在不同地理族群之間有很不一樣的發育性狀。例如在泰國大型雄蟲會長很長的胸角並且頭角有個突起（siamensis form），在越南及中國南方和中國西南等地的大型雄蟲則沒有這些性狀，但是通常會長很粗的胸角，頭角前方會往回勾（翹角姬兜tonkinensis form）。這些姬兜基本上都是同一種，生殖器完全相同，因此在Rowland 2011年的名錄裡面整個中國西南（不包含西藏）及中南半島（不包含馬來半島）就只有一種姬兜X. siamensis。

那這時候就有一個很有趣的問題了！這兩群長得很不一樣的X. siamensis到底是怎麼演變出來的呢？


根據生態棲位模擬的結果，現今中國西南，印度阿薩姆及中南半島都是X. siamensis適合的棲地：

但是在上一次大冰河期，X. siamensis只有越南北部及泰國西南部保有兩個主要的合適棲地：


那根據過去的地理分佈變化能不能解釋我們現在看到兩群型態很不一樣的地理族群呢？
現在我手上剛好有三個基因的資料，我把不一樣採樣點的個體歸納到四個地理範圍之中：
Pop1：中國東南及越南，Pop2：雲南（型態像翹角姬兜），Pop3：緬甸靠泰國面的山區（長得像泰國姬兜），Pop4：泰國。

我可以大略得馬上想到三個不一樣的歷史模式來解釋現在看到的兩個地理群：
（一）頭胸角有長有短的是祖先型（來自泰國南部的冰河期避難所），冰河期結束後往東北擴散，然後新族群失去長角特徵。
（二）短角無頭角突的是祖先型（越南北部避難所），冰期結束後往西南擴散，新族群演化出長角特徵！
（三）一個冰河期前廣布的祖先被迫在冰期撤退進入兩個避難所，不同避難所族群分化出不同角特徵，冰期過後再一同擴散進入中國西南及緬甸山區。


下一步就是用DIYABC來分析以上三個不一樣的模型哪一個對我的分子資料有最好的解釋度！
用DIYABC開一個新的MSS的資料夾，之後所有結果都會存在裡面：


然後就是導入sequence alignment的資料。這大概是使用這個程式最麻煩的地方，他有自己的format：
第一行是說雌雄比例為1:1
然後有三個基因 h3，its2，及co1。
h3及its2為核基因 A 代表autosomal
co1為粒線體 M 代表mitochondrial
之後用Pop來區分不同族群的個體，個體名稱與sequences之間用 , 分開。
不同基因需放在<>裡面相互隔開。同基因的不同alleles放在[]裡面相互隔開。


如果你的format都正確，倒入mss file之後會進入下面這一個畫面：
我通常先開始設定歷史模型，把historical model下面的set按下去就好。

右邊的add scenario按鈕可以讓您輸入增加更多的歷史模型，但是在這裡我玩三個就好。
這裡的Scenario 1, 2, 3剛好對應我上面想到的三個歷史模型（如果想知道N1, N2 varNe等argument的意思請參照DIYABC的手冊）：

如果你不確定自己arguments寫得對不對，可以按右邊的check and draw scenario鍵來看看：

在確定了這三個不一樣的分化過程就是對應我的三個不同歷史模型之後，下一步就是設定priors：
這一步很重要，如果用不對的族群量大小或分化時間的範圍來做電腦模擬，那結果就不可信，切記！
prior setting是ABC分析裡面最重要的一環，設錯了就全錯了！！！

ABC的分析絕對沒有跑一次就可以了。我的例子用default setting的effective pop prior，但是初步跑完後都需要再改：


Historical model跟priors設定好了之後就要換到設定Genetic data的選項：

我首先把我的分子資料分成兩個群，一群核基因（h3 + its2）跟另一群粒線體基因（co1）：
之後需要分別設這兩個群的分子演化模型（modeltest的結果），即決定要選用哪些population genetics的summary statistics來協助MCMC search（檢測電腦模擬的結果是不是好）。

以粒線體co1這個群組為例，我設定了TN+I+G的演化模型，然後summary statistics我選了tajima's D跟各族群間兩兩比較的FST：


當然你可以選很多很多的summary satistics來假裝你的電腦simulation一定要很接近的結果你才要。但是要記住這樣的話你的電腦模擬可以永遠都得不到任何接近真實data的結果，或是你要simulate非常非常久才能得到一點點可以用的資料。決定權在個人，你有辦法說服reviewers就好。


該設定的東西都設定好之後就按下validation 鍵，沒問題就會回到主畫面：

下一步就是決定MCMC要跑多久，我的例子就照這個program給我的建議好3百萬代（每個historical model各1百萬代MCMC）：
program人很好的跟我說大概要跑14分鐘。



跑完MCMC模擬三百萬代之後終於可以分析資料啦～終於可以按下Analyses這一個tab：

首先根據我們真實得到的資料跟電腦模擬資料計算出來的summary statistics跑個PCA：

我的結果還好，模擬資料的結果（小點點）都還在真實資料（黃色大點點）的附近，不會差太遠：


接下來就是做model comparison！！！我們最想知道的結果：
根據一定的誤差範圍（自己定可以接受的範圍），歷史模型（二）最常產生接近真實資料的結果！



到這裡，很高興了！得到了結果顯示（二）是三個歷史模型中最可能的model！！
但是這樣就結束了嗎？就可以說小角角會生出大角角的後代的嗎？
當然不是...............

來看看族群參數的估計結果吧（紅色是prior綠色是posterior，只用model 2的模擬結果）：
有效族群量：

朔源時間：


有效族群量N1 N2等等的posterior都爛斃了...prior設太小了....
想要發表的話就在重跑一次吧，但是uniform prior的上限要開大一點。
其實可以先用 IM或Migrate-N等等program先估計好這些族群遺傳的參數（就是Ne跟tau），然後再根據這些估計結果設定prior，會比較快，不用像我這個例子一樣亂亂試。




就醬，abc的分析已經不難做了！大家對親緣地理，族群分化，種化等等議題有興趣的話，試試ABC吧！

Von guten Mächten
--by Dietrich Bonhoeffer

Von guten Mächten treu und still umgeben,
（良善的力量圍繞著我們，堅定且安祥）
behütet und getröstet wunderbar,
（受到保護且美好的安撫）
so will ich diese Tage mit euch leben,
（所以我想在這些日子與祢一同生活）
und mit euch gehen in ein neues Jahr.
（並且與祢的同在一起進入新的一年）

Noch will das alte unsre Herzen quälen,
（想起過去發生過的事情依然會讓內心煎熬）
noch drückt uns böser Tage schwere Last.
（未來也一定會持續有不好的事情發生壓迫我們）
Ach Herr, gib unsern aufgeschreckten Seelen das Heil,
（主啊，請醫治我們受傷的靈魂）
für das Du uns geschaffen hast.
（那是你在創造我們之初就給予我們的東西）

Und reichst Du uns den schweren Kelch, den bittern, des Leids,
（當祢遞給我們生命的苦杯，充滿苦及悲傷難過的生活）
gefüllt bis an den höchsten Rand,
（請把這杯子裝到最滿的邊緣）
so nehmen wir ihn dankbar ohne Zittern
（我們將充滿感謝的把這杯子接過，毫不畏懼）
aus Deiner guten und geliebten Hand.
（因它來自你良善及充滿關愛的手）

Doch willst Du uns noch einmal Freude schenken
an dieser Welt und ihrer Sonne Glanz,
（如果祢想再次用這完美的世界和耀眼的陽光來關照我們）
dann woll'n wir des Vergangenen gedenken,
（我們當紀念之前困苦的時光）
und dann gehört Dir unser Leben ganz.
（然後我們的生活將完全屬於祢）

Laß warm und hell die Kerzen heute flammen
（讓今天的燭火燒得既溫暖又明亮）
die Du in unsre Dunkelheit gebracht,
（在我們的黑暗之中，祢為我們帶來這燭火）
führ, wenn es sein kann, wieder uns zusammen!
（如果祢願意，請再次領導我們）
Wir wissen es, Dein Licht scheint in der Nacht.
（我們知道，祢的光能照亮黑夜）

Wenn sich die Stille nun tief um uns breitet,
（當祢將寧靜深深地散播到我們身上）
so laß uns hören jenen vollen Klang der Welt,
（讓我們可以聽到這世界上所有的聲響）
die unsichtbar sich um uns weitet,
（這些肉眼見不到的東西在我們的周遭擴散）
all Deiner Kinder hohen Lobgesang.
（祢的所有子民高聲讚頌）

Von guten Mächten wunderbar geborgen
（從良善的力量之中，美好的事物會被發現）
erwarten wir getrost, was kommen mag.
（我們安心得等待將來）
Gott ist bei uns am Abend und am Morgen
（我主日夜與我同在）
und ganz gewiß an jedem neuen Tag.
（毫無疑問的每個新的一天都將與我同在）

不同類群生物之間的遺傳交流發生的頻率及發生的前提是什麼？
這一直都是演化學研究上未有定論的議題。比如說（一）雜交可能會導致物種的特性消失，最後導致兩個不一樣的生物種混成一種：因此不同生物種間雜交應該是很少發生的事件。但是（二）不同種之間雜交可以讓具有適應性的等位基因（allele）在種間相互交換，可能促進適應性等位基因接受者的類群（allele-receivor lineage）有辦法發育新的表型（phenotype）或進入並使用原本無法使用的生態棲位（niche）：因此種間雜交可能會有適應性演化上的好處。論述（二）很好的一個例子就是龍鱗展鬥：人類用雜交的方法把龍鱗這個特徵從另外一個Betta的種轉到splendens這個種裡面，導致他有極高的適存性（fitness），因為人類會繁殖很多有這個特徵的魚體的子代。

但是討論這些問題之前，其實我們並不了解生物種間雜交發生的前提是什麼，或是有哪些因子可能會限制種間雜交發生。今年我用長戟大兜蟲這個類群來專門探討這一類型問題：（一）地理分佈與種間遺傳滲入的關聯性，（二）親緣關係遠近與遺傳滲入發生的關聯性，及（三）型態適應與遺傳滲入的關聯性。

研究結果顯示不同物種之間是否發生過遺傳滲入（至少在長戟大兜蟲這個類群）是不受親緣關係遠近或是適應性型態相似程度限制的：親緣關係差很遠的類群（例如長戟septentrionalis跟白兜maya彼此之間的最近共祖要追朔到4百萬年前）一樣可以有遺傳滲入，長相跟體色（與環境適應有關的型態）差很多的不同物種之間也可以發現有遺傳滲入。
種間遺傳滲入的唯一前提就是兩個種必須要有地理分佈上的重疊（很像是廢話，但沒有人真的去檢測過這個），但是地理分佈上有重疊的物種並不一定會有遺傳滲入發生。

說了這麼多的結果，其實就只是顯示遺傳滲入發生的原因跟限制是什麼還是不清楚的...但是在長戟的例子裡有一個很有趣的發現，那就是在中美洲有兩種共域生活的物種（septentrionalis跟maya），他們之間有很顯著的遺傳滲入發生，且遺傳滲入是單方向的由septentrionalis進入maya。maya這個物種的體色與其他典型白兜很不一樣，但是卻與長戟沒有統計上顯著的差別。這個發現指向很可能maya及他的姐妹種moroni的特殊體色（在熱帶地區有適應性的功能：camouflage）可能是藉由雜交從septentrionalis偷過來的。當然在沒有全基因組的資料及功能性的研究之前，這都還只是推測，期待之後的基因體資料可以有更好的分析來回答這個適應性遺傳滲入（adaptive introgression）的問題。

Mysterium Fidei.（信德的奧跡）
Mortem tuam annuntiamus, Domine,
（基督，我們傳報祢的聖死[祢的死被宣告，我主]）
et tuam resurrectionem confitemur,
（我們歌頌祢的復活[及祢的復活被承認]）
donec venias.
（我們期待祢光榮的來臨[直到祢再來]）

如何客觀量化地辨識生物種是生物學上最基本的問題。因為生物學上絕大部分的問題都必須仰賴清楚的分類結果（例如保育政策制定，生態學及演化學的研究都會因為分類學的結果不一樣造成不同的結果）。但是因為生物種的定義百百種，所以既使近年來量化辨識生物種的統計方法發展很多（絕大多數都以分子分類為基礎），拿到同樣的採樣甚至是完全一樣的資料還是可能的到不一樣的的生物種辨識結果。

分類學上並不是只有在定義及辨識”生物種“上分歧很大，還有另一個更大的問題存在叫做“亞種”。傳統上不同亞種是用來命名一個多態型（polytypic）生物種的不同形態群（大多數例子存在于地理隔離的不同形態群）。但是大部份的亞種分類辨識的過程跟結果是非常不客觀也不是根據可以量化分析的結果（有興趣可以看Mayr及Ashlock的書）。那我們現在很多資料庫裡面的不同亞種跟不同種真的代表不一樣的東西嗎？
這問題很重要，因為有些資料庫會把亞種當成同物異名：如果完全不一樣的東西被當成是亞種，那我們可能實際上少計算也少關心了很多特別的生物。

我今年剛有一篇文章發表就是用一群甲蟲（Dynastes屬）來檢測被定義成不同種跟不同亞種的分類群是不是真的有生物學上相對應的不同關係。
結果顯示並沒有。

這篇研究主要的發現及討論有：
1. 就統計學上來說，不同種之間的差異並沒有量化的比不同亞種間的差異大。
2. 不同的分類群之間的關係常常都代表著不同量的種化程度（different positions along the speciation continuum）：例如兩個現被定義為不同種的Dynastes granti跟D. hyllus只有分子生物學的資料可以顯示他們是不同東西，但是形態學跟生態學上這兩者並沒有顯著的差異。另一方面現被當作不同亞種的D. hercules lichyi跟D. h. ecuatorianus不止分子生物學資料顯示兩者完全不同，他們在形態上也有顯著的差異並住在不一樣的森林類型之中。所以不同分類群（taxa）之間可以有很多也可能只有一種資料顯示他們顯著分化。
不過很重要的一點就是：在Dynastes屬甲蟲，地理分佈的重疊只出現在分子及形態學（或生態學）資料上都顯示已經完全顯著分化的不同分類群之間。沒有全部都分化的不同分類群在地理上都是異域（allopatric）分佈的。
3. 生物種界定必須要納入所有可以得到的資料一起討論。如果我這篇文章只跟大部分其它近期的文章一樣只針對使用分子資料（2004年至今），那我就完全沒辦法得到點2.裡面的資訊了（i.e., speciation is a continuous process）。
****我真的超肚爛純粹使用分子資料做分類學的人****
4. 是時候該開始客觀量化的探討比較不同類群到底是把什麼樣的東西當成不同種了。例如鳥類學家認為的不同種跟魚類學家認定的不同種的原則是一樣的嗎？如果不一樣，差別在哪？
5. 系統分類學者不能只是給給結果然後打嘴炮，必須要把結果寫成大家都能懂的東西才能對自己有興趣的分類群真的有貢獻。很多文章title是"molecular species delimitation"的最後都只有提出對現在分類結果的建議但是沒有真的做出重新整理分類地位，這樣就結束研究對其他生物學的研究（如保育生物學及生態學）是很難有影響的。如果得到有用的資料，就應該要好好的將生物誌或分類回顧寫好刊出來當一個負責任的生物學家。因此，根據我這一篇文章的結果，Dynastes屬的生物誌也將被發表（正在被review中），其重點會是正式的把所有的Dynastes屬內的亞種全提升為不同種。

如果有興趣，我的文章可以在下面的網址找到（也可以寫email跟我要）：
http://sysbio.oxfordjournals.org/content/early/2015/12/16/sysbio.syv119.short?rss=1
所有的分子形態及分佈地點（生態學及地理分佈重疊）的資料都可以在Dryad資料庫免費下載：
http://www.datadryad.org/resource/doi:10.5061/dryad.8p6m0

生物有很多不一樣的性狀，有些性狀可以很簡當的量化，例如體重或是體長等等可以輕易的測得單一的值，然後可以進行所有母數類型的統計分析。如果這些類型的性狀是研究員有興趣的，那這個研究員真的是很幸運，可以很直接的測得想要的資料（當然也有個前提就是要有足夠的樣本數）。
在生物學（尤其是演化及系統分類學），很多有趣的性狀其實都是沒辦法直接測得可以分析的“值”。例如不同隻鳥有不一樣形狀的鳥喙，有的是尖尖長長的，有得卻很像勾勾。其他又像葉子的形狀可以很不規則，那我們到底要怎樣比較這些差異呢？

Kuhl與Giardina在1982年發表的一篇文章，大大的改變生物學上對形態分析的方法：本來只能用文字表達的個體差異現在可以轉化成數值做統計分析。例如本來你只能說種A百分之90是正三角形百分之10是等腰三角形，種B是矩形跟直角三角形之類的，現在都可以轉換成主成分分析。
R有一個很棒的package叫做Momocs可以很方便地做elliptical Fourier（橢圓傅立葉轉換？）類型的分析（此外Momocs也可以做tangent angle及radii variation Fourier行的分析）。這個package雖然不像geomorph可以直接處理圖檔，但是處理圖檔並不困難，下面就用最容易的powerpoint來介紹怎麼把jpg檔的圖檔裡面我們要的shape變成XY-coordinates給Momocs分析。

首先先打開PPT並畫一個黑色的正圓（這應該很簡單）：

用這個目的是讓之後的每個圖都可以在同一個XY-dimension中。黑色的圓可以當所有形狀的同源中心：大家都有這個一樣的地方，之後做elliptical Fourier分析需要對齊，有一樣的XY象限，又有一樣的point of origin，才能確定我們分析的真的是形狀上的不一樣，而不是因為不同角度造成的不一樣。

畫完元之後就把你的照片貼到PPT上（把你有興趣的地方疊到黑色的圓上面）：


然後把這個檔用PPT直接存成jpg檔之後用imageJ打開，然後用process裡面的binary把你的jpg檔變成binary format（黑白檔）：

然後在toolbar裡面按下這一個仙女棒：

你就可以發現你的圖的otline被選起來了(黃色的線）！
選起來之後就把這個檔用imageJ存成XY-coordinates就完畢了，Momocs就可以吃這個檔案了：


在做elliptical Fourier之前一定要做的就是讓所有檔案的point of origin變成一樣的及center the file。Momocs把txt檔（你剛剛存的XY-coordinate）裡面的第一個XY值當成是the point of origin，所以我們需要手動改一下所有的txt檔。還記得所有的檔都有同樣的圓嗎？我們可以在圓上面隨便選一個點當the pont of origin，因為所有的圖都會有那一個相同的點：

你可以用各種不一樣的text file editor來改txt檔，我用wrangler來改我的檔案，讓所有txt檔的第一個XY-coordinate都是696 260：



改完txt檔之後就可以用Momocs把所有的XY-coordinates都讀到R裡面去了（我隨便選了四個不同種的長戟，各有2到4個大型雄蟲個體的頭角形狀）：
#################
記得要center記得要center記得要center！！
#################


檔案讀進去之後下一步就是要決定用多少harmonics（等距離分佈的點）來描述你的形狀！這很重要，用越少的harmonics分析越快，但是形狀跟原本的形狀就越不像，下面用四個例子來看4，8，16，及32個harmonics對不一樣的outlines的描述：

用眼睛看起來32個harmonics比其他數量的harmonics對輸入的outlines的描述好很多。
但是除了用眼睛看之外，其實還是可以用統計方法來看，下面我就分享兩個Momocs裡面常用的方法（上圖看如果你在新的outline上面隨機選幾個點，他會與原本的outline相差多遠，下圖就是算harmonic power）：

32個harmonics真的比其他好很多，所以就先用32來玩！（harmonics沒辦法比你輸入的XY-coordinates多就是了，最多就是用你的original XY-coordinates來當harmonics）。

先用eFourier這個function可以把存進去的outlines一起做elliptical Fourier的分析：



得到harmonics的coefficients之後就可以根據這些值來做主成分分析了：


是不是很方便呢？除此之外還可以看每個PC的描述範圍，看在那個PC上面的不一樣outlines可以差多少：


如果真的想做形質（多變量）分析也可以用Manova之類的分析，或是就用hierarchical clustering來做80年代分類學家會做的nummerical taxonomy：



現在很多博物館都把館藏標本數位化！這代表著既使我們沒辦法採集到一些樣本，我們還是可以做很多形態演化的分析的！這個連結裡面有一個很有趣的研究，他們就是用很多小提琴的照片（很多是不一樣的館藏經典不同年代的小提琴設計）來分析不同時代設計家演奏員到底是用什麼不一樣的小提琴（小提琴形態的演化）。

博物館館藏數位化是很重要的工作，很多博物館都已經分享的大量的資料。但是除了大家分享資料之外，我們也要想辦法讓這些資料變成研究交流的重點，而不是放在那裡沒人用。Momocs就是一個可以活用大家已經分享的數位資料的package。

大尺度演化學的研究中一個很有趣的問題類型就是：
（一）為什麼某些時期系群的分化速度會變快？又有些時候會變慢？
（二）為什麼某些系群分化的速度會比其他系群快？

為了回答這些問題，首先需要確認的就是：（一）不同時期的分化速度真的有顯著的改變嗎？（二）分化速度在不同系群真的有顯著的不一樣嗎？
雖然很多時候重建好一棵親緣關係樹之後很常發現其中不同的系群（lineages）常常有不一樣的物種數。或是生物的分歧（lineage splits）似乎都集中在某些特定的時候。但是這些都只是基於看圖說故事，這些故事真的有其支持度嗎？

就拿我研究的分類群Dynastes屬甲蟲的物種關係樹來做例子：
很明顯得Dynastes亞屬比Theogenes亞屬多了很多種生物（假設所有的物種這棵親緣樹都有），並且大部份的分化都是在紅色那一條線的時間點之後才發生的！那在這個分類群中是不是有shifts in diversification rates呢？


一個最早被拿來用的分析方法是用R的package laser做lineage through time plot（LTT plot）。
以下我先很單純的根據pure birth/Yule及birth death models來看一樣16個物種如果沒有rate shift的話，在一樣的分化時間內會有怎樣的lineage through time plots。

黑色的線條是真正第一個圖Dynastes species tree的LTT plot，不一樣的顏色部分代表根據一個model（左：pure birth，右：birth death）預期的LTT plot會長什麼樣。
這裡的結果大概就是說如果是pure birth model，現在的Dynastes tree的diversification pattern沒有辦法被單一一個diversification rate fit很好。Birth death model似乎有fit比較好一點，但是還是在時間點-2以內沒辦法fit很好。

然後呢？既然單一一個rate不能fit真的data很好，那我們就試試看fit這個Dynastes LTT plot多個rates，看看會不會比較好（以下例子都是用Yule model）！

很簡單的就比較不同models的AIC values，我們發現fit三個不一樣的rates（兩個rate shifts）有最佳（最小）的AIC，那大概就是最好的model！


LTT是個很簡單的model，基本上就是看時間對bifurcation event發生的線性回歸。到底是單一一條直線回歸就可以描述data，或是需要在不同的時間內用不一樣的斜率？但是回來看例子的Dynastes tree，Theogenes只有兩個物種，哪裡有足夠的data可以讓model算rate呢？那Theogene與Dynastes亞屬之間的分化速度到底要怎麼比呢？

Modeling Evolutionary Diversification Using Stepwise AIC (MEDUSA)就是一個這樣的model來偵測一個phylogeny上面有沒有很顯著的diversification rate shift。
MEDUSA的原則就是先給整個Dynastes tree fit單一一個diversification rate然後計算這個單一rate model的AIC value。之後隨機的挑選某個node（或是stem），讓那個node之前及之後fit不一樣的diversification rates，然後再計算AIC value。如果fit多個rate的model的AIC明顯的比fit比較少rate的model的AIC好（用AICc threshold），那就留顯著較好的model，然後再繼續嘗試在不同的nodes前後fit不一樣的rates。這樣的步驟一直重複直到較複雜的model的AIC value不再明顯的比簡單的model來的好就停止。

MEDUSA的分析現在可以使用R的package geiger就好，但是我習慣用turboMEDUSA。Joseph的GitHub site有詳細的介紹如何操作，所以以下我就直接分享MEDUSA分析Dynastes tree的結果（假設有完整的sample沒有missing data）：

Theogenes lineage可能發生過rate shift，而且是rate decrease！

turboMEDUSA的好處就是不只可以比較一棵親緣樹。如果你有很多彼此之間沒有統計上顯著差異的樹（例如Bayesian analysis的post-burnin之後的樹），全部都可以拿來一起分析。如此就可以了解如果加入phylogenetic reconstruction的uncertainty，Dynastes甲蟲的分化過程到底有沒有發生過diversification rate shift。
另外一個用MEDUSA分析的好處就是可以加入missing data的計算。舉例來說，在我的Dynastes屬甲蟲的研究中，我並沒有採集到D. hercules tuxtlas，D. hercules takakuwai，及D. neptunus rouchei的樣本。如果這三個族群是不一樣的independently evolving lineages，那現在這個Dynastes的species tree就其實有三個missing data。MEDUSA的分析可以加入一個richness的data，用來表示某些類群可能有多於一個種。比如說我雖然不知道D. hercules tuxtlas是不是一個獨立的種，但是我知道他的關係與D. hercules septentrionalis很相近是姐妹種，那我就可以在richness的file裡面把Dhs的richness設為2來表示那個地方應該其實有兩個種。MEDUSA分析的時候就會在Dhs那個lineage上隨機加上一個split來表示那裡有多一個分化事件。Richness file的例子如下：


下面這個圖就是分析8000棵post-burnin species tree的結果：

node上的數字代表多少百分比的post-burnin trees支持那個node發生rate shift。左邊的圖是假設complete taxon sampling，右邊的圖是加入假設有三個missing data的訊息。
大部份的tree是支持在Theogenes發生分化速度減緩，不過也有很多是支持在Dynastes發生分化速度加快。但不管如何這裡得到很明顯的訊息就是Dynastes亞屬比Theogenes的分化速率快（至於是減速還是加速就各有支持度）。



看到了這裡，其實很多研究需要的分析統計支持都可以得到了，那還有什麼可以改進的地方嗎？當然有的！親緣關係分析方法中最常見的likelihood vs. Bayesian在這個時候就出現了。
Bayesian Analysis of Macroevolutionary Mixtures（BAMM）就是一個使用reverse jump的方法讓MCMC search可以在不同rate shifts的models之間來回估計不同models的可能性。Dan的網頁上有簡短的解釋這樣的分析與根據delta AIC value差異選擇model的差異性：簡短的來說，很多不一樣的model算出來的likelihood（或換算的AIC）可能非常的接近，只選擇其中一個來表示整個dataset不一定有考慮到所有的可能性。BAMM可以把所有equally likely的models整合起來一起看，比較客觀不失公平。

Dan的網頁有很詳細的介紹怎麼跑分析，需要的東西就只有去下載的主程式bamm，要被分析的親緣關係樹（nwk format），還有一個control file。control file可以在網頁上下載，下載之後只需要自己決定MCMC要跑多久及MCMC的sample frequency。我初步對於Dynastes species tree的BAMM分析的設定如下：

跑完之後用R的BAMMTool package檢視結果，MCMC幾乎馬上就plateau了：

把前20%的samples去掉burnin，基本上估計的參數都有很好的effective sample size：

大部份的MCMC samples要不是支持一個shift的model（45%）就是支持沒有shift的model（30%）。

然後看一下所有post burnin MCMC samples平均估計出來每個lineage的分化速率：

是不是跟MEDUSA的結果很像呢？Dynastes亞屬比Theogenes快。但這裡還指出另外一個有趣的點：長戟大兜這個clade似乎分化速度又比白兜clade要來得快。

最後看一下Macroevolutionary cohort analysis：

這分析主要是看兩倆物種之間的分化速率比較。如果兩個物種之間有較高機率出現一樣的分化速度，他們就會顯示比較紅色的顏色。在Dynastes這一棵species tree上很明顯的同一個lineage裡面的物種估計出來的分化速度也比較相近。


如果lineage specific diversification rate或是diversification rate shift會是你很有興趣的pattern，MEDUSA跟BAMM基本上都是非常合適好用的分析方法。需要煩惱的就只剩下likelihood vs. Bayesian的事，那就是統計學與哲學上的問題了。
我可以講的就是MEDUSA就是純粹看diversification rate，很清楚，而且非常容易可以accommodate missing data（如果你知道有missing data）。壞處就是一般likelihood analyses都會有的壞處：慢，會stock在某個local maximum。BAMM非常的快（因為是C），可以連Trait data一起分析（例如其實有的時候你覺得某個key innovation導致分化速度變快，但是實際上可能只是因為其中某個lineage intrinsic的分化速度變快也可以被當假說檢測）。缺點就是：要跑多少MCMC及priors要怎麼設定（所有Bayesian分析都有的問題）？

一個特徵（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造成在不同類群之間有不同程度的分化！

除了可以研究不同生物之間的親緣關係之外，不同生物之間的形態差異很多時候還代表了不同生物在分化之後在不一樣的環境下面適應的結果。所以如何量化分析形態上的差異，一直都是演化學和系統分類學研究的主要研究方向之一。其中一個方法可以將不同形狀量化成值的做法就是利用Geometric morphometric analysis。

geomorph是一個R用來做geometric morphometrics分析的package，只要研究員有已經照好的照片（以2d的data為例，但是這個package也可以做3d的影像處理跟分析），就可以一條龍的從選擇landmarks到Procrustes analysis一直到PCA全部完成（link）。非常的方便！
因為我用的背景是Mavericks，想在這個背景之下執行geomorph需要的tcltk已經沒有OSX裡面，所以如果你是跟我一樣喜歡使用Xcode 5或是不小心已經update變成Mavericks，你需要先去下載安裝新的X11（XQuartz）才可以跑geomorph這個package。
安裝geomorph的方法就跟其它安裝R package的方法一樣，只要在R 輸入：
install.packages('geomorph');

接下來在R用Misc裡面的changing working directory找到放圖檔的資料夾就可以開始把選好的Landmarks存成.tps檔了！
先把geomorph這個package的所有functions叫到R裡面：
library(geomorph);
之後在R裡面輸入下面這一個指令：
digitize2d('DhsCR2.jpg', 9, 10, tpsfile = 'Dynastes.tps');

digitize2d是分析圖檔用的function。
filename.jpg是告訴R你要用哪一個圖檔分析。我的例子是DhsCR3.jpg這個圖檔。
第三個數字9是告知這個function這個圖中我要存9個landmarks。
第四個數字10是說我會告知這個程式10 minimeters在圖片裡面是多長。
最後一個是跟程式說執行完了之後把結果存成.tps檔，名稱叫做Dynastes.tps。
（我的例子總共分析15個.jpg檔，全部一起存在Dynastes.tps中）

一切順利的話圖檔就會在R開啓（如下圖）：

接著我們要在那條黑色的scale bar前後各按一下（用滑鼠左鍵按，沒錯就是這麼簡單），告知R這兩個紅色圈圈之間代表長度10mm（R prompt出現scale bar 正不正確？你就按y之後按return跟它說對）：

接下來就是一個landmark一個landmark的按。每按一個之後R prompt都會問你一次你點的landmark是不是你要的（這個例子都是）：

另外一個例子，用不一樣的種類（DheE2）：

全部都按完了之後，資料夾裡面就會出現儲存好的.tps檔了，可以用Wrangler打開來看裡面長什麼樣子：


等到將所有的.jpg圖檔都用digitize2d將landmarks的資訊存成.tps檔之後，就可以先對這些landmarks做procrustes analysis。首先用readland.tps這一個function將Dynastes.tps這個檔案叫進去R：
lms <- readland.tps('Dynastes.tps', specID = 'ID');
lms是之後要用的所有landmarks檔案我給它的名字，specID是跟R說每個個體的種名用tps檔案裡面的ID表示。

再來就是generalized procrustes analysis：
lms.gpa <- gpagen(lms);
下圖就會出現在R console：灰色的點是原始資料，黑色的點則是不同個體之間procrustes analysis之後的那個landmark的anchor：

設一下依照順序出現的species的顏色：
spid <- c('gold','gold','gold','gold','gold','gold','brown','brown','brown','brown','green','green','green','green','green');
這個例子用了6個ecu（金色gold），4個occ（棕色brown），及5個sep（綠色）物種的image data。
然後使用plotTangentSpace這個function就可以得到大家喜歡的geometric morphometric analysis的圖了！
lms.pca <- plotTangentSpace(lms.gpa$coords, label = NULL, group=spid);


在R console就打lms.pca，下面的表就會出現就可以知道每個PC解釋Data的百分比有多少：

相當的方便！PC1就解釋五十多百分比的變量，PC2二十多百分比！

最後就是對data做anova的分析看看有沒有顯著差異在各個不同種間：
lms.2d <- two.d.array(lms.gpa$coords);
procD.lm(lms.2d~spid);
下面的結果就會出現在R console（根據999個permutatios，如果想用不一樣的permutation number就在function裡面加上iter = ###的只是就好，不同種類有顯著的不同）：

然後再做pairwise的比較：
lms.pd <- pairwiseD.test(lms.2d~spid);
就可以得到兩兩相比的結果（1: brown, 2: gold, 3: green：注意這邊結果會是alphabetical的用123排列）：

ecu（gold）跟其它兩個種類很不一樣（P = 0.058, 0.015），occ跟sep就分不開來（P = 0.642）。

很方便的一個package！
當然還有更多新的分析可以參考作者自己的網頁：
http://www.geomorph.net/

生物適應不同環境的過程時常就像學生參加考試一樣：當遇到一個特定的題目出現的時候（Environmental Challenge），不一樣的人會有不一樣的方法來解題。（一）有些人考試之前已經做過很多考古題，什麼不一樣的題目出來都看過了就用原本準備好的方法解題（Standing Genetic Variation）。（二）有些人就是慧根高，遇到沒看過的題目還是可以想到新的解題法回答問題（New Mutation）。以上這兩種方法最常見，也最常被拿來互相比較優劣及對環境改變時適應速度的影響。但是其實還有另一種不常被拿來一起討論的方法，那就是（三）偷看別人的答案（Introgression）。因為人類社會通常不喜歡作弊的行為，而動物之間長久以來（近來有改變的想法）很多學者都假設分化久遠的物種很少會相互雜交，偷看答案（introgression）vs自己做答（standing variation and/or new mutation）很少被拿來一起比較。

最近因為二代定序技術的進展以及價錢的降低，還有基因體分析技術的發展，越來越的例子顯示其實不同動物物種之間的遺傳交流是會發生的，而且有的時候會有適應性的基因在不同物種之間傳布。大概就像現在大家都知道三星手機跟頻果手機還有其他牌子都手機設計，硬體，及軟體都互相抄來抄去一樣，這個解決問題的方法可能是很常見的。大家也就越來越常研究討論這類型的問題。

我的研究物種之中也存在一個系群它的兩個物種身體的顏色與它的姐妹種（ty: tityus）非常得不一樣，但是與這兩個種（下圖的ma & mo）共域的其它種類卻與它有相似的體色（長戟，下圖用sep為例子:septentrionalis）。因為這類甲蟲的體色被當成是保護色可以融入森林背景顏色之中，當ma （Dynastes maya）被發表為新種的時候，其發表者就討論到ma可能是熱帶雨林性長戟類及亞熱帶北美白兜的中間型，也許有可能是有雜交也說不定（Hardy 2003）。
對我的研究來說，這就是一個很有趣的問題了：
共域生存（坐隔壁的同學）的兩個分化已久的物種會相互交換基因嗎（抄解答嗎）？


回答這個問題的方法其中之一就是用Patterson's D。它主要用來偵測兩個不同物種之間是否有過高的shared serived alleles（兩個學生之間是不是有過高的一樣問題解法跟答案）。如下圖所示，如果一切分化過程（解題過程）兩個物種（2 & 3）都是獨立進行的，把這兩個物種與第三個物種（1）及外群（解答）比較，2與3使用相同的新穎解法（與解答本的解法不同）的量（ABBA）應該會和1與3使用一樣新穎解法的數量（BABA）一樣多。但是如果2與3之間有發生introgression（互相交換偷看答案），則ABBA的量就會明顯的多過BABA的量。


這一個方法後來再被擴充成partitioned D。如果在species 3(s31)的姐妹群也有得到需要的分析資料（s32），那我們還可以推估introgression（我想中文是翻譯成遺傳滲入）的方向！
假設遺傳交流方向是由s31進入s2，我們應該會得到顯著結果的D1（excessive shared derived alleles between s31 & s2），但是不會有顯著的D2（that between s32 & s2）。除此之外，s31及s32這個lineage特有的derived alleles也應該會在s2被發現很多（顯著的D12），因為s31遺傳滲入s2的時候會把來自s31與s32共通祖先的unique alleles一起帶過去。



因為遺傳滲入的alleles可能在接收的物種內造成很多異形合子（heterozygous），我需要有很高可信度的定序資料才能做這個分析（如果只sequence兩個reads 然後是不一樣的sequences就把它們稱作heterozygous，在現今Next-Gen定序的錯誤率之下，我是不敢相信），我就把我的library資料重新用pyRAD整理過留下quality很高的結果而已（totally 32320 loci）。拿來分析的文件一小部分如下（pyRAD的.loci output）：

每個block都是一個不同的locus，每個都可能有一些missing data：上圖例子中的locus有外群，有我有興趣的那兩個中美洲的中間型（hypothesized allele receiving species）也有與它共域的物種（allele donor）的資料。但是呢！！！這個locus剛好沒有姐妹種們的資料（殘念），所以不會被拿來分析。就將，一定都要大家都有被sequenced到的loci才會拿來分析。


分析完的結果如下圖（共三個）：

先看Patterson's D，Z score是4點多，統計上非常的顯著，所以ma+mo一定跟sep有著正常關係以外的特別關係（八千多個loci的資料）！
然後看partitioned D，由sep進入ma+mo(2)及相反方向由ma+mo進入sep(3)都有顯著的D12及D1但是D2都不顯著（七千多個loci的資料）！所以中美洲長戟sep跟中美洲白兜ma+mo之間很可能是有雙向的遺傳滲入的，而不是只有單方向給introgression。


最後來試著回答Dynastes屬研究員們長久以來的問題，這些住在臨近的地方而且又長得很像的物種之間很可能真的是有introgression發生的！下一步需要釐清的就是（一）這樣的遺傳滲入真的是造成convergent adaptation的原因嗎？還是（二）introgression歸introgression，身體顏色適應上的區同演化不一定是由introgression帶進來的？
這兩個問題就不是現在我有的資料可以回答的了.....