許多基因組的文章里都會提到使用PAML進行基因的正選擇分析(positive selection), 網(wǎng)上也有一些教程介紹如何用PAML進行分析饿幅。無論是文章忱叭，還是教程玫鸟，大多只介紹了過程梢卸，讀完之后，是能夠做相應(yīng)的分析了摔笤，但卻不知道為什么要這樣子做够滑，這篇教程就做這一方面的補充。
我們應(yīng)該知道組成蛋白序列的氨基酸對應(yīng)的核苷酸突變可以分為兩類吕世，同義置換(synonymous substitution)和非同義置換(nonsynonymous substitution),通過計算非同一置換速率和同義置換速率的比值, omega=dN/dS, 我們可以衡量蛋白的選擇壓力彰触。如果選擇不影響物種適應(yīng)環(huán)境，那么比值兩者的速率應(yīng)該相等命辖，因此比值為1况毅，如果非同義突變會降低物種的適應(yīng)性，那么dN < dS, 因此比值小于1尔艇， 如果非同義突變讓提高物種的適應(yīng)性尔许，那么dN > dS, 則比值大于1. 于是乎，非同義突變率顯著高于同義突變即為蛋白質(zhì)適應(yīng)性進化的證據(jù)终娃。
接下來味廊，我們以書中被子植物光敏色素(phy)的適應(yīng)性進化作為例子結(jié)合理論來講解，數(shù)據(jù)方式如下
# PAML的安裝和配置不再贅述
# codeml的配置文件
wget http://abacus.gene.ucl.ac.uk/ziheng/data/phyACF.codeml.ctl
# 15個物種的codon聯(lián)配結(jié)果, 包含gap
wget http://abacus.gene.ucl.ac.uk/ziheng/data/phyACF.txt
# 15個物種的系統(tǒng)發(fā)育樹
wget http://abacus.gene.ucl.ac.uk/ziheng/data/phyACF.trees
下載的配置文件phyACF.codeml.ctl信息顯示如下
seqfile = phyACF.txt
treefile = phyACF.trees
outfile = mlc * main result file name
noisy = 3 * 0,1,2,3,9: how much rubbish on the screen
verbose = 0 * 1: detailed output, 0: concise output
runmode = 0
seqtype = 1 * 1:codons; 2:AAs; 3:codons-->AAs
CodonFreq = 2 * 0:1/61 each, 1:F1X4, 2:F3X4, 3:Fcodon
aaDist = 0 * 0:equal, +:geometric; -:linear, 1-6:G1974,Miyata,c,p,v,a
model = 2
NSsites = 2
icode = 0 * 0:universal code; 1:mammalian mt; 2-10:see below
fix_kappa = 0 * 1: kappa fixed, 0: kappa to be estimated
kappa = 5 * initial or fixed kappa
fix_omega = 0 * 1: omega or omega_1 fixed, 0: estimate
omega = 0.1
fix_alpha = 1 * 0: estimate gamma shape parameter; 1: fix it at alpha
alpha = .0 * initial or fixed alpha, 0:infinity (constant rate)
ncatG = 3 * # of categories in dG of NSsites models
clock = 0 * 0:no clock, 1:global clock; 2:local clock; 3:TipDate
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)?
* ndata = 3
method = 0 * 0: simultaneous; 1: one branch at a time
雖然參數(shù)很多棠耕，但正選擇分析上, 我們所需要修改的參數(shù)是 model, NSsites, fix_omega, omega 這四項余佛，來決定codeml的分析模式。
首先窍荧， 修改配置文件中的 model = 0 和 NSsites = 0. 此時, codeml會計算全局omega.
使用codeml phyACF.codeml.ctl 運行辉巡，輸出結(jié)果在當(dāng)前目錄的mlc文件中。里面信息很多蕊退，我們重點關(guān)注如下幾項
...
lnL(ntime: 27 np: 29): -29984.121043 +0.000000
...
Detailed output identifying parameters
kappa (ts/tv) = 1.98351
omega (dN/dS) = 0.08975
lnL是似然值(likeilhood value)的自然對數(shù)郊楣，之所以是負數(shù)憔恳，是因為計算出似然值是一個非常小的小數(shù)，如果不取對數(shù)痢甘，結(jié)果顯示就是0喇嘱，難以使用。
從結(jié)果來看塞栅，我們算出的全局omega非常小者铜，約為0.09。這很容易理解放椰，因為我們都知道一個蛋白序列作烟，保守的位點肯定遠遠多于不保守的位點，那么平均下來砾医，整體的值就會很小拿撩。因此正選擇通常分析的是系統(tǒng)發(fā)育關(guān)系的特定譜系或者是蛋白質(zhì)的某幾個位點。
接下來如蚜，我們修改 model = 2 和 NSsites = 0, 此時codeml會分析我們提供系統(tǒng)發(fā)育樹中某個分支(foreground)相對于其他分支(background)是否處于正選擇压恒。 問題來了，codeml如何判斷哪個是foreground错邦，哪些是background呢探赫？此時需要看下 phyACF.trees.
1
((C.Sorg:1.414715, (F.Tom:1.174355, C.Arab:1.734907):0.401510):7.949045 #1, (Oat3:0.217161, (A.Rice:0.255094, (A.Zea:0.084488, A.Sorg:0.041302):0.239315):0.038530):1.329584, ((A1.Pea1:0.304507, A.Soy:0.370169):0.329161, ((A.Pars:0.912283, (A.Tob:0.154040, (A.Tom:0.051204, A.Pot:0.054514):0.100673):0.456387):0.182221, (A.Zuc:0.729785, A.Arab:0.792669):0.136483):0.130833):0.547791);
我們可以使用iTOL這個網(wǎng)頁工具對這個樹進行展示，樹形如下撬呢。
不難從圖中發(fā)現(xiàn)伦吠，有一個明顯和其他格格不入的分支，其中包含C. Sorg(高粱,PhyC)和 C.Arab(擬南芥,PhyC), F.Tom(番茄,PhyF)魂拦。我們想要檢驗這一支毛仪，是不是受到了正選擇。 為了強調(diào)這一支芯勘，我們需要在樹結(jié)構(gòu)對應(yīng)的位置上加上#1箱靴，表示foreground, 而這在我們下載的.tree文件中已經(jīng)有了，你可以檢查下荷愕。
我們運行codeml phyACF.codeml.ctl衡怀，關(guān)注輸出文件mlc里的如下內(nèi)容
...
lnL(ntime: 27 np: 30): -29983.513876 +0.000000
...
Detailed output identifying parameters
kappa (ts/tv) = 1.99551
w (dN/dS) for branches: 0.08998 0.03881
在dN/dS這一行有兩個omega值，第一個是background, 第二個是foreground路翻。雖然后者比前者大，但是也沒有超過1茄靠。另外茂契，似然比檢驗(LRT)也顯示，這兩個模式?jīng)]有明顯差異慨绳。
這里掉冶，我們提到了似然比檢驗(likelihood-ratio test, LRT), 這個概念很重要真竖，我們需要稍稍展開說明下。它指的是厌小，根據(jù)兩個競爭的統(tǒng)計模型的似然值的比值恢共，評估兩者的擬合度，其中一個是最大化整個參數(shù)空間璧亚，另一個則是做一些限制讨韭。如果限制條件（零假設(shè)）被觀測數(shù)據(jù)所支持，那么兩者的似然值的差異不會超過抽樣誤差癣蟋。因此透硝，LRT檢驗的是，比值是不是和1有顯著區(qū)別疯搅，或者說比值的自然對數(shù)和0有顯著區(qū)別濒生。
通常似然比檢驗的統(tǒng)計量表現(xiàn)為兩個對數(shù)似然值的差值
這個統(tǒng)計量，我們可以使用卡方檢驗(chi2)來分析它的顯著性幔欧。
# 用R算出
abs(-2*(-29983.513876 - -29984.121043))
# PAML計算chi2, 自由度設(shè)置見最后的補充
chi2 1 1.27
df = 1 prob = 0.270475480 = 2.705e-01
對比我們model=2和model=0, 顯著性0.27 > 0.05, 不足以拒絕零假設(shè)罪治，即我們檢測分支的omega相對于全局的omege沒有明顯差異。
接下來礁蔗，我們來介紹在所有教程都會提到的branch-site model. 也就是設(shè)置model=2 NSsites=2. 它會分析目標(biāo)分支里的位點是不是受到了正選擇觉义。
我們需要建立兩個假設(shè)，分別是零假設(shè)和備擇假設(shè)
- 零假設(shè): 檢驗的分支里的位點不受選擇瘦麸，我們設(shè)置參數(shù)fix_omega=1, omega = 1
- 備擇假設(shè): 檢驗的分支里的位點受到正選擇谁撼，我們設(shè)置fix_omege=0, omege = 1.1
在零假設(shè)時，輸出結(jié)果的內(nèi)容如下滋饲，記錄lnL= -29704.738847
lnL(ntime: 27 np: 31): -29704.738847 +0.000000
...
Detailed output identifying parameters
kappa (ts/tv) = 2.16177
MLEs of dN/dS (w) for site classes (K=4)
site class 0 1 2a 2b
proportion 0.77433 0.07298 0.13953 0.01315
background w 0.07767 1.00000 0.07767 1.00000
foreground w 0.07767 1.00000 1.00000 1.00000
...
在備擇假設(shè)時厉碟，輸出結(jié)果的內(nèi)容如下，記錄lnL= -29694.784206
lnL(ntime: 27 np: 31): -29694.784206 +0.000000
...
Detailed output identifying parameters
kappa (ts/tv) = 2.18201
MLEs of dN/dS (w) for site classes (K=4)
site class 0 1 2a 2b
proportion 0.81323 0.07539 0.10194 0.00945
background w 0.07958 1.00000 0.07958 1.00000
foreground w 0.07958 1.00000 17.59530 17.59530
計算統(tǒng)計量
abs(-2 *(-29694.784206 - -29704.738847 ) )
# 19.90928
卡方檢驗(關(guān)于自由度是1, 見最后的補充)
chi2 1 19.9
df = 1 prob = 0.000008160 = 8.160e-06
p遠遠小于0.01, 為正選擇提供了強有力的證據(jù)屠缭。此時檢查mlc輸出文件的如下內(nèi)容箍鼓，我們還能夠確定被選擇的位點有哪些。 *表示顯著性呵曹。
Bayes Empirical Bayes (BEB) analysis (Yang, Wong & Nielsen 2005. Mol. Biol. Evol. 22:1107-1118)
Positive sites for foreground lineages Prob(w>1):
17 T 0.539
33 G 0.634
34 D 0.672
35 S 0.750
43 E 0.621
55 R 0.982*
61 I 0.643
66 H 0.924
71 K 0.701
102 T 0.974*
104 V 0.932
105 S 0.994**
115 D 0.875
117 P 0.990*
120 G 0.819
130 T 0.976*
....
最后補充下為什么自由度是1, 一方面這是PAML的文檔中提到的款咖。。
另一方面奄喂，根據(jù)定義铐殃，自由度指的是統(tǒng)計量計算時不受限制的變量個數(shù)。在我們檢驗分支的時候跨新，只有一個foreground w是需要自由富腊，所以df=1; 在我們對比branch-sites model時，備擇假設(shè)相對零假設(shè)域帐，只有一個w2是自由(fix_omega = 0), 所以df=1
以上是我在學(xué)習(xí)正選擇分析時的整理結(jié)果赘被，由于數(shù)學(xué)功底太弱是整，有些疑問我還沒有結(jié)果，比如說似然值為什么那么忻窦佟浮入？ML模型參數(shù)是如何確定的？輸出結(jié)果中Bayes Empirical Bayes是怎么運算的羊异？這些還需要不斷的學(xué)習(xí)事秀。
