最近被數(shù)據(jù)分析虐的死去活來(lái)牵现，好難。
參考文獻(xiàn)：
Likelihood ratio tests for detecting positiveselection and application to primate Lysozyme Evolution
User guide of Paml version 4.9
理論部分
dS ：synonymous substitution 同義替換
dN ：non-synonymous substitution 非同義替換
同義替換指：由于密碼子具有簡(jiǎn)并性嗽桩，密碼子中一些核苷酸的替換并不會(huì)導(dǎo)致編碼氨基酸的變化嗅骄。非同義替換指：密碼子中核苷酸的替換導(dǎo)致了編碼氨基酸的變化床绪。
dN/dS：非同義替換和同義替換速率的比值
Darwinian selection may have a non-synonymous/synonymousrate ratio (dN/dS) that is different fromthose of other lineages or greater than one. (positive selection)
基因的變異只有表達(dá)出來(lái)（編碼蛋白發(fā)生改變）机杜，自然選擇才能起作用帜讲。
如果一個(gè)分支的基因變異受到自然選擇的正選擇（positive selection）（正選擇指的是自然選擇支持變異，變異是有利的）椒拗，那么該分支的dN/dS應(yīng)該與其他分支不同似将，或者大于1获黔。
需要注意的是某基因在某一分支的dN/dS大于1來(lái)證明該基因在該分支上受正選擇是很保守的，dN/dS很少會(huì)大于1在验，如果某一個(gè)分支的dN/dS值顯著大于其他分支玷氏，也能說(shuō)明問(wèn)題。
如果分支的某基因選擇上中性（neutral selection）腋舌，即無(wú)害也無(wú)利预茄，那么理論上dN/dS=1
但大多數(shù)情況下，dN/dS<1侦厚，即非同義 替換率小于同義替換率
The basic model for the likelihood analysisis a version of the codon-substitution model of Goldman and Yang (1994) andaccounts for the genetic code structure, transition/transversion rate bias, anddifferent base frequencies at codon positions.
codeml里用于最大似然分析的基礎(chǔ)模型（basal model）基于一個(gè)密碼子替換模型，該模型考慮了遺傳序列的結(jié)構(gòu)拙徽、轉(zhuǎn)換/顛換率的偏差刨沦、以及密碼子不同位置替換的基礎(chǔ)頻率。
從密碼子i到密碼子j的瞬時(shí)替換率為：
K代表轉(zhuǎn)換率/顛換率（transition/transversion）
ω 代表dN/dS
The ω ratio is a measure of natural selection acting on the protein. Simplisitically, values of ω<1, =1, and >1 means negative purifying sleection, neutral selection, and positive selection.
πj 代表密碼子j的平衡頻率膘怕，由密碼子三個(gè)位點(diǎn)的核苷酸頻率估算得到
由這一個(gè)模型衍生出不同的模型想诅，這些模型具有不同水平的dN/dS異質(zhì)性。
However, the ratio averaged over all sites and all lineages is alomost never >1, since positive selection is unlikely to affect over all sites and all lineages over prolonged time. Thus, interest has been focused on detecting positive selection that affects only some lineages or some sites.
分支模型branch model
The branch models allow the?ω ratio to vary among branches in the phylogeny and are useful for detecting positive selection acting on particular lineages.
分支模型允許ω在系統(tǒng)發(fā)育樹(shù)上的不同分支上存在變異岛心，適用于檢驗(yàn)編碼基因在特定分支是否存在正選擇来破。
最簡(jiǎn)單的模型是one-ratio模型：假設(shè)所有的分支都有共同的dN/dS
另一個(gè)極端是Free-ratio模型：每個(gè)分支都有獨(dú)立的dN/dS
還有很多中間類型的模型，比如two-ratio模型：假設(shè)感興趣的分支具有的非同義替換和同義替換比率（ω1）和背景的非同義替換和同義替換比率（ω0）不同忘古。
The above models can be compared using thelikelihood ratio test to examine interesting hypothesis.
上述模型可通過(guò)似然比檢驗(yàn)進(jìn)行比較徘禁，驗(yàn)證感興趣的假說(shuō)。
比如one-ratio模型和two-ratio模型就可以一起比較髓堪，驗(yàn)證某一分支的dN/dS是否和其他分支不同送朱，這里two-ratio模型就是備擇假設(shè)(alernative hypothesis),one-ratio模型就是零假設(shè)（null hypothesis）
位點(diǎn)模型site models
下面不會(huì)講到位點(diǎn)模型的實(shí)際操作，這里簡(jiǎn)單介紹一下干旁。
The site models allow the ω ratio to vary among sites (among codons or amino acids in the protein).
位點(diǎn)模型允許不同位點(diǎn)（密碼子之間或蛋白的氨基酸之間）上的ω存在差異驶沼。
分支-位點(diǎn)模型brach-site models
The branch-site models allow ω to vary among sites in the protein and across branches on the tree and aim to detect positive selection affecting a few sites along particular lineages.
分支-位點(diǎn)模型允許編碼基因的不同位點(diǎn)在系統(tǒng)發(fā)育樹(shù)的不同分支的ω存在變異,目的是檢驗(yàn)正選擇是否作用于特定分支的一些位點(diǎn)。
Easycodeml里的分支位點(diǎn)模型用到的是下面這篇文章里提到的Test2：
Evaluation of an improved branch-site likelihoos method for detecting positive selection at the molecular level.
Test2, also known the branch-site test of positive selection, compares the modified model A with the corresponding null model with ω2=1 fixed (fix_omege=1 and omega=1)
Test2里用于似然比檢驗(yàn)的模型就是model A和對(duì)應(yīng)的null model争群。
先來(lái)介紹下model A：
modelA在似然比檢驗(yàn)中是備擇假設(shè)回怜。
在modelA里，系統(tǒng)發(fā)育樹(shù)被預(yù)先分為前景枝和背景枝换薄，只有前景枝才經(jīng)歷正選擇玉雾。模型假設(shè)有四個(gè)類型的位點(diǎn)（site class）：
class 0：包括在整個(gè)系統(tǒng)發(fā)育樹(shù)中保守的密碼子，dN/dS估算值:0<ω0<1
class 1:包括在整個(gè)系統(tǒng)發(fā)育樹(shù)中演化上中性的密碼子专控，dN/dS估算值:ω1=1
class 2a and 2b:包括在背景枝保守(2a抹凳，0<ω0<1)，在前景枝受正選擇（dN/dS估算值:ω2>1）的密碼子
class 2b:包括在背景枝演化上中性（2b伦腐，ω1=1）赢底，在前景枝受正選擇（dN/dS估算值:ω2>1）的密碼子
Test2中，和modelA進(jìn)行似然比檢驗(yàn)的零模型（null model）就是將modelA中ω2值固定，ω2=1幸冻。
在modelA里面粹庞，一些位點(diǎn)（2a和2b）的前景枝的ω值預(yù)設(shè)為ω2>1,因此Test2也是對(duì)前景枝某些位點(diǎn)是否經(jīng)歷正選擇的直接檢驗(yàn)。
進(jìn)行似然比檢驗(yàn)的方法是卡方分布洽损。
實(shí)操部分
實(shí)操部分介紹兩個(gè)方法做分支模型的檢驗(yàn)庞溜，
一個(gè)是windows系統(tǒng)下用可視化軟件EasyCodeML，
一個(gè)是在linux系統(tǒng)下用paml的codeml碑定。
1 EasyCodeML
更詳細(xì)的介紹可以參考高芳鑾老師的這篇文章：
https://blog.sciencenet.cn/blog-460481-1163040.html
EasyCodeML：A visualtool for analysis of selection using CodeML
Paml: a program package for phylogeneticanalysis by maximum likelihood
Paml 是使用最大似然法進(jìn)行系統(tǒng)發(fā)育分析的程序包流码，包括的程序有baseml、codeml延刘、evolver漫试、basemlg等等。
Paml：http://abacus.gene.ucl.ac.uk/software/paml.html
而EasyCodeML是使用codeml進(jìn)行選擇分析的可視化工具碘赖，也就是說(shuō)EasyCodeML可以實(shí)現(xiàn)Paml里面codeml程序的功能驾荣，而且提供可視化操作界面。
EasyCodeML：https://github.com/BioEasy/EasyCodeML/
下面是利用EasyCodeML中的分支模型（branch model）檢驗(yàn)感興趣的分支是否具有和其他分支不同的非同義替換和同義替換比率(dN/dS)的實(shí)踐普泡。
Easycodeml
1 ) 安裝軟件
軟件的下載地址：https://github.com/BioEasy/EasyCodeML/
下載的壓縮包解壓之后播掷，就可以直接雙擊打開(kāi)其中的應(yīng)用程序使用。
2 ) 準(zhǔn)備文件
需要準(zhǔn)備的文件包括序列文件和樹(shù)文件撼班，如果是在paml里運(yùn)算歧匈，還需要一個(gè)控制文件，但在EasyCodeML里不需要
序列文件就是不同分類群的要檢測(cè)的基因序列权烧，進(jìn)行translation align（這種比對(duì)方法是先將堿基翻譯成密碼子眯亦，然后將密碼子作為一個(gè)整體比對(duì)）。比對(duì)后般码，刪去終止密碼子妻率。如果有g(shù)ap，我采用的方法是直接刪去比對(duì)矩陣中有g(shù)ap的整列板祝。
上方工具欄Tools-Seqformat Convertor:將比對(duì)序列的.fasta或其他格式的文件轉(zhuǎn)化為.paml文件
直接將文件拖入original file框中宫静，右邊三角下來(lái)菜單中選擇轉(zhuǎn)換前的文件格式，生成的.paml文件會(huì)出現(xiàn)在轉(zhuǎn)換前的文件所在的文件夾中券时。
樹(shù)文件就是要研究的分類群的系統(tǒng)發(fā)育樹(shù)文件孤里，需要newick格式。分析時(shí)用到的樹(shù)的拓?fù)浣Y(jié)構(gòu)橘洞，而不用枝長(zhǎng)信息捌袜。要求樹(shù)的拓?fù)浣Y(jié)構(gòu)能反映真實(shí)的系統(tǒng)發(fā)育關(guān)系，因此這個(gè)系統(tǒng)發(fā)育樹(shù)不一定是基于要研究的基因建立的樹(shù)炸枣。
The tree topology, but not the branch length is used to fit different models.
需要注意的是序列文件和樹(shù)文件中的分類單元名稱要一致虏等，不要出現(xiàn)空格或分號(hào)弄唧。
3 ) 運(yùn)行EasyCodeML
runing model
preset是預(yù)設(shè)模式，適合不太熟悉分析的新手霍衫。
custom適合熟悉分析候引，希望自己設(shè)置參數(shù)的用戶。
Setup
Model Selection 根據(jù)目的選擇敦跌，本例這里選branch model
working directionary 指定生成文件所在路徑
If the option “Clean data” was selected, all sites with ambiguity characters and alignment? gaps will be removed from the aligned sequence file. Note that alignment gaps are treated as missing data.
clean data 如果選擇clean data選項(xiàng)澄干，對(duì)比序列中存在的gaps和ambiguity characters會(huì)被移除掉。
Aligned sequence file in palm format 指定paml格式的序列文件所在的路徑柠傍，
icode下拉選項(xiàng)中可以選擇編碼基因的序列類型麸俘，比如這里用的是昆蟲的線粒體基因，選擇4 invertebrate mt.
Tree File in Newick Format:指定newick格式的樹(shù)文件所在的路徑惧笛，點(diǎn)擊Check疾掰，檢查樹(shù)文件和序列文件的分類群名稱是否一致，點(diǎn)擊label對(duì)樹(shù)文件中感興趣的分支進(jìn)行標(biāo)記：
選中感興趣的分支徐紧，標(biāo)記成功后該分支會(huì)變黃，并帶有#1炭懊，也就是前景枝并级。這里只有一個(gè)前景枝，如果有兩個(gè)侮腹，選中第二個(gè)分支嘲碧，點(diǎn)擊2nd，以此類推父阻。其他的分支稱為背景枝
在自己的樹(shù)文件所在的文件夾下愈涩，會(huì)自動(dòng)生成一個(gè)在原樹(shù)文件名后帶后綴labeled的樹(shù)文件（如果有需要，可用于paml程序包的分析）
最后點(diǎn)擊Run CodeML
運(yùn)算過(guò)程可以在Runing status中查看加矛，如果遇到錯(cuò)誤也可以根據(jù)報(bào)錯(cuò)結(jié)果解決履婉。
4 ）結(jié)果解讀
程序運(yùn)行結(jié)束后，會(huì)彈出提示任務(wù)完成的小窗口斟览。
Summary of results處點(diǎn)擊export毁腿，給出結(jié)果(excel表格)的輸出路徑。結(jié)果呈現(xiàn)如下圖：
在整個(gè)運(yùn)算過(guò)程中苛茂，用到了兩個(gè)branch model已烤，
一個(gè)是model0，也就是之前提到過(guò)的one-ratio模型：假設(shè)所有的分支都有共同的dN/dS
一個(gè)是two ratio model 2:假設(shè)感興趣的分支和其他分支有不同的dN/dS
在two ratio model 2中妓羊，ω1代表的是前景枝的dN/dS的估算值胯究，ω0代表的是背景支的dN/dS的估算值
在Model 0中，ω代表的是在假設(shè)整個(gè)樹(shù)只有一個(gè)dN/dS的前提下躁绸，dN/dS的估算值
likelihood ratio test結(jié)果：
lnL代表的是模型的最大似然值的log值：
以上圖的結(jié)果為例裕循，Two-ratio model的lnL值為-5820.843753臣嚣，one-ratio model的lnL值為 -5822.141900，二者差值的二倍（1.20x2）與自由度為1（df=1）的卡方分布進(jìn)行比較费韭。
自由度的值為兩個(gè)模型參數(shù)個(gè)數(shù)（np, number of parameters）的差值.
p-value<0.05（顯著significant）或 p-value<0.01（極顯著extremely significant）茧球，說(shuō)明相比于one-ratio model，Two-ratio model更符合數(shù)據(jù)星持，選擇備擇假設(shè)（Two-ratio model）抢埋，
原始結(jié)果文件.mlc在指定的工作路徑下查看。
2 paml
以下操作都是在linux系統(tǒng)下進(jìn)行：
1）安裝paml
conda安裝：
conda install -c bioconda palm
2) 準(zhǔn)備文件
序列文件的準(zhǔn)備可以參考paml的案例文件督暂，這是我用的序列文件（我偷懶了揪垄，用Easycodeml里面的Seqformat Convertor轉(zhuǎn)化.fasta文件為.paml文件，再將后綴改為.nuc逻翁，順便說(shuō)下linux系統(tǒng)批量改文件后綴的命令：rename .paml .nuc? *.paml）饥努。
最上方的數(shù)字：59代表分類單元數(shù)，120代表比對(duì)序列長(zhǎng)度（核苷酸長(zhǎng)度）
樹(shù)文件的話需要準(zhǔn)備兩個(gè)八回，都為.newick文件酷愧，一個(gè)是沒(méi)有標(biāo)注感興趣分支的，一個(gè)是標(biāo)注感興趣分支的（我又偷懶了缠诅，標(biāo)注的樹(shù)文件也是通過(guò)easycodeml標(biāo)注的）
控制文件的內(nèi)容如下溶浴，本例中我們需要兩份控制文件，分別是one ratio model的控制文件和two ratio model的控制文件管引。
最上方的seqfile士败、treefile、outfile分別是告訴程序你的序列文件褥伴、樹(shù)文件和最后結(jié)果文件的絕對(duì)路徑
需要注意的幾個(gè)參數(shù)是：
runmode = 0谅将，意思是所用的樹(shù)文件是用戶提供的樹(shù)。（因?yàn)榛谒治龌驑?gòu)建的系統(tǒng)發(fā)育樹(shù)往往不能反映真實(shí)的系統(tǒng)發(fā)育關(guān)系）
icode = 4 代表分析的invertebrate mt（無(wú)脊椎動(dòng)物的線粒體基因）
fix_kappa = 0重慢， 代表kappa根據(jù)數(shù)據(jù)進(jìn)行估算（kappa to be estimated） ,kappa代表轉(zhuǎn)換率/顛換率（transition/transversion）
mode = 0饥臂， 代表的是該控制文件是運(yùn)算one ratio model的控制文件，即假設(shè)所有的分支都有共同的dN/dS
mode = 2似踱，代表該控制文件是運(yùn)算具有多個(gè)ratio model的控制文件擅笔，即假設(shè)分支上具有兩個(gè)及以上的dN/dS
3）運(yùn)行paml中的Codmel程序
我是在控制文件所在的路徑下運(yùn)行該程序。
#查看codeml的路徑：
which codeml
#運(yùn)行程序：one ratio model
~/miniconda3/envs/paml/bin/codeml ATP8_M0_codeml.ctl
在上面的命令中屯援，~/miniconda3/envs/paml/bin/codeml是codeml程序路徑烹俗，ATP8_M0_codeml.ctl是控制文件
#運(yùn)行程序：two ratio model
~/miniconda3/envs/paml/bin/codeml ATP8_BM_codeml.ctl
運(yùn)算結(jié)束后衡蚂，結(jié)果會(huì)出現(xiàn)在控制文件指定的輸出路徑下
4）結(jié)果解讀
在結(jié)果文件中找到以下數(shù)據(jù)(這里以one ratio model 運(yùn)算結(jié)果為例)：
模型的lnL值，omega值（ω）。自己可以進(jìn)行最大似然比檢驗(yàn)循衰，或者用Easycodeml上方工具欄中Tools下面的最大似然比檢驗(yàn)工具LRT's calculator
其解讀與上述Easycodeml部分相同轧邪。
