最近被數(shù)據(jù)分析虐的死去活來，好難。
參考文獻(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 非同義替換
同義替換指：由于密碼子具有簡并性，密碼子中一些核苷酸的替換并不會導(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á)出來（編碼蛋白發(fā)生改變），自然選擇才能起作用。
如果一個(gè)分支的基因變異受到自然選擇的正選擇（positive selection）（正選擇指的是自然選擇支持變異，變異是有利的），那么該分支的dN/dS應(yīng)該與其他分支不同，或者大于1。
需要注意的是某基因在某一分支的dN/dS大于1來證明該基因在該分支上受正選擇是很保守的，dN/dS很少會大于1，如果某一個(gè)分支的dN/dS值顯著大于其他分支，也能說明問題。
如果分支的某基因選擇上中性（neutral selection），即無害也無利，那么理論上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ā)育樹上的不同分支上存在變異，適用于檢驗(yà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.
上述模型可通過似然比檢驗(yàn)進(jìn)行比較，驗(yàn)證感興趣的假說。
比如one-ratio模型和two-ratio模型就可以一起比較，驗(yàn)證某一分支的dN/dS是否和其他分支不同，這里two-ratio模型就是備擇假設(shè)(alernative hypothesis),one-ratio模型就是零假設(shè)（null hypothesis）
位點(diǎn)模型site models
下面不會講到位點(diǎn)模型的實(shí)際操作，這里簡單介紹一下。
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ā)育樹的不同分支的ω存在變異,目的是檢驗(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和對應(yīng)的null model。
先來介紹下model A：
modelA在似然比檢驗(yàn)中是備擇假設(shè)。
在modelA里，系統(tǒng)發(fā)育樹被預(yù)先分為前景枝和背景枝，只有前景枝才經(jīng)歷正選擇。模型假設(shè)有四個(gè)類型的位點(diǎn)（site class）：
class 0：包括在整個(gè)系統(tǒng)發(fā)育樹中保守的密碼子，dN/dS估算值:0<ω0<1
class 1:包括在整個(gè)系統(tǒng)發(fā)育樹中演化上中性的密碼子，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也是對前景枝某些位點(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)行選擇分析的可視化工具，也就是說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/
下載的壓縮包解壓之后，就可以直接雙擊打開其中的應(yīng)用程序使用。
2 ) 準(zhǔn)備文件
需要準(zhǔn)備的文件包括序列文件和樹文件，如果是在paml里運(yùn)算，還需要一個(gè)控制文件，但在EasyCodeML里不需要
序列文件就是不同分類群的要檢測的基因序列，進(jìn)行translation align（這種比對方法是先將堿基翻譯成密碼子，然后將密碼子作為一個(gè)整體比對）。比對后，刪去終止密碼子。如果有g(shù)ap，我采用的方法是直接刪去比對矩陣中有g(shù)ap的整列。
上方工具欄Tools-Seqformat Convertor:將比對序列的.fasta或其他格式的文件轉(zhuǎn)化為.paml文件
直接將文件拖入original file框中，右邊三角下來菜單中選擇轉(zhuǎn)換前的文件格式，生成的.paml文件會出現(xiàn)在轉(zhuǎn)換前的文件所在的文件夾中。
樹文件就是要研究的分類群的系統(tǒng)發(fā)育樹文件，需要newick格式。分析時(shí)用到的樹的拓?fù)浣Y(jié)構(gòu)，而不用枝長信息。要求樹的拓?fù)浣Y(jié)構(gòu)能反映真實(shí)的系統(tǒng)發(fā)育關(guān)系，因此這個(gè)系統(tǒng)發(fā)育樹不一定是基于要研究的基因建立的樹。
The tree topology, but not the branch length is used to fit different models.
需要注意的是序列文件和樹文件中的分類單元名稱要一致，不要出現(xiàn)空格或分號。
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)，對比序列中存在的gaps和ambiguity characters會被移除掉。
Aligned sequence file in palm format 指定paml格式的序列文件所在的路徑，
icode下拉選項(xiàng)中可以選擇編碼基因的序列類型，比如這里用的是昆蟲的線粒體基因，選擇4 invertebrate mt.
Tree File in Newick Format:指定newick格式的樹文件所在的路徑，點(diǎn)擊Check，檢查樹文件和序列文件的分類群名稱是否一致，點(diǎn)擊label對樹文件中感興趣的分支進(jìn)行標(biāo)記：
選中感興趣的分支，標(biāo)記成功后該分支會變黃，并帶有#1，也就是前景枝。這里只有一個(gè)前景枝，如果有兩個(gè)，選中第二個(gè)分支，點(diǎn)擊2nd，以此類推。其他的分支稱為背景枝
在自己的樹文件所在的文件夾下，會自動(dòng)生成一個(gè)在原樹文件名后帶后綴labeled的樹文件（如果有需要，可用于paml程序包的分析）
最后點(diǎn)擊Run CodeML
運(yùn)算過程可以在Runing status中查看，如果遇到錯(cuò)誤也可以根據(jù)報(bào)錯(cuò)結(jié)果解決。
4 ）結(jié)果解讀
程序運(yùn)行結(jié)束后，會彈出提示任務(wù)完成的小窗口。
Summary of results處點(diǎn)擊export，給出結(jié)果(excel表格)的輸出路徑。結(jié)果呈現(xiàn)如下圖：
在整個(gè)運(yùn)算過程中，用到了兩個(gè)branch model，
一個(gè)是model0，也就是之前提到過的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è)樹只有一個(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），說明相比于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，順便說下linux系統(tǒng)批量改文件后綴的命令：rename .paml .nuc? *.paml）。
最上方的數(shù)字：59代表分類單元數(shù)，120代表比對序列長度（核苷酸長度）
樹文件的話需要準(zhǔn)備兩個(gè)，都為.newick文件，一個(gè)是沒有標(biāo)注感興趣分支的，一個(gè)是標(biāo)注感興趣分支的（我又偷懶了，標(biāo)注的樹文件也是通過easycodeml標(biāo)注的）
控制文件的內(nèi)容如下，本例中我們需要兩份控制文件，分別是one ratio model的控制文件和two ratio model的控制文件。
最上方的seqfile、treefile、outfile分別是告訴程序你的序列文件、樹文件和最后結(jié)果文件的絕對路徑
需要注意的幾個(gè)參數(shù)是：
runmode = 0，意思是所用的樹文件是用戶提供的樹。（因?yàn)榛谒治龌驑?gòu)建的系統(tǒng)發(fā)育樹往往不能反映真實(shí)的系統(tǒng)發(fā)育關(guān)系）
icode = 4 代表分析的invertebrate mt（無脊椎動(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é)果會出現(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部分相同。
