分歧時間是目前進化分析的一個熱點，以某幾個特定類群的化石時間作為校正，然後通過基因序列間的分歧程度以及分子鐘來估算物種間的分歧時間，同時估算系統發育樹上其它節點的發生時間，進而推斷相關類群的起源和不同類群的分歧時間。目前，采用依據BI樹和ML樹估計物種分歧時間的程式很多，例如R8S、MCMCTREE、MULTIDIVTIME、BEAST、MEGA等，通過不同的政策将化石時間資訊整合到一個系統發育樹中，進而計算得到Divergence time Tree。
進化速率（Evolutionary rate）指分子的進化速率（蛋白質或核酸等大分子中的氨基酸或核苷酸在一定時間内的替換率），也可指某一類群的物種分化速率，甚至可指特定性狀的進化速率。是以，籠統地說，進化速率指在進化過程中，機關時間内發生的改變。分子鐘的研究與分子進化速率相關。
分子進化速率有針對病毒進化速率研究或非編碼序列的nucleotide substitution速率，也有針對編碼序列的密碼子同義替換速率與非同義替換速率。
一、基本概念
1.分子鐘概念
分子鐘的概念源于對蛋白質的研究。Zuckerkandl and Pauling (1962)比較幾種動物的血紅蛋白、細胞色素C的序列後注意到：這些蛋白質的氨基酸取代速率在不同的種系間大緻相同，即分子水準的進化存在恒速現象。Zuckerkandl and Pauling (1965)提出将1962年提出的概念命名為分子鐘，分子水準進化存在一個“時鐘“， 也即進化速率是近似恒定的。是以，分子鐘假設的成立條件是對于給定的任意大分子（DNA序列或者蛋白質序列）在所有進化譜系中有近似恒定的進化速率，如在一個進化分支上所聚集的突變數與該分支的獨立進化時間長度成正比，其替代速率在進化過程中近似保持一個恒定的數值。
2.分子鐘應用面臨的問題
分子鐘假說成立的條件DNA或者蛋白質序列的替代速率是恒定的。20世紀80年代以來，随着DNA序列資料快速積累，大量的證據表明：在長期進化過程中，很多類群的絕大多數基因或蛋白質的序列替換速率根本不符合分子鐘假說。對于蛋白質序列，在物種适應輻射過程中，其進化速度可能會大大加快。是以，以蛋白質為基礎的恒定進化速率并非理想的分子鐘；對于核酸分子，不同基因的分子鐘速率不同；并且同一基因在不同的生物類群間可能有顯著差異，是以同一基因的分子異速進化現象是顯而易見的。目前分子鐘面臨的一些挑戰也主要與分子異速進化相關，由于分子異速進化的存在，特别是同一級基因在不同生物類群的進化速率可能有顯著的差異，這給應用同一分子鐘來重建物種系統發育關系及估算物種分歧年代帶來了困難，這是分子鐘在應用上面臨的有一個挑戰。
3.分子鐘的類型
在BEAST2中，分子鐘模型的模型可選項有四個包括：-Strict Clock(嚴格分子鐘); -Relaxed Clock Exponential(松散指數分子鐘); -Relaxed Clock Log Normal(松散對數分子鐘); -Random Local Clock(随機本地分子鐘)。
4.分子鐘校正
分子進化是生物分子層次上的進化，分子系統學是從生物大分子（蛋白質、核算）的資訊推斷生物進化曆史，或者說重建系統發育關系，并以系統樹的形式表示出來。系統發生樹的枝長僅代表每位點的堿基替代數或者遺傳距離。基于生物遺傳距離所建立的分子進化速率，使得分子鐘方法在分子系統學中的應用使得利用片段估算物種的起源和分化時間成為可能，并可用于進一步推測不同生物類群在進化曆史上的分歧時間。在速率恒定的假設下，遺傳距離是時間的線性函數，為了将遺傳距離轉化為分歧時間，至少需要一個能夠提供時間資訊的标定點（calibration point）。常用的校準資訊可以分為：(1)已知的堿基替代速率；(2)化石校準點；(3)生物地理事件校準點；(4)二次校準點
5.中性理論
受達爾文自然選擇學說影響，當分子時鐘剛提出的時候，人們認為決定氨基酸替換或者堿基替換的關鍵性因素是自然選擇。如果分子時鐘成立，那麼也就意味着适應性替換在物種的發生是恒定的。但是适應性變異的發生和固定受到很多因素的影響，比如選擇壓力大小、群體效應大小、突變率等等，這些因素在各個物種中差異很大。但是這些因素作用下的适應性變異為什麼會是很定的呢？為此，Kimura在1968年提出了另一個理論，認為如果一個突變沒有帶來表型差異，那麼它在該群體中的命運則完全受各種随機因素的影響。也就是”中性理論“。
二、mcmctree估算物種分歧時間
1.第一步建構物種樹
使用單拷貝直系同源基因四倍簡并位點或者線粒體基因，通過IQTREE或者Mrbayes等軟體建構物種的系統發育樹。這裡就不詳細講建樹流程了。
2. 擷取化石分歧時間
這一步還是很關鍵的，隻有根據标定的化石點才可以估算出樹中每個物種的分化時間。可以根據已發表的文章獲得化石點，可以标定2-3個化石點，也可以5-10個化石點，根據物種樹的規模來看。
假如在已發表的文章中沒有找到可用的化石點，那麼這裡有一個很好用的物種分化時間查詢工具 - TimeTree
TimeTree通過彙總大量的原始文獻，包括大量的人工核實，估計出不同物種的最近共同祖先，并建構出整個所有可及物種的進化樹。目前TimeTree來到了第五個版本，包含了4000多個相關的研究，以及近14萬個物種。通過網頁版或者手機APP版可以友善查詢任意兩個物種的分化時間，或者多個物種的進化樹。我們以黑腹果蠅Drosophila melanogaster, 甘比亞按蚊Anopheles gambiae為例。
網站連結：http://www.timetree.org/
（1）單個物種的進化過程
我們在網頁中輸入黑腹果蠅：
然後我們可以看到整個黑腹果蠅的進化曆程，包括細胞的形成，細胞核的形成……雙翅目的形成……果蠅亞組的形成，從數十億年前到今天的整個過程。
有意思的是，整個時間線不僅包括了物種的進化和形成，還包括了非生物因素的變化，比如日照強度，二氧化碳水準，氧氣水準。同時，還有一些地球（宇宙）事件的影響，比如彗星的撞擊等。比如點選65mya左右的大紅圈，會顯示墨西哥的彗星撞擊。
（2）兩個物種的分化時間
在如下網頁中輸入黑腹果蠅和甘比亞按蚊。
根據12個不同的研究，估計出黑腹果蠅和甘比亞按蚊的中位分化時間大約在2.43億年前， 95%CI（1.96, 2.80）。
（3）多個物種的進化樹
可以輸入某個物種群組，把該群組中的所有物種的關系都顯示出來。或者上傳一個自定義的物種清單。比如我們此處随意輸入4個物種：Drosophila melanogaster （黑腹果蠅） Anopheles gambiae （甘比亞按蚊）Aedes aegypti （埃及伊蚊） ulex pipiens （淡色庫蚊）
3.估算物種分歧時間
step 1: 估計替換速率
輸入檔案：newick格式的系統發育樹（不帶分支長度）、多序列比對檔案、baseml控制檔案
13 1
(((((((Ba,((Hm,(Mc,Kin)),Dp)),(Cc,Cn)),Pr),La),Pp),(Bm,Ha)),Dm)'@2.72'; seqfile = allSingleCopyOrthologsAlign.4Dsite.fas
treefile = species.tree0
outfile = mlb * main result file
noisy = 3 * 0,1,2,3: how much rubbish on the screen
verbose = 1 * 1: detailed output, 0: concise output
runmode = 0 * 0:user tree; 1:semi-automatic; 2:automatic
* 3:StepwiseAddition; (4,5):PerturbationNNI
model = 7 * 0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85, 5:T92, 6:TN93, 7:REV
* 8:UNREST, 9:REVu; 10:UNRESTu
Mgene = 0 * 0:rates, 1:separate; 2:diff pi, 3:diff kappa, 4:all diff
clock = 1 * 0:no clock, 1:clock; 2:local clock; 3:CombinedAnalysis
fix_kappa = 0 * 0: estimate kappa; 1: fix kappa at value below
kappa = 2 * initial or fixed kappas
fix_alpha = 0 * 0: estimate alpha; 1: fix alpha at value below
alpha = 0.5 * initial or fixed alpha, 0:infinity (constant rate)
ncatG = 5 * # of categories in the dG, AdG, or nparK models of rates
fix_rho = 1 * 0: estimate rho; 1: fix rho at value below
rho = 0 * initial or fixed rho, 0:no correlation
Malpha = 0 * 1: different alpha's for genes, 0: one alpha
nparK = 0 * rate-class models. 1:rK, 2:rK&fK, 3:rK&MK(1/K), 4:rK&MK
getSE = 1 * 0: don't want SEs of estimates, 1: want SEs
RateAncestor = 0 * (0,1,2): rates (alpha>0) or ancestral states
method = 0 * Optimization method 0: simultaneous; 1: one branch a time
Small_Diff = 0.5e-6
cleandata = 0 * remove sites with ambiguity data (1:yes, 0:no)?
fix_blength = 0 * 0: ignore, -1: random, 1: initial, 2: fixed
運作
baseml baseml.ctl 輸出檔案中，會得到替換速率，用于mcmctree運算。
step 2: 估計Gradient and Hessian of Branch Lengths
輸入檔案：newick格式的系統發育樹（不帶分支長度，帶有多個化石校正時間）、多序列比對檔案、mcmctree控制檔案
13 1
(((((((Ba,((Hm,(Mc,Kin))'B(58,84)',Dp)),(Cc,Cn)),Pr)'B(79,118)',La),Pp)'B(76,146)',(Bm,Ha)'B(100,122)'),Dm)'B(217,314)'; seed = -1
seqfile = allSingleCopyOrthologsAlign.4Dsite.fas
treefile = species.tree1
outfile = out_usedata3
ndata = 1
usedata = 3 * 0: no data; 1:seq like; 2:normal approximation
clock = 3 * 1: global clock; 2: independent rates; 3: correlated rates
* RootAge = '<10' * constraint on root age, used if no fossil for root.
model = 7 * 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 = 1 2.288 * gamma prior for rate for genes ### 1/替換率
sigma2_gamma = 1 4.5 * gamma prior for sigma^2 (for clock=2 or 3)
finetune = 1: 0.06 0.5 0.006 0.12 0.4 * times, rates, mixing, paras, RateParas
print = 1
burnin = 500000
sampfreq = 5000
nsample = 20000
*** Note: Make your window wider (100 columns) when running this program.
運作
mcmctree mcmctree3.ctl step 3: 估算物種分歧時間
輸入檔案：newick格式的系統發育樹（不帶分支長度，帶有多個化石校正時間）、多序列比對檔案、mcmctree控制檔案、in.BV（由上一步分析獲得）
mv out.BV in.BV seed = -1
seqfile = allSingleCopyOrthologsAlign.4Dsite.fas
treefile = species.tree1
outfile = out_usedata2
ndata = 1
usedata = 2 * 0: no data; 1:seq like; 2:normal approximation
clock = 3 * 1: global clock; 2: independent rates; 3: correlated rates
* RootAge = '<10' * constraint on root age, used if no fossil for root.
model = 7 * 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 = 1 2.288 * gamma prior for rate for genes ### 1/替換率
sigma2_gamma = 1 4.5 * gamma prior for sigma^2 (for clock=2 or 3)
finetune = 1: 0.06 0.5 0.006 0.12 0.4 * times, rates, mixing, paras, RateParas
print = 1
burnin = 500000
sampfreq = 5000
nsample = 20000
*** Note: Make your window wider (100 columns) when running this program. 運作
nohup mcmctree mcmctree2.ctl >mcmctree.ctl.log 2>&1 &; 檢測運作結果
最直接的檢測方法是：分别使用不同的Seed值進行mcmctree或infinitesites進行兩次或多次分析，然後比較兩個結果樹是否一緻，實際就是比較樹檔案中各内部節點的Height值（分歧時間 / Posterior time）。計算各枝長總的偏差百分比，當偏差百分比低于0.1%，則認為兩次結果非常吻合，差異低于0.1%，認為達到收斂。此外，還可以使用Tracer分析mcmc.txt檔案，檢測其ESS值，一般認為該值高于200，則可能達到收斂。該方法可用于輔助檢測。最後，若不收斂，則需要提高burnin、nsample值，重新運作程式。
三、MEGA-X估算物種分歧時間
MEGA-X中我們可以通過「分子鐘」來追溯生物進化的曆程,接下來我們來試一下用MEGA-X的RelTime方法來推斷一下進化時間，首先在推斷前需要準備的檔案有以下幾個：可用于推斷分化時間的系統發育樹（注意這裡的發育樹上需要有外群以及校準點）的序列檔案以及Newick格式的樹檔案。
第一步 ：打開MEGA-X ------ CLOCKS ------ Compute Timetree ------ RelTime-ML
第二步 ：先設定前三個輸入（1）序列檔案fasta格式 (2) 物種樹檔案 .nwk或者.tre格式 （3）設定外群
第三步 ：标定化石點
在CALIBRATE NODES這，點Add Constraints就會進入下面的界面
這裡可以添加多個化石點，根據Timetree查詢到的化石點就可以用，添加好了後點ok就行。
第四步 ：設定可選參數
根據ML樹計算分歧時間，一般預設參數就行，直接點選OK
第五步 ：點選分析
這裡就會得到一個帶有分歧時間的樹，也可以顯示每個節點的時間。
四、分子進化速率計算
進化速率，實際存在兩種不同的東西
第一種就是對非編碼區域來說，核酸的替代的速率。當然編碼區域也可以用核酸替代速率來表示進化速率。這一種算法相對簡單，機關就是nucleotide substitution per Myr。在病毒進化速率計算，非編碼序列的進化速率計算時通常使用這種方法。
第二種就是對蛋白質編碼基因來說，一般是密碼子的同義替換率與非同義替換率來表示進化速率，以及通過dN/dS來表示選擇壓力的大小。可以通過PAML 中的codeml來計算。
1.核酸的替代的速率
這個計算起來很簡單，根據這篇文章中的方法來計算
By dividing the pairwise nucleotide distances between Prolaupala and Laupala elements within each cluster by 10 million years (Myr) (twice the divergence time between Laupala from Prolaupala estimated from biogeography) (18), we obtain an estimate of the absolute rate of point substitution. The average over all clusters is 3.8 * 10-3 nucleotide substitution per Myr, which is almost four times lower than the 15 * 10-3 nucleotide substitution per Myr estimated in Drosophila (19).
公式就是 進化速率 =（P-distance） / (twice the divergence time)
那麼具體怎麼計算呢，不知道有沒有專門的軟體，這個可以自己寫個代碼算。
你先要獲得一個P-distance的矩陣，表示兩兩物種間的分歧距離。MEGA就可以得到，多序列比對完有個計算。如下
可以把這個p-distance矩陣用python讀到資料框中
dt = pd.read_excel("pdistance.xls",sheet_name=[0],index_col=0) 然後需要一個兩兩物種分歧時間的矩陣也讀到資料框中
dt1 = pd.read_excel(input1,sheet_name=[0],index_col=0) 這個是兩兩物種分歧時間機關：百萬年
這個是兩兩物種間分歧距離也是兩兩物種間的堿基替換率。
接下來下寫了個python腳本根據上面進化速率的公式計算即可
output_file = open(output, "w")
k = 0
f = 0
allrate = 0
for i in list:
k = k + 1
list1 = list[k:]
for z in list1:
linename = i + '_' + z
Mya =df.loc[z, i]
if math.isnan(Mya):
Mya = df.loc[i, z]
Pdistance = df1.loc[z, i]
if math.isnan(Pdistance):
Pdistance = df1.loc[i, z]
print( str(Mya) + '_' + str(Pdistance))
rate = float(Pdistance) / (float(Mya)*2)
output_file.write(linename + "\t" + str(rate) + '\n')
f = f + 1
allrate = allrate + rate
meanrate = float(allrate) / f
output_file.write('mean' + "\t" + str(meanrate) + '\n')
output_file.close() 最後就會得到一個平均進化速率。
2.密碼子同義替換與非同義替換速率
使用PAML 中baseml來計算
第一步安裝
wget http://abacus.gene.ucl.ac.uk/software/paml4.9j.tgz
tar xf paml4.9j.tgz
rm bin/*.exe
cd src
make -f Makefile
rm *.o
mv baseml basemlg codeml pamp evolver yn00 chi2 ../bin 第二步 準備輸入檔案
paml的輸入檔案有三個，seq檔案，tree檔案，另外一個是配置檔案，裡面就是調整我們輸入資料，輸出資料，使用模型，進化速率等等的調整。
第三步 修改配置檔案baseml.ctl
seqfile = stewart.aa * sequence data filename
treefile = stewart.trees * tree structure file name
outfile = mlc * main result file name
noisy = 9 * 0,1,2,3,9: how much rubbish on the screen
verbose = 1 * 0: concise; 1: detailed, 2: too much
runmode = 0 * 0: user tree; 1: semi-automatic; 2: automatic
* 3: StepwiseAddition; (4,5):PerturbationNNI; -2: pairwise
seqtype = 2 * 1:codons; 2:AAs; 3:codons-->AAs
CodonFreq = 2 * 0:1/61 each, 1:F1X4, 2:F3X4, 3:codon table
* ndata = 10
clock = 0 * 0:no clock, 1:clock; 2:local clock; 3:CombinedAnalysis
aaDist = 0 * 0:equal, +:geometric; -:linear, 1-6:G1974,Miyata,c,p,v,a
aaRatefile = dat/jones.dat * only used for aa seqs with model=empirical(_F)
* dayhoff.dat, jones.dat, wag.dat, mtmam.dat, or your own
model = 2
* models for codons:
* 0:one, 1:b, 2:2 or more dN/dS ratios for branches
* models for AAs or codon-translated AAs:
* 0:poisson, 1:proportional, 2:Empirical, 3:Empirical+F
* 6:FromCodon, 7:AAClasses, 8:REVaa_0, 9:REVaa(nr=189)
NSsites = 0 * 0:one w;1:neutral;2:selection; 3:discrete;4:freqs;
* 5:gamma;6:2gamma;7:beta;8:beta&w;9:betaγ
* 10:beta&gamma+1; 11:beta&normal>1; 12:0&2normal>1;
* 13:3normal>0
icode = 0 * 0:universal code; 1:mammalian mt; 2-10:see below
Mgene = 0
* codon: 0:rates, 1:separate; 2:diff pi, 3:diff kapa, 4:all diff
* AA: 0:rates, 1:separate
fix_kappa = 0 * 1: kappa fixed, 0: kappa to be estimated
kappa = 2 * initial or fixed kappa
fix_omega = 0 * 1: omega or omega_1 fixed, 0: estimate
omega = .4 * initial or fixed omega, for codons or codon-based AAs
fix_alpha = 1 * 0: estimate gamma shape parameter; 1: fix it at alpha
alpha = 0. * initial or fixed alpha, 0:infinity (constant rate)
Malpha = 0 * different alphas for genes
ncatG = 8 * # of categories in dG of NSsites models
getSE = 0 * 0: don't want them, 1: want S.E.s of estimates
RateAncestor = 1 * (0,1,2): rates (alpha>0) or ancestral states (1 or 2)
Small_Diff = .5e-6
cleandata = 1 * remove sites with ambiguity data (1:yes, 0:no)?
* fix_blength = 1 * 0: ignore, -1: random, 1: initial, 2: fixed, 3: proportional
method = 0 * Optimization method 0: simultaneous; 1: one branch a time
* Genetic codes: 0:universal, 1:mammalian mt., 2:yeast mt., 3:mold mt.,
* 4: invertebrate mt., 5: ciliate nuclear, 6: echinoderm mt.,
* 7: euplotid mt., 8: alternative yeast nu. 9: ascidian mt.,
* 10: blepharisma nu.
* These codes correspond to transl_table 1 to 11 of GENEBANK. seqfile = test.codon * sequence data filename
treefile = test.tree * tree structure file name
outfile = test.rlt * main result file name
noisy = 0 * 0,1,2,3,9: how much rubbish on the screen
verbose = 0 * 0: concise; 1: detailed, 2: too much
runmode = -2 * 0: user tree; 1: semi-automatic; 2: automatic
* 3: StepwiseAddition; (4,5):PerturbationNNI; -2: pairwise###此處模式選擇一般是0或者-2,0表示使用你生成的tree，-2表示兩兩比較
seqtype = 1 * 1:codons; 2:AAs; 3:codons-->AAs
CodonFreq = 2 * 0:1/61 each, 1:F1X4, 2:F3X4, 3:codon table
* ndata = 5504 ###此處即為多少配對的數目，就我們上面統計的
clock = 0 * 0:no clock, 1:clock; 2:local clock; 3:CombinedAnalysis
aaDist = 0 * 0:equal, +:geometric; -:linear, 1-6:G1974,Miyata,c,p,v,a
aaRatefile = dat/jones.dat * only used for aa seqs with model=empirical(_F)
* dayhoff.dat, jones.dat, wag.dat, mtmam.dat, or your own
model = 0
* models for codons:
* 0:one, 1:b, 2:2 or more dN/dS ratios for branches
* models for AAs or codon-translated AAs:
* 0:poisson, 1:proportional, 2:Empirical, 3:Empirical+F
* 6:FromCodon, 7:AAClasses, 8:REVaa_0, 9:REVaa(nr=189)
NSsites = 0 * 0:one w;1:neutral;2:selection; 3:discrete;4:freqs;
* 5:gamma;6:2gamma;7:beta;8:beta&w;9:betaγ
* 10:beta&gamma+1; 11:beta&normal>1; 12:0&2normal>1;
* 13:3normal>0
icode = 0 * 0:universal code; 1:mammalian mt; 2-10:see below
Mgene = 0
* codon: 0:rates, 1:separate; 2:diff pi, 3:diff kapa, 4:all diff
* AA: 0:rates, 1:separate
fix_kappa = 0 * 1: kappa fixed, 0: kappa to be estimated
kappa = 2 * initial or fixed kappa
fix_omega = 0 * 1: omega or omega_1 fixed, 0: estimate
omega = .4 * initial or fixed omega, for codons or codon-based AAs
fix_alpha = 1 * 0: estimate gamma shape parameter; 1: fix it at alpha
alpha = 0. * initial or fixed alpha, 0:infinity (constant rate)
Malpha = 0 * different alphas for genes
ncatG = 8 * # of categories in dG of NSsites models
getSE = 0 * 0: don't want them, 1: want S.E.s of estimates
RateAncestor = 1 * (0,1,2): rates (alpha>0) or ancestral states (1 or 2)
Small_Diff = .5e-6
cleandata = 1 * remove sites with ambiguity data (1:yes, 0:no)?
* fix_blength = 1 * 0: ignore, -1: random, 1: initial, 2: fixed, 3: proportional
method = 0 * Optimization method 0: simultaneous; 1: one branch a time
* Genetic codes: 0:universal, 1:mammalian mt., 2:yeast mt., 3:mold mt.,
* 4: invertebrate mt., 5: ciliate nuclear, 6: echinoderm mt.,
* 7: euplotid mt., 8: alternative yeast nu. 9: ascidian mt.,
* 10: blepharisma nu.
* These codes correspond to transl_table 1 to 11 of GENEBANK. 可以參考上面這個配置檔案，修改好了之後就可以運作baseml
###配置完成後，開始運作
baseml baseml.ctl 結果後就可以得到一下結果檔案
其中包含dN,dS,rst這些為結果檔案，根據需要去畫圖。
以上就是物種分歧時間、分子進化速率的所有内容了，包括了兩種方法推斷物種分歧時間（mcmctree 和 MEGA reltime），TimeTree查詢化石節點，計算分子進化速率（核酸替換速率、密碼子同義替換速率與非同義替換速率）。
如果感覺有幫助的話，下面來個硬币吧！
參考連結
- http://www.chenlianfu.com/?p=2974
- https://www.jianshu.com/p/f9e5fe95478d
- http://www.fish-evol.org/mcmctreeExampleVert6/text1Eng.html
- http://abacus.gene.ucl.ac.uk/software/MCMCtreeStepByStepManual.pdf
- https://dosreislab.github.io/2017/10/24/marginal-likelihood-mcmc3r.html
- http://nebc.nerc.ac.uk/bioinformatics/documentation/paml/doc/MCMCtreeDoc.pdf
https://www.jianshu.com/p/b12e058c6597
https://zhuanlan.zhihu.com/p/105159767
https://www.sci666.com.cn/20394.html
以上内容僅為個人總結，如有錯誤之處敬請批評！
聯系我們：
生信藝術家
點選下方即可關注我們！