介紹

還是接上之前Orthofinder分析后，很多老師和同學(xué)都會(huì)選擇做時(shí)間分歧分析。

Orthofinder下篇

Orthofinder上篇

"時(shí)間分歧" 簡(jiǎn)單來(lái)說(shuō)，就是確定物種分化的大概時(shí)間區(qū)域。

分歧時(shí)間區(qū)域計(jì)算原理：基于已知的時(shí)間節(jié)點(diǎn)（已被報(bào)道的分支時(shí)間節(jié)點(diǎn)/具有化石時(shí)間）和分子鐘計(jì)算，推算系統(tǒng)發(fā)育樹的其他節(jié)點(diǎn)發(fā)生的時(shí)間。

簡(jiǎn)單補(bǔ)充分子鐘原理：任意兩個(gè)物種，自從分化成兩個(gè)物種后，他們之間的遺傳差異的進(jìn)化速率應(yīng)該與分化的時(shí)間保持相對(duì)穩(wěn)定。舉個(gè)例子：人與青蛙的分化點(diǎn)肯定比人與猴子的早，那么人與猴子的序列差異理論上會(huì)比人與青蛙的序列差異要小。

目前，采用估算物種分歧時(shí)間的程序常見有：mcmctree、r8s..

這次推送主要分享mcmctree的使用方法~

運(yùn)行

PAML安裝

mcmctree主要是在PAML軟件包。

conda安裝~

兩個(gè)方法下載的paml的軟件包都是4.9版本的~

主要用到paml軟件包的mcmc，其中學(xué)習(xí)手冊(cè)的介紹：Bayesian estimation of species divergence times incorporating uncertainties in fossil calibrations (mcmctree).

paml的學(xué)習(xí)手冊(cè)：http://abacus.gene.ucl.ac.uk/software/pamlDOC.pdf

準(zhǔn)備文件

我們?cè)谏洗瓮莆囊呀?jīng)完成了orthofinder的分析，所以從其分析結(jié)果獲mcmctree運(yùn)行所需的文件。

主要需要兩個(gè)文件：

1.單拷貝直系同源基因數(shù)據(jù)集合的基因ID?

2. 物種的系統(tǒng)發(fā)育樹文件

1. 準(zhǔn)備序列比對(duì)文件

2. 準(zhǔn)備樹文件

在進(jìn)行分歧時(shí)間預(yù)測(cè)，我們至少知道一個(gè)節(jié)點(diǎn)的大概分歧時(shí)間

TIMETREE數(shù)據(jù)庫(kù)：http://timetree.org/?是一個(gè)記錄多個(gè)植物分歧進(jìn)化時(shí)間的公共數(shù)據(jù)庫(kù)。而且提供相關(guān)文獻(xiàn)，我們可以進(jìn)一步進(jìn)行確定

這邊我們暫定葡萄與擬南芥/柑橘的分歧點(diǎn)約為大于 100

參考o(jì)rthofinder的系統(tǒng)進(jìn)化關(guān)系

稍微解釋一下如何編寫這個(gè)物種樹，借鑒Orthofinder的構(gòu)建好的樹結(jié)構(gòu)，補(bǔ)充了葡萄分歧時(shí)間大于100MYB。

也可以直接使用?'L(1.0)'、'>1.0'

如果是在某個(gè)時(shí)間段可以：'B(1.0,1.17)'、'>1.0<1.17'

如果是小于某個(gè)時(shí)間節(jié)點(diǎn)可以：'U(1.17)'、'<1.17'

而括號(hào)后的兩個(gè)參數(shù)分為 p=0.1, c=0.2, 由于物種數(shù)量比較小而且距離還相對(duì)較遠(yuǎn)，增加兩個(gè)內(nèi)置函數(shù)更嚴(yán)格限制已知的分歧點(diǎn)。

更詳細(xì)請(qǐng)參考PAML User Guide的 第9節(jié)mcmctree的Fossil calibration 部分。

mcmctree運(yùn)行

簡(jiǎn)單說(shuō)一下，最主要填寫的項(xiàng)：

seqfile = all.singlecopy.phy

treefile = treefile

outfile = out *輸出前綴

ndata = 30

建議將隨機(jī)采取的比對(duì)序列情況設(shè)置大一點(diǎn)；因?yàn)槲疫@里六個(gè)物種單拷貝基因數(shù)量太多了，運(yùn)行起來(lái)耗時(shí)多，所以我僅采用30個(gè)情況進(jìn)行時(shí)間分歧計(jì)算。

usedata

設(shè)置是否利用多序列比對(duì)的數(shù)據(jù)：0，表示不使用多序列比對(duì)數(shù)據(jù)，則不會(huì)進(jìn)行l(wèi)ikelihood計(jì)算；1，表示使用多序列比對(duì)數(shù)據(jù)使用Felsenstein（1981）的似然計(jì)算，也是MCMC常使用的參數(shù)，不過計(jì)算速度較慢; 2，進(jìn)行approximation likelihood分析，此時(shí)不需要讀取多序列比對(duì)數(shù)據(jù)，直接讀取當(dāng)前目錄中的in.BV文件。該文件是使用usedata = 3參數(shù)生成的out.BV文件重命名而來(lái)的。

clock = 3

model = 4

0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85 說(shuō)明書列舉了四種情況，其中JC69是適合比較簡(jiǎn)單的關(guān)系，HKY85是相對(duì)比較遠(yuǎn)源的關(guān)系（更復(fù)雜的模型，配合alpha=0.5使用）。

alpha = 0.5

首先設(shè)置usedate=3

獲得 out.BV文件

usedate=2

查看FigTree.tre

在Figtree可以可視化，后期導(dǎo)出SVP。

至于如何更美化這個(gè)分歧樹... 就各自修行。

mcmctree.ctl配置

seed = -1 ? ? ? seqfile = all.singlecopy.phy ? ? ?treefile = treefile ? ? ? outfile = out ? ? ? ? ndata = 20 ? ? ? seqtype = 0 ?* 0: nucleotides; 1:codons; 2:AAs ? ? ? usedata = 2 ? ?* 0: no data; 1:seq like; 2:use in.BV; 3: out.BV ? ? ? ? clock = 3 ? ?* 1: global clock; 2: independent rates; 3: correlated rates ? ? ? RootAge = '<1.0' ?* safe constraint on root age, used if no fossil for root. ? ? ? ? model = 4 ? ?* 0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85 ? ? ? ? alpha = 0.5 ? ?* alpha for gamma rates at sites ? ? ? ? ncatG = 5 ? ?* No. categories in discrete gamma ? ? cleandata = 0 ? ?* remove sites with ambiguity data (1:yes, 0:no)? ? ? ? BDparas = 1 1 0 ? ?* birth, death, sampling ? kappa_gamma = 6 2 ? ? ?* gamma prior for kappa ? alpha_gamma = 1 1 ? ? ?* gamma prior for alpha ? rgene_gamma = 2 2 ? * gamma prior for overall rates for genes ?sigma2_gamma = 1 10 ? ?* gamma prior for sigma^2 ? ? (for clock=2 or 3) ? ? ?finetune = 1: 0.1 ?0.1 ?0.1 ?0.01 .5 ?* auto (0 or 1) : times, musigma2, rates, mixing, paras, FossilErr ? ? ? ? print = 1 ? ? ? ?burnin = 2000 ? ? ?sampfreq = 10 ? ? ? nsample = 20000