防杠疊甲：下述所有內(nèi)容均為本人理解，歡迎友好的學(xué)術(shù)討論。覺得看不進去的，寫的不好的請windows用戶移步右上角×標(biāo)記，mac用戶移步左上角，手機用戶請關(guān)機。
持續(xù)更新中
將依據(jù)[Beginner’s Guide on the Use of PAML to Detect Positive Selection]介紹PAML的正選擇有關(guān)內(nèi)容
一 不同的密碼子替換模型檢測適應(yīng)性進化
- branch model
branch model 假設(shè)系統(tǒng)發(fā)育樹上的不同枝有不同的 ω 值。它們可以用來檢測作用于特定譜系的正選擇，而不需要對整個系統(tǒng)發(fā)生樹的 ω 比率進行平均。例如，在檢測基因加倍時間后的正選擇是useful的，在這種情況下，其中一份拷貝可能獲得了新功能，因此有一個更快的進化速率。 -
site model
site model 將基因中任何位點(密碼子)的 ω 比率視為來自統(tǒng)計分布的隨機變量，因此允許 ω 在密碼子之間變化。ω > 1的密碼子被定義為正選擇。
Likelihood ratio test(LRT)被用來比較零假設(shè)（ω < 1）和備擇假設(shè)（ω > 1）。
可以使用Na?ve Empirical Bayes (NEB) method 和Bayes empirical Bayes (BEB)檢測哪個密碼子受到正選擇。
圖1 - Branch-site models
Branch-site models 旨在檢測預(yù)先設(shè)置的譜系上的一少部分（only a few sites）的正選擇。受到正選擇檢測的值叫前景枝（foreground branches，ω2 > 1）, 其他的枝叫背景枝（background branches）。背景枝有兩類位點，保守的位點（0 < ω0 < 1）和中性的位點（ω1 = 1）?？ǚ綑z測鑒定模型選擇，BEB檢測受選擇位點。
下面將會介紹CODEML的四種分析模式：
1）所有位點、所有枝系具有相同的ω
2）檢測正選擇進化下一部分位點的正選擇（site model）
3）檢測系統(tǒng)發(fā)育樹種一個特定的枝或一些特定的枝的正選擇（branch model）
4）特定枝一些位點的正選擇（branch site model）
Protocol
The Control File
control file包含了COMEML運行需要的參數(shù)，所有的可以參考PAML的說明，這里只提及分析有關(guān)的參數(shù)。
control file的一些規(guī)則：
1）一行一行讀，如果有相同的兩行，第二行的內(nèi)容會覆蓋第一行。
2）一些內(nèi)容對不同的model都適用。
3）空行會被程序忽略，*之后的內(nèi)容會被認(rèn)為是注釋。
The Input/Output Files
control file的第一塊指定了輸入文件的路徑：
1)序列比對文件(seqfile)，PHYLIP格式，每一個基因的比對文件都有一個開頭(header, 例如"12 1989")，如果一個文件有多個基因比對，需要多個header。
2)系統(tǒng)發(fā)育樹(treefile)，Newick格式，需要刪除bootstrap值，枝長信息最好刪除?？赡苄枰粋€header（自身體會好像不是剛需），如果有多個樹，header是必要的。
3)主要結(jié)果儲存的輸出文件(outfile)
Screen Output
輸出到屏幕和輸出文件的內(nèi)容總量是由noisy和verbose控制。
Input Sequence Data
用來推斷的數(shù)據(jù)類型（e.g., nucleo-tides, amino acids, or codons that are to be translated into amino acids by CODEML）是由seqtype控制。
ndata: 基因或者alignments的數(shù)量
icode: genetic code
cleandata: if sites with ambiguity data and alignment gaps should be kept (cleandata = 0) or removed (cleandata = 1)
在使用CODEML之前，序列是要對齊的，內(nèi)含子，非編碼區(qū)以及終止子要去除。為了保持閱讀框架（reading frame），一個有用的策略是首先對齊蛋白質(zhì)序列，然后相應(yīng)地構(gòu)建密碼子對齊序列。推薦刪除主要是gap的列或者區(qū)域，然后使用cleandata=0
Substitution Model
幾個參數(shù)用于指定分析中使用的模型。第一，ω隨著譜系變化（model）,隨著位點變化（NSsites）,或者兩者都有（model and NSsites）。第二，密碼子頻率（CodonFreq）是可以相同的（1/61）或者是不同的來定義密碼子bias。FmutSel模型（CodonFreq=7）定義每個密碼子的fitness都是60（=61-1）。FmutSel0模型（CodonFreq=6）是FmutSel的特殊情況，對于同義密碼子的fitness是一樣的，因此只使用19（=20-1）個氨基酸fitenss參數(shù)。該模型假設(shè)氨基酸頻率是由蛋白質(zhì)的功能需求決定的，但同義密碼子的相對頻率完全由突變偏好參數(shù)決定。在這些突變選擇模型下，estFreq使用從數(shù)據(jù)中觀察到的密碼子頻率或最大似然估計密碼子頻率。取決于用這些參數(shù)定義的模型，ω值可以通過fix_omega和omega來指定的或由軟件計算。clock指定分子鐘(i.e., rate constancy among lineages）。因為我們在這里討論的模型都是時間可逆的，因此我們不假設(shè)分子鐘（clock=0），并且系統(tǒng)發(fā)育樹是無根樹（unrooted）。Note，幾乎所有的系統(tǒng)發(fā)育樹構(gòu)建軟件例如RAxML, IQ-TREE和MrBayes都是無根樹。
One Ratio With Homogeneous ω Across Lineages and Sites
CODEML中最簡單的模型就是M0（one-ratio）模型,其假定多有譜系的所有位點只有一個ω值（圖3A）。大多數(shù)情況下，這種假設(shè)是不現(xiàn)實的，因為大多數(shù)位點都被期望在ω < 1的約束下。因此，嘗試在M0的假設(shè)下尋找正選擇的證據(jù)是會失敗的。事實上，當(dāng)我們想鑒定正選擇時，最好的方式是在后續(xù)的描述中使用site或branch-site模型計算ω。然而，M0模型給出了一個零假設(shè)，可以作為與其他模型比較的參考，以確定更復(fù)雜的模型是否能更好地擬合數(shù)據(jù)。M0模型下計算出的ω也給出了在位點和譜系上平均的選擇性約束的總體度量。以 codeml codeml-M0.ctl運行后，一些內(nèi)容會打印在屏幕上，一些內(nèi)容會保存到輸出文件（out_M0.txt）內(nèi)，輸出文件可以分為5個部分。
1. 輸入比對文件的位點模式總結(jié)Summary of the Site Patterns in the Input Alignment
輸入的比對序列和相對應(yīng)的簡略比對在輸出文件的頂部被輸出，每一個位點模式只被呈現(xiàn)一次，相對應(yīng)的頻率（pattren counts）在下方的空白顯示。如果使用cleandata=1，刪除gap和ambiguous位點前后的比對情況都會被打印。注意，如果gap被刪除后，位點會重新編號，這會影響site和branch-site模型下的輸出。
2. 輸入比對和選擇模型總結(jié)Summary of the Input Alignment and the Model Selected
CODEML的版本和比對文件的名字會一并輸出。后面會是control file里指定模型的具體信息和物種個數(shù)以及密碼子個數(shù)。
3. 核酸和密碼子頻率的總結(jié)Summary of Nucleotide and Codon Frequencies
每個序列中觀察到的核苷酸和密碼子頻率以及它們在所有序列中的平均值被打印在輸出文件中，然后是模型下的密碼子頻率，可以使用諸如evolver之類的軟件進行模擬。
4. 樹得分和模型參數(shù)值計算總結(jié)Summary of Tree Scores and Estimated Model Parameter Values
某一模型下的log-likelihood，free parameters的數(shù)量，樹的枝長和其他的參數(shù)在前四行輸出。 Free parameters包括枝長，equilibrium frequencies，轉(zhuǎn)換和顛換的比率（κ）和omega分布參數(shù)。有n個物種的無根樹有2×n-3的枝長，3個mutation-bias parameters （four frequencies with the sum to be 1）。
因為指定的是M0模型(model = 0 and NSsites = 0)，每一個枝有相同的ω 值(0.3357) 。The tree length is dN = 1.1386 at the nonsynonymous sites and dS = 3.3921 at the synonymous sites, with ω = dN/dS = 0.3357，意味著該基因總體受到凈化選擇，即非同義突變的固定概率(即ω = dN/dS = 0.3357)平均為同義突變的三分之一。
在下一個模塊中，我們聚焦三個鑒定正選擇的模型，site，branch，branch-site模型
位點ω不一致的Site模型
在這個板塊中，在允許ω沿密碼子變化的site模型下運行CODEML（圖3B）。在這個例子中，運行CODEML基于M0模型（類似之前）和下面的位點不一致模型：M1a,M2a,M7和M8（table1）。編輯配置文件，并重命名為codeml-sites.ctl，并運行CODEML。
輸出文件包含控制文件中定義的五個模型中的每個模型的一個部分，其中提供了log-likelihood 和總參數(shù)的數(shù)量，后者用于確定每個模型比較的自由度。
為了比較nested位點模型，使用LRT來進行檢驗，其定義為零假設(shè)和備擇假設(shè)的log-likelihood差值的2倍，2Δ? = 2(?1 ? ?0), where ?0 is the log-likelihood score for the null model, whereas ?1 is the log-likelihood under the alternative model。LRT統(tǒng)計量與χ2分布進行比較，其自由度等于兩個模型之間自由參數(shù)(在每個模型的輸出文件中給出)數(shù)量的差異。結(jié)果總結(jié)在表3中。
M0（one-ratio）和M1a（nearly netural）模型可以使用LRT相比。這是對氨基酸位點之間的選擇壓力變異性的測試，而不是對正選擇的測試。相比于M0，M1a更適合示例數(shù)據(jù)，with 2Δ? = 559.26，表明ω所反映的選擇壓力在不同位點之間差異很大。
與M1a相比，M2a增加了一些受到正選擇（ ω2 > 1）的位點。這并沒有提高模型的適合度因為2Δ? = 0（table3）。
同時，做了另一個正選擇測試通過比較M7(beta, null model)和M8 (beta&ω, alternative model)。M8更適合數(shù)據(jù)，with 2Δ? = 12.54，意味著存在帶ω >1的正選擇位點。陽性選擇位點的測試與M1a-M2a和M7-M8比較給出了相互矛盾的結(jié)果。當(dāng)正選擇的證據(jù)存在但不是很強時，注意到M1a-M2a檢驗更為嚴(yán)格，因為弱正選擇的位點往往被集中到ω1 = 1的位點類別中。
site模型下的參數(shù)的最大似然估計（MLE）列在table4?？梢允褂胋ash腳本從輸出文件中提取這些值?；蛘?，我們在GitHub存儲庫中的分步教程中包含代碼片段。
https://github.com/abacus-gene/paml-tutorial/tree/main/positive-selection/00_ data
盡管M0假定在一個基因中所有的密碼子有相同的ω值（model = 0 and NSsites = 0），位點模型假定了幾種位點類型（圖3A、B）。例如，在M8模型下，94.2%的位點（p0）具有ω的beta分布，然而5.8%的位點的ω = 1.841，表明存在一小部分氨基酸殘基在正選擇下。這些信息能在輸出文件中找到，其位置在MLE開頭的下一行。再一次運行中如果指定了多個位點模型，每一個模型在輸出文件中有自己的輸出信息塊，在這里也會顯示估計的參數(shù)。
在M2a和M8模型下，使用BEB（Bayes Empirical Bayes）方法計算來自不同位點類別的每個位點的后驗概率。在正選擇類別下具有高后驗概率的位點更可能受到正選擇。M8模型下的輸出如下所示：
在每一行，第一列顯示位點的位置（例如，10,25,108，123），后面是氨基酸的類型。第三列（Pr (w > 1)）顯示該位點來自正選擇類的后驗概率。最后一列顯示了該位點ω分布的后驗均值和標(biāo)準(zhǔn)差。注意，BEB計算僅在正選擇模型(即M2a和M8)下進行，而不是在零模型(即M1a或M7)下進行。
在該例子中，14個位點具有大于50%的概率有正選擇（ ω > 1），上述只列出4個。例如，第一個序列的位點10是絲氨酸，是正選擇類的概率是0.54，ω后驗概率的分布的平均值是1.468，標(biāo)準(zhǔn)差是0.634。在這里，BEB的計算只是提供改位點受到正選擇的微弱證據(jù)，因為ω > 1的后驗概率很低并且ω的分布在每個位點是彌漫分布。同時，M1a-M2a的LRT比較不顯著并且M7-M8的比較只在5%的水平下顯著?？偟膩碚f，這些結(jié)果顯示在這個基因上一些位點受到了正選擇，但強度不夠強。
沿分枝的不同ω的Branch模型
在這個模塊，將展示使用CODEML運行branch模型（圖3C），其假定ω沿樹的不同分枝變化。分支模型是通過在樹文件中使用標(biāo)記（tags）標(biāo)記分支來指定的：#0（默認(rèn)），#1等。至多允許8個不同分枝的ω。在這里我們只考慮兩種類型。前景枝（foreground）用#1標(biāo)記（ ω1），其他的枝作為背景枝且標(biāo)記為#0（默認(rèn)，ω0）。注意，#0不用在樹文件中標(biāo)記?？梢栽贜ewick樹中手動添加前景枝標(biāo)記，或者使用GitHub中的腳本(https://github.com/abacus-gene/paml-tutorial/tree/main/positive-selection/00_data)，或者使用PhyloTree和EasyCodelML。
這里，我們介紹四種測驗，設(shè)置了不同的枝作為前景枝：（1）chicken，（2）duck，（3）chicken和duck，（4）chicken和duck以及他們的共同祖先。零假設(shè)在所有的測驗里都是M0，其假設(shè)所有枝ω相同。除了測驗4，其他測驗均使用無根樹。
在第一個測驗里，我們想知道chicken是否受到正選擇。樹文件如下所示：
然后，在其他三種不同的檢驗下運行相同的分析。配置文件和之前相同，除了以下的部分：
通過使用model=2，我們指定與背景枝相比，前景枝可能在不同的ω下進化。運行其他檢驗時，注意更改相應(yīng)的輸出文件名。
輸出文件格式與之前相同，唯一的區(qū)別是輸出文件中有一個額外的塊，其中包含輸入樹的每個分支的dN/dS比率。
以chicken為例，輸出結(jié)果如下：
在此樹中，每個分支的估計ω比顯示在“#”之后。這里，所有背景枝的ω0 = 0.293，而前景枝的ω1 = 1.171。
注意到以duck為前景枝時，ω1 = 999。999是數(shù)值的上限意味著無窮，如果在有關(guān)分支上缺乏同義替換，則可能發(fā)生這種極端估計。在這種情況下，雖然不能準(zhǔn)確的估計ω1，但LRT依然有效。檢驗3時，ω0 = 0.287，ω1 = 0.711，與檢驗4一致。結(jié)果顯示，與其他枝相比chicken和duck有更高的ω值，暗示了可能的正選擇。
可以使用LRT檢驗假設(shè)的顯著性。根據(jù)LRT，我們得出結(jié)論，分支模型比M0模型更適合所有測試假設(shè)的數(shù)據(jù)。換言之，不同檢驗得到的ω值與背景枝的ω比有顯著不同。
在這里，我們展示了(1)如何在分枝模型下執(zhí)行CODEML，對樹上不同的進化譜系或分枝假設(shè)不同的選擇壓力(即不同的ω比)，以及(2)如何使用M0作為零模型進行影響預(yù)先指定分支的正選擇的LRT。我們進行了四個測試，部分是為了說明，如果模型為根周圍的兩個分支指定了不同的進化過程，則可以使用有根樹。注意，測試中前景分支的識別取決于被測試的生物學(xué)假設(shè)，這應(yīng)該是先驗的。如果在沒有任何先驗假設(shè)的情況下，將樹的每個分支依次指定為前景進行分支測試，則需要對多重測試進行校正。
沿分枝和位點的不同ω的Branch-site模型
在本節(jié)中，我們在Branch-site模型下運行CODEML，其假設(shè)ω在譜系和位點之間都是變化的（圖3D）。該模型可用于檢測沿預(yù)先指定的前景枝影響特定氨基酸位點的正選擇。
樹文件和branch模型的格式相同，前景枝用#1做標(biāo)簽。我們使用之前branch模型的樹文件和按如下修改branch模型的配置文件：
輸出結(jié)果包括了一個模塊，其中列出了branch-site模型A中假設(shè)的k = 4個位點類的比例估計值以及背景和前景分支的ω值（table 6）。site class 0 從0-1（0 < ω0 < 1），對于前景枝和背景枝，這一類的位點經(jīng)歷純化選擇。在site class 1，ω值固定為1。這些位點受到中性選擇。在site class 2a 或 2b，前景枝經(jīng)歷正選擇（ω2 ≥ 1），背景枝受到純化選擇（ 0 < ω0 < 1，site class 2a）或者經(jīng)歷中性進化（ ω1 = 1，site class 2b）。
計算出的ω顯示（table 6）在雞和鴨譜系中存在正選擇的位點。事實上，正選擇似乎影響了鳥類進化中的所有三個分支。
