Hi PAML Users,
Welcome to the PAML discussion group. This site is for posting questions and discussions of the PAML software package. PAML stands for Phylogenetic Analysis by Maximum Likelihood.
----Ziheng
PAML是一個用最大似然法來對DNA和蛋白質序列進行系統(tǒng)發(fā)育分析的軟件包,此軟件由楊子恒教授開發(fā)断医,最新版本為v. 4.9h孤澎。
此筆記主要用于記錄PAML分析中遇到的一些問題窄陡,并長期更新垂睬。
本文已遷移至個人博客
PAML主頁
PAML手冊
PAML FAQs
PAML discussion group
一. PAML可以做什么鹏秋?
- 系統(tǒng)發(fā)育樹的檢驗與比較;
- 復雜替代模型中參數(shù)的估計川慌;
- 對不同假設模型進行似然率假設性檢驗;
- 基于全局或局部分子鐘模型估算分歧時間祠乃;
- 使用核酸梦重,氨基酸和密碼子模型重構祖先序列;
- Monte Carlo模擬生成核酸亮瓷,密碼子和氨基酸序列數(shù)據(jù)集琴拧;
- 估算同義和非同義替代率和檢測蛋白編碼DNA序列中的正選擇;
- 貝葉斯估算物種的分歧時間來合并化石尺度上的不確定性嘱支;
PAML理論基礎原則 蛋白編碼基因(protein coding sequence)的自然選擇壓力水平可以通過dN/dS(ω)值的大小來衡量蚓胸,其中,dS代表同義替換率(synonymous rate)除师,dN代表非同義替換率(non-synonymous rate)沛膳。在沒有受到選擇壓力時,同義替換率和非同義替換率相等汛聚,此時dN/dS = 1锹安;當受到負選擇或凈化選擇壓力時,自然選擇會阻止氨基酸發(fā)生改變倚舀,同義替換率會大于非同義替換率叹哭,即dN/dS < 1;當受到正選擇壓力時,氨基酸的置換率會受自然選擇的青睞痕貌,即蛋白功能可能會發(fā)生適應性改變风罩,此時dN/dS > 1。
PAML簡易流程
二. 如何運行PAML
要正常運行paml軟件舵稠,需要四個標準文件:
- tree文件:Newick格式的樹文件超升,主要包括分析涉及物種的系統(tǒng)發(fā)生關系,物種名必須與序列比對里面的物種名保持一致柱查,否則會報錯廓俭。
- fasta文件:fasta格式的多序列比對文件,通過要經過mafft唉工、prank等軟件進行比對研乒,在經過gblocks、trimal等質量過濾軟件處理后獲得淋硝,序列長度必須為3的倍數(shù)且不允許存在終止密碼子雹熬。
- codeml.ctl:paml軟件的參數(shù)控制文件宽菜,通過修改該文件參數(shù),實現(xiàn)模型切換
-
codeml.exe:該軟件的可執(zhí)行文件
多序列比對與質量過濾在正選擇分析以及系統(tǒng)發(fā)育研究中尤其重要竿报,關乎到所有分析結果的正確與否铅乡。如何正確的處理編碼序列,在之前的博文中已經詳細說明烈菌,此處不在贅述阵幸。這里,主要強調一下控制文件的設置:
*輸入輸出參數(shù):
seqfile = input.txt *設置輸入的多序列比對文件路徑芽世。
treefile = input.trees *設置輸入的樹文件路徑挚赊。
outfile = mlc *設置輸出文件路徑。
noisy = 9 *設置輸出到屏幕上的信息量等級:0济瓢,1荠割,2,3旺矾,9蔑鹦。
verbose = 0 *設置輸出到結果文件中的信息量等級:0,精簡模式結果(推薦)箕宙;1嚎朽,輸出詳細信息,包含堿基序列扒吁;2火鼻,輸出更多信息。
getSE = 0 *設置是否計算并獲得各參數(shù)的標準誤:0雕崩,不需要魁索;1,需要盼铁。
RateAncestor = 1 *設置是否計算序列中每個位點的替換率:0粗蔚,不需要;1饶火,需要鹏控。當設置成1時,會在結果文件rates中給出各個位點的堿基替換率肤寝;同時也會進行祖先序列的重構当辐,在結果文件rst中有體現(xiàn)。
*數(shù)據(jù)使用說明參數(shù):
runmode = 0 *設置程序獲取進化樹拓撲結構的方法:0鲤看,直接從輸入文件中獲取進化樹拓撲結構(推薦)缘揪;1,程序以輸入文件中的多分叉樹作為起始樹,并使用星狀分解法搜索最佳樹找筝;2蹈垢,程序直接以星狀樹作為起始樹,并使用星狀分解法搜索最佳樹袖裕;3曹抬,使用逐步添加法搜索最佳樹;4急鳄,使用簡約法獲取起始樹谤民,再進行鄰近分枝交換尋找最優(yōu)樹;5攒岛,以輸入文件中的樹作為起始樹赖临,再進行鄰近分枝交換尋找最優(yōu)樹;-2灾锯,對密碼子序列進行兩兩比較并使用ML方法計算DnDs,或對蛋白序列進行兩兩比較計算ML距離嗅榕,而不進行其它參數(shù)(枝長和omega等)的計算顺饮。
* fix_blength = 0 *設置程序處理輸入樹文件中枝長數(shù)據(jù)的方法:0,忽略輸入樹文件中的枝長信息凌那;-1兼雄,使用一個隨機起始進行進行計算;1帽蝶,以輸入的枝長信息作為初始值進行ML迭代分析赦肋,此時需要注意輸入的枝長信息是堿基替換率,而不是堿基替換數(shù)励稳;2佃乘,不需要使用ML迭代計算枝長,直接使用輸入的枝長信息驹尼,需要注意趣避,若枝長信息和序列信息不吻合可能導致程序崩潰;3新翎,讓ML計算出的枝長和輸入的枝長呈正比伞芹。
seqtype = 1 *設置輸入的多序列比對數(shù)據(jù)的類型:1笋敞,密碼子數(shù)據(jù);2,氨基酸數(shù)據(jù)纸镊;3,輸入數(shù)據(jù)雖然為密碼子序列抗俄,但先轉換為氨基酸序列后再進行分析裹刮。
CodonFreq = 2 *設置密碼子頻率的計算算法:0,表示除三種終止密碼子頻率為0,其余密碼值頻率全部為1/61枫攀;1括饶,程序分別計算三個密碼子位點的四種堿基頻率,再算四種堿基頻率在三個位點上的算術平均值来涨,計算61種密碼子頻率時图焰,使用該均值進行計算,這61種密碼子頻率之和不等于1蹦掐,然后對其數(shù)據(jù)標準化技羔,使其和為1,得到所有密碼子的頻率卧抗;2藤滥,程序分別計算三個密碼子位點的四種堿基頻率,計算61種密碼子頻率時社裆,使用相應位點的堿基頻率進行計算拙绊,這61種密碼子頻率之和不等于1,然后對其數(shù)據(jù)標準化泳秀,使其和為1标沪,得到所有密碼子的頻率;3嗜傅,直接使用觀測到的各密碼子的總的頻數(shù)/所有密碼值的總數(shù)金句,得到所有密碼子的頻率。一般選擇第三種方法吕嘀,設置CodonFre的值為2违寞。第一種方法讓所有密碼子頻率相等,不符合實際輕卡偶房;第二種方法沒有考慮到三種位點各堿基頻率的差異趁曼,得到的密碼子頻率不準確;第三種方法能比較準確計算出各密碼子的頻率蝴悉;第四種方法由輸入數(shù)據(jù)得到真實的各密碼子頻率彰阴,但有時候有些非終止密碼子的頻率為0,可能不利于后續(xù)的計算拍冠。
cleandata = 1 *設置是否移除不明確字符(N尿这、?庆杜、W射众、R和Y等)或含以后gap的列后再進行數(shù)據(jù)分析:0,不需要晃财,但在序列兩兩比較的時候叨橱,還是會去除后進行比較典蜕;1,需要罗洗,**有可能會導致運行失敗**愉舔。
* ndata = 1 *設置輸入的多序列比對的數(shù)據(jù)個數(shù)。
clock = 0 *設置進化樹中各分支上的變異速率是否一致伙菜,服從分子鐘理論:0轩缤,變異速率不一致,不服從分子鐘理論(當輸入數(shù)據(jù)中物種相差較遠時贩绕,各分支變異速率不一致)火的;1,所有分枝具有相同的變異速率淑倾,全局上服從分子鐘理論馏鹤;2,進化樹局部符和分子鐘理論娇哆,程序認為除了指定的分支具有不同的進化速率湃累,其它分枝都具有相同的進化速率,這要求輸入的進化樹信息中使用#來指定分支碍讨;3脱茉,對多基因數(shù)據(jù)進行聯(lián)合分析。
Mgene = 0 *設置是否有多個基因的多序列比對信息輸入垄开,以及多各基因之間的參數(shù)是否一致:0,輸入的多序列比對文件中僅包含一個基因時税肪,或多個基因具有相同的Kappa和Pi參數(shù)溉躲;1,輸入文件包含多個基因益兄,這些基因之間是相互獨立的(這些基因之間具有不同的Kappa和Pi值锻梳,且其進化樹的枝長也不相關);2净捅,輸入文件包含多個基因疑枯,這些基因具有相同的Kappa值,不同的Pi值蛔六;3荆永,輸入文件包含多個基因,這些基因具有相同的Pi值国章,不同的Kappa值具钥;4,輸入文件包含多個基因液兽,這些基因具有不同的Kappa值和不同的Pi值骂删。當值是2、3或4時,多個基因的進化樹枝長雖然長度不一樣宁玫,但是呈正比關系粗恢。這點和參數(shù)值等于1是不一樣的。
icode = 0 *設置遺傳密碼欧瘪。其值1-10和NCBI的1-11遺傳密碼規(guī)則對應:0眷射,表示通用的遺傳密碼。
Small_Diff = .5e-6 *設置一個很小的值恋追,一般位于1e-8到1e-5之間凭迹。推薦檢測設置不同的值在比較其結果,需要該參數(shù)值對結果沒什么影響苦囱。
*位點替換模型參數(shù):
model = 1 *若輸入數(shù)據(jù)是密碼子序列嗅绸,該參數(shù)用于設置branch models,即進化樹各分枝的omega值的分布:0撕彤,進化樹上所有分枝的omega值一致鱼鸠;1,對每個分枝單獨進行omega計算羹铅;2蚀狰,設置多類omega值,根據(jù)樹文件中對分枝的編號信息來確定類別职员,具有相同編號的分枝具有相同的omega值麻蹋,沒有編號的分枝具有相同的omega值,程序分別計算各編號和沒有編號的omega值焊切。
*若輸入數(shù)據(jù)是蛋白序列扮授,或數(shù)據(jù)是密碼子序列且seqtype值是3時,該參數(shù)用于設置氨基酸替換模型:0专肪,Poisson刹勃;1,氨基酸替換率和氨基酸的觀測頻率成正比嚎尤;2荔仁,從aaRatefile參數(shù)指定的文件路徑中讀取氨基酸替換率信息,這些信息是根據(jù)經驗獲得芽死,并記錄到后綴為.dat的配置文件中乏梁。這些經驗模型(Empirical Models)文件位于PAML軟件安裝目錄中的dat目錄下,例如收奔,Dayhoff(dayhoff.dat)掌呜、WAG(wag.dat)、LG(lg.dat)坪哄、mtMAN(mtman.dat)和mtREV24(mtREV24.dat)等质蕉;
aaRatefile = dat/wag.dat *當對蛋白數(shù)據(jù)進行分析势篡,且model = 2時,該參數(shù)生效模暗,用于設置氨基酸替換模型禁悠。
aaDist = 0 *設置氨基酸之間的距離。
NSsites = 0 *輸入數(shù)據(jù)時密碼子序列時生效兑宇,用于設置site model碍侦,即序列各位點的omega值的分布:0,所有位點具有相同的omega值隶糕;1瓷产,各位點上的omega值小于1或等于1(服從中性進化neutral);2枚驻,各位點上的omega值小于1濒旦、等于1或大于1(選擇性進化selection);3再登,discrete尔邓;4,freq锉矢;5:gamma梯嗽;6,2gamma沽损;7灯节,beta;8绵估,beta&w显晶;9,beta&gamma壹士;10,beta&gamma+1偿警;11躏救,beta&normal>1;12螟蒸,0&2normal>1盒使;13,3normal>0七嫌。
*可以一次輸入多個模型進行計算并比較少办,其結果輸出的rst文件中。
fix_alpha = 1 *序列中不同的位點具有不同的堿基替換率诵原,服從discrete-gamma分布英妓,該模型通過alpha(shape參數(shù))和ncatG參數(shù)控制挽放。該參數(shù)設置是否給定一個alpha值:0,使用ML方法對alpha值進行計算蔓纠;1辑畦,使用下一個參數(shù)設置一個固定的alpha值。
*對于密碼子序列腿倚,當NSsites參數(shù)值不為0或model不為0時纯出,推薦設置fix_alpha = 1且alpha = 0,即不設置alpha值敷燎,認為位點間的變異速率一致暂筝,否則程序報錯。若設置了alpha值硬贯,則程序認為不同密碼子位點的變異速率不均勻焕襟,且同時所有位點的omega值一致,當然各分枝的omega值也會一致澄成,這時要求NSsites和model參數(shù)值都設置為0(這一般不是我們需要的分析胧洒,它不能進行正選擇分析了)。
alpha = 0 *設置一個固定的alpha值或初始alpha值(推薦設置為0.5)墨状。該值小于1卫漫,表示只有少數(shù)熱點位置的替換率較高;該值越小肾砂,表示位點替換率在各位點上越不均勻列赎;若設置fix_alpha = 1且alpha = 0則表示所有位點的替換率是恒定一致的。
Malpha = 0 *當輸入的多序列比對結果中有多基因時镐确,設置這些基因間的alpha值是否相等:0包吝,分別對每個基因單獨計算alpha值;1源葫,所有基因的alpha值保持一致诗越。
ncatG = 8 *序列中不同位點的變異速率服從GAMMA分布的,ncatG是其一個參數(shù)息堂,一般設置為5嚷狞,4,8或10荣堰,且序列條數(shù)越多床未,該值設置越大。
*對于密碼子序列振坚,當NSites設置為3時薇搁,ncatG設置為3;當NSites設置為4時渡八,ncatG設置為5啃洋;當NSites值設置>=5時传货,ncatG值設置為10。
fix_kappa = 0 *設置是否給定一個Kappa值:0裂允,通過ML迭代來估算Kappa值损离;1,使用下一個參數(shù)設置一個固定的Kappa值绝编。
kappa = 2 *設置一個固定的Kappa值僻澎,或一個初始的Kappa值。
fix_omega = 0 *設置是否給定一個omega值:0十饥,通過ML迭代來估算Kappa值窟勃;1,使用下一個參數(shù)設置一個固定的omega值逗堵。
omega = .4 *設置一個固定的omega值秉氧,或一個初始的omega值。
method = 0 *設置評估枝長的ML迭代算法:0蜒秤,使用PAML的老算法同時計算所有枝長汁咏,在clock = 0下有效;1作媚,PAML新加入的算法攘滩,一次對一個枝長進行計算,該算法僅在clock參數(shù)值為1纸泡,2或3下工作漂问。
樹文件設置:
(((1,2),3),4,5);
((((1,2),3),4),5);
(((1:0.1, 2:0.2):0.8, 3:0.3):0.7, 4:0.4, 5:0.5);
(((Human:0.1, Chimpanzee:0.2):0.8, Gorilla:0.3):0.7,Orangutan:0.4, Gibbon:0.5);
需要注意的是,Newick格式的樹尾部一定要有分號女揭,沒有的話程序無法正常運行蚤假。
三.PAML分析四大模型
在paml分析里面,主要涉及四種模型:位點模型(site-model)吧兔、支模型(branch model)磷仰、支位點模型(branch site model)以及進化支模型。位點模型通常適用于檢測某一支系普遍性境蔼、廣泛性的正選擇芒划,這種正選擇是由于位點持續(xù)改變所引起的,例如欧穴,適應多種病原體;支模型主要是檢測某一支系是否存在快速進化泵殴、選擇壓力約束以及正選擇涮帘,但卻無法檢測到正選擇位點;支位點模型比較準確與穩(wěn)定笑诅,適用于檢測某一支系斷點性的正選擇事件调缨,此結果是由于適應某一時期環(huán)境改變所引起的,通常會保留在后代中弦叶;進化支模型(clade model)伤哺,則主要是判別不同物種之間是否存受分化選擇壓力作用燕侠,不局限于正選擇,它一次可以標記多個分支進行比較绢彤。
1. 位點模型
M0:假設所有位點具有相同的dN/dS值蜓耻;
M1a:假設存在兩類位點—保守位點0 < dN/dS < 1茫舶,中性進化位點dN/dS = 1,并且估算這兩類位點的比率(p0刹淌,p1)和ω值(ω0,ω1)疹启;
M2a:假設存在三類位點——純化選擇位點dN/dS < 1柠衅,中性進化位點dN/dS = 1與正選擇位點dN/dS > 1菲宴,并估算三類位點的比率(p0,p1势誊,p2)谣蠢;
M3: 離散模型,假設所有位點的ω值呈離散分布挤忙;
M7:假設所有位點0 < ω <1且呈現(xiàn)beta分布谈喳;
M8:在M7模型的基礎上婿禽,增加一類正選擇位點(ω >1 );
M8a:與M8類似大猛,只不過是將新增的ω固定為1淀零;
位點模型的Codeml.ctl參數(shù)設置:
model = 0
NSsites = 0 1 2 3 7 8
即分別比較如下模型(* 為正選擇模型,前提是需要LRTs顯著):
M3(分散比率) vs. M0(負選擇)
M2a* (正選擇)vs. M1a (選擇約束放松)
M8*(正選擇) vs. M7 (選擇約束放松)
M8* (正選擇)vs. M8a(選擇約束放松)
2.支模型
one rario:假設所有的進化譜系都具有相同的ω值唉堪;
free ratio:假設所有的支系都具有獨立的ω值哀卫;
two ratio:假設前景支與背景支的ω不同
分支模型的Codeml.ctl參數(shù)設置:
two ratio: model = 2, NSsites = 0
one ratio: model = 0, NSsites = 0
free ratio: model = 1, NSsites = 0
one ratio(負選擇) vs. free ratio(自由比率)
one ratio(負選擇) vs. two ratio (正選擇)
two ratio natural(中性進化) vs. two ratio (正選擇)
3.支位點模型
假定位點間的ω值是變化的此改,同時也假定支系間的ω值是變化的共啃。該模型主要用于檢測前景支中正選擇作用對部分位點的影響。
modelA null(零假設):ω值設定為固定值1
modelA(備擇假設):估算其ω值是否大于1
背景支與前景支具有相同的位點ω值:
K0:前景支與背景支中的位點受到純化選擇0 < ω < 1究珊;
K1:前景支與背景支中的位點處于中性進化0 < ω = 1;
背景支與前景支具有不同的位點ω值:
K2a:前景支處于中性進化剿涮,而背景支處于純化選擇攻人;
K2b:前景支受到正選擇壓力( ω>1 ),而背景支處于中性進化瞬浓;
支位點模型的的Codeml.ctl參數(shù)設置:
ModelA:
model = 2, NSsites = 2, fix_omega = 0, omega = 1.5
ModelAnull:
model = 2, NSsites = 2, fix_omega = 1, omega = 1
即比較如下模型(* 為正選擇模型猿棉,前提是需要LRTs顯著):
ModelA*(正選擇) vs. ModelAnull (中性進化)
4.進化支模型
與枝位點模型類型屑咳,能同時檢測多個進化枝(Clade)兆龙,但是該模型并沒有將背景支的dN/dS值約束在(0,1)。
M2a_model: model = 0 NSsites = 22
Clade Model C: model = 3 NSsites = 2
Clade Model D: model = 3 NSsites = 3
Clade Model C vs. M2a_model
四.PAML FAQ
1. 有根樹與無根樹
簡單的區(qū)別方法就是,有根樹的祖先節(jié)點是二叉樹坝橡,而無根樹為三叉樹计寇,舉個例子:
((1,2),3); # 有根樹
(1,2,3); #無根樹
(((1,2),3),4); # 有根樹
((1,2),3,4); #無根樹
(1,2,(3,4)); #無根樹
在使用codeml時,如果沒有指定有根樹參數(shù)卻使用了有根樹作為輸入元莫,那么在輸出結果中將會得到這樣的報錯信息:"This is a rooted tree. Please check!"
踱蠢。對于大多數(shù)模型棋电,即使使用有根樹赶盔,其該模型似然值仍然是正確的,但是root周圍兩個分支的長度不穩(wěn)定撕攒,因為它們的和是估計值烘浦。對于其他模型谎倔,似然估計和參數(shù)估計都是不正確的片习。因此,分析時確實應該注意到這一信息状知,并盡可能使用一棵無根樹孽查。我們可以使用R包ape將有根樹轉換為無根樹:
setwd("C:/Users/lenovo/Desktop/")
library(ape)
tr <- read.tree("tree.txt")
unrooted_tr <- unroot(tr)
write.tree(unrooted_tr, "tree_unrooted.txt")
** 如果你有多個樹需要轉換,可以使用lapply函數(shù)操作**
write.tree(lapply(trees, unroot), "trees_unrooted.tr")
2.是否刪除多序列比對中的gap與模糊字符
在多序列比對過程中瓣铣,對齊gap是極其困難的贷揽,paml軟件包沒有辦法無法處理gap禽绪。因此,我們可以通過設置cleandata = 1
來去除gap循捺;此外从橘,還可以將gap當作為模糊字符進行處理柠衍。但是珍坊,這都不是最好的解決辦法,這兩種策略都低估了序列差異驻民。就個人而言履怯,我認為除了一個或兩個序列之外叹洲,大多數(shù)序列都有序列信息的位點也許應該保留运提,而除了一個或兩個序列之外,所有序列都有對齊間隙的位點最好被移除癣丧。因此胁编,選擇合適的多序列比對軟件以及過濾軟件,尤其重要早直。
3.當前景支的dN/dS值顯著大于背景支時市框,該如何解釋
如果檢測前景支時拾给,其dN/dS > 1時蒋得,我們可以認為它受到正選擇作用乒疏。但是怕吴,如果其dN/dS <1但大于背景支時,就不能認為它是受正選擇作用的伟件,選擇壓力約束放松可能是較為合理的解釋斧账。此外煞肾,在Ohta's
的微有害突變假說下籍救,凈化選擇在大種群中比在小種群中更有效,因此不同譜系的種群規(guī)模的差異提供了另一個相容的假設闪萄。如果氨基酸的變化是稍微有害的桃煎,我們預計在大群體中它們從群體中移除的速度會比在小群體中更高大刊。因此,即使兩系在選擇壓力或基因功能上沒有差異搜锰,我們也期望在一個大群體中看到一個較小的dN/dS比值蛋叼,例如剂陡,許多核基因的dN/dS比值在嚙齒類動物中低于靈長類或偶蹄類鸭栖。
4.不同的模型鑒定出不同的位點,該選擇相信哪一種模型松却?
通常晓锻,如果一個位點在一種模型下出現(xiàn)在列表中飞几,那么在另一種模型下也會有相當大的概率屑墨。如果你以這種方式看待它們绪钥,結果可能沒有太大不同。確定位點的問題很困難匣吊,而且容易出錯色鸳。這種情況類似于從一個班級中得到排名前幾的學生见转。列表包含的內容越多斩箫,質量就越差。因此狐血,我們通常會認為后驗概率大于95%或者99%的位點匈织,是比較可信的缀匕。
5.如何標記前景支
如果你是關注某一個節(jié)點,則可以使用"#"岳链;如果是關注某一類群,則可以使用"$"零远。對于位點模型以及free ratio分析厌蔽,無需標注前景支奴饮;對于支模型以及進化支模型戴卜,單次分析可以標注多個前景支;但支位點模型單次分析僅能標注一個前景支师脂。
(((rabbit, rat) $1, human), goat_cow, marsupial)
(((rabbit #1, rat #1) #1, human), goat_cow, marsupial)
上面的兩個例子其實是等同的吃警,因此$1可以標記大的分支酌心,包括其祖先節(jié)點以及現(xiàn)生節(jié)點挑豌。而#1僅僅是表示末端枝或者祖先節(jié)點。
這里有一些關于嵌套進化枝標簽的規(guī)律泰鸡。符號#的優(yōu)先級高于$盛龄,樹頂端的進化枝標簽優(yōu)先于接近根處的祖先節(jié)點的進化枝標簽芳誓。下面的兩棵樹也是一樣锹淌。第一個樹中$1適用于整個胎盤哺乳動物的進化枝(除人類世系外)赂摆,$2適用于兔與鼠的進化枝烟号。
((((rabbit, rat) $2, human #3), goat_cow) $1, marsupial);
((((rabbit #2, rat #2) #2, human #3) #1, goat_cow #1) #1, marsupial);
使用TreeView軟件可以很方便的創(chuàng)建樹文件并且檢查樹和標簽是否正確。上面作為例子的樹都可以被TreeView識別达传。在TreeView X中宪赶,你必須使用單引號來標記搂妻。如下:
((((rabbit '#2', rat '#2') '#2', human '#3') '#1', goat_cow '#1') '#1', marsupial);
此外叽讳,還可以同時標注多個支為統(tǒng)一前景支坟募,例如
(((a #1, b) c #1 d) e)
6.基因復制后選擇壓力檢測
鏈接1
鏈接2
[圖片上傳失敗...(image-4ab1a5-1584449462345)]
某一基因懈糯,在動物祖先發(fā)生基因復制赚哗,分化為A,B兩個不同支系渐逃。
7. dN/dS代表進化速率而非突變率
較高的dN/dS值可以解釋為正選擇或者快速進化民褂。盡管突變本身也處于選擇壓力之中(大多數(shù)為凈化選擇)赊堪,但不可以解釋為“基因A增加了突變的選擇壓力”哭廉,因為突變是隨機發(fā)生的。原則上講辽幌,突變率會同時影響dN乌企、dS,但通常dN/dS不受突變率影響梁剔。
dN/dS是一種進化速率荣病,但它不是突變速率渗柿,因為同義和非同義替代率具有不同的選擇約束水平朵栖。 選擇壓力測試的基本原理是:它假設同義替換是中性進行陨溅,也就是說门扇,它們大多是在遺傳漂移下進化的偿渡。如果這是真實情況溜宽,那么dS可以用作(中性)突變率的替代适揉。然而涡扼,非同義替代率總是處于凈化選擇壓力吃沪,在正選擇下程度較小票彪。因此不狮,dN/dS是中性偏離的度量摇零。所以驻仅,dN > dS噪服,既dN/dS > 1,則受到正選擇粘优;如果dN小于dS,則dN/dS < 1丹墨,則是凈化選擇带到。選擇壓力測試的關鍵就是它通過特定基因的同義替換的“中性”進化速率歸一化非同義替換的速率揽惹。
8.基因樹or物種樹
在任何情況下搪搏,使用代表真實演化歷史的基因樹都是最好的疯溺。但是囱嫩,有時可能無法輕易判斷是否符合真實演化歷史墨闲,那么你可以選擇物種樹作為代替鸳碧≌袄耄基因組水平分析,那么推薦使用物種樹推励。你可以可以使用基因樹和物種樹進行數(shù)據(jù)穩(wěn)健性測試验辞。
9. 祖先節(jié)點有正選擇信號而整個枝系卻沒有檢測到如何解釋
如果以哺乳動物祖先為前景受神,假設該基因在共同祖先中存在適應性,可能是由于獲得了新的功能财著,但隨后該基因在凈化選擇下得到了保守進化联四。如果你把整個clade作為前景,那么假設該基因在整個哺乳動物中的所有分支都承受著持續(xù)的壓力來改變或多樣化撑教,如果該基因涉及到防御或免疫朝墩,則可能是這種情況。
對于到底是檢測祖先伟姐,還是整個枝系收苏,取決于生物學問題亿卤。例如鹿霸,溶菌酶在所有colobine猴子中都應該具有相同的功能排吴,因此蛋白質在clade內被期望受到選擇性約束,但在colobine clade的分支上懦鼠,該酶顯然獲得了一個新的功能钻哩,其正選擇驅動了氨基酸的變化。有了這個假設肛冶,你就應該把分支的祖先貼上clade的標簽街氢,而不是那些在clade中的分支埋涧。
10.Clade模型判別背景支是否同樣受正選擇
經常遇到審稿人的審稿意見是這樣的:前景分支上正選擇的基因的重要支持并不意味著在背景分支上沒有正選擇界酒,這些基因可能仍在許多(即使不是全部)背景分支上處于正選擇狀態(tài)。為了檢驗原假設(基因僅在前景支受到正選擇)是否正確绒极,可以進一步通過Clade模型檢驗扣泊,因為進化枝模型允許在不限制背景dN/dS小于1的約束下估計前景支與背景支dS/dN的比率近范。
11.基因受到正選擇但無正選擇位點
一種可能性是,對整個基因有積極選擇的證據(jù)延蟹,但每個單獨位點的信息或證據(jù)太弱评矩。 可以查看rst文件,該文件具有所有站點的后驗概率阱飘,以查看是否是這種情況斥杜,mlc文件只列出后驗概率高于0.5的文件。
12. 有正選擇位點但是通過序列比對卻找不到
可能是因為codeml在刪除gap或模糊字符的列沥匈,之后重新編號位點了(cleandata =1)蔗喂。
13. dN/dS(ω)值過大,結果是否可信高帖?
遇到這種極大值的dN/dS時缰儿,比如ω = 999,首先要確保你的序列是否正確散址;其次乖阵,該位置的dn,ds是否遠小于0.0001预麸,枝長是否過小瞪浸。顯然,高度相似的序列和非常發(fā)散的序列都不是信息豐富的吏祸,很難指定確切的值对蒲。為了避免此類問題的出現(xiàn),可以先通過M0模型獲得枝長,再將具有分支長度的進化樹應用到codeml蹈矮,并在ctl中設置FIX_blength=2砰逻。
14. Branch-site檢測適應性趨同進化
如圖所示,紅色分支代表表型趨同的進化分支泛鸟,如果你想通過分支位點模型來檢測適應性趨同進化诱渤,應該將所有紅色分支統(tǒng)一設置為前景分支。當然谈况,前提是假設所有前景支均具有相同的位點受到正選擇勺美。 對于背景支是否同樣存在類似的適應性趨同,可以通過clade進化支模型進行檢測碑韵。
15. 若p2(正選擇位點概率)接近于0
p0/ω0,p1/ω1,p2=(1-p0-p1)/ω2:在正選擇分析的備擇假設結果文件中赡茸,通常會獲得三個p值,其中祝闻,p0代表處于凈化選擇下的位點概率占卧;p1代表處于中性進化下的位點概率;p2則代表處于正選擇下的位點概率联喘。
16.分支模型檢測正選擇)
可以使用two ratio備擇假設(fix_omega = 0 omega = 1)與two ratio零假設(fix_omega = 1 omega = 1)進行統(tǒng)計假設檢驗华蜒。
17. 選擇約束放松
codeml檢測選擇約束放松可以通過2步完成:首先,確定dN/dS顯著增加的情況(由于正選擇或者選擇約束放松)豁遭;然后叭喜,過濾掉顯著正選擇的情況。
18. 不同模型下的起始ω
針對于CladeC和CladeD模型蓖谢,一般要設置幾個不同的起始ω捂蕴,測試其lnL值是否穩(wěn)定(ω=0.001,ω=0.01闪幽,ω=0.1啥辨,ω=0.5,ω=1盯腌,ω=1.5)溉知,最終一般會選取lnL較大的值作為最終值。
19. 枝長會影響PAML結果
正常分析腕够,應該先用M0估計樹的枝長级乍,Kappa值,然后用跑出來的樹作為初始樹燕少,并設置fix_blength = 2卡者。
20. 基因家族選擇壓力分析
21. CladeC檢測正選擇
CladeC經常用于檢測蒿囤,不同分支的分化選擇壓力客们,但有時候檢測到前景枝dN/dS>1,此時我們需要進一步利用CladeC的零假設進一步測試(fix_omega=1,omega=1),或者利用branch-site進一步驗證一致性底挫。
22. Free ratio檢測結果是否可以作為正選擇的證據(jù)
free ratio模型估計通常會造成較大的抽樣誤差恒傻,例如,較短的分支通常會具有較大的dS/dN建邓。所以盈厘,一般對于free ratio出現(xiàn)dN/dS > 999或者dN、dS < 0.0005結果不建議采用官边。
23. dS>5
24. CladeC模型檢測多組支系選擇壓力差異
在dataset中沸手,有三個clade:A, B, C,那么在clade 1和clade 2之間是否存在顯著差異注簿?
假設契吉,前景枝的標注如下所示:
((A1,A2)$1 , (B1,B2)$2 , (C1,C2)$0);
首先,為了測試clades之間的顯著差異诡渴,您可以比較CladeC與M2a_rel捐晶。M2a_rel假設2妄辩、$0等都是在相同的選擇壓力下進化的惑灵,這個測試應該有兩個自由度。
(A1眼耀,a2)$1英支,(b1,b2)$1哮伟,(c1潭辈,c2)$0);
其次,為了測試clade A和B之間的顯著差異澈吨,同時允許clade C是不同的把敢,您可以比較使用上面提供的樹運行CMC與使用更簡單的樹運行CMC的契合程度。在這種情況下谅辣,更簡單的樹會將clade A和B分配給同一個組修赞,這個測試應該有一個自由度。如下所示:
25. 多個前景支的支位點模型檢測
當數(shù)據(jù)集中有多個前景支時:1)進行多個測試桑阶,然后在每個測試中設置一個感興趣的分支作為前景分支柏副;2)只進行一次測試,在此測試中中將所有感興趣的分支設置為前景分支蚣录。那么割择,此時又會出現(xiàn)另一個問題就是進行多個測試時其他感興趣分支是否應該被去除?這也許應該取決于具體的生物學問題萎河。
https://groups.google.com/g/pamlsoftware/c/aVj2opOg7PA
26.p值校正
27. root-to-tip dN/dS
使用free ratio結果可能會產生較大誤差https://groups.google.com/g/pamlsoftware/c/2DrYS0ff7_o
28. paml輸出結果中正選擇位點
正選擇位點的狀態(tài)是多序列比對中第一條參考序列的狀態(tài)蕉饼,而不是前景支的序列狀態(tài)。此外玛歌,還有注意cutdata是否設置為1昧港。
https://groups.google.com/g/pamlsoftware/c/ZnPaysiZKbI