『壹』 paml計算 KaKs值
PAML 可實現系統發育樹的構建，祖先序列估計,進化模擬和 KaKs 計算等功能。其中分支及 位點 KaKs 的計算是本軟體包的特色功能世或。
此次用到的是codeml
所需文件：
本次使用我最近的數據，來源於2個物種的同源基因對，具體如何得到基因對請挪步 python版的MCScan繪圖 。
使用
上述腳本有疑問請挪步 Kaks_calculator計算ka/ks 值
上述得到一paml_result文件夾，每個同源基因對兒形成一個單獨的以*.paml結尾的文件
可獲得共有27個同源基因對
將所有*paml文件合並為paml的輸入文件
關於樹文件，可參考paml安裝目錄下*.trees格式
其中3表示，3個物種，4表示樹的個數；
在本次我只有兩個物種，所以得到如下樹的輸入文件
可將paml安裝目顫汪錄下baseml.ctl 拷貝到自己所需目搜洞伍錄下即可進行修改
即可得到相應的Ka，Ks
『貳』 PAML-discussion-group
PAML 是一個用最大似然法來對DNA和蛋白質序列進行系統發育分析的軟體包，此軟體由楊子恆教授開發，最新版本為v. 4.9h。
此筆記主要用於記錄PAML分析中遇到的一些問題，並長期更新。
PAML主頁
PAML手冊
PAML FAQs
PAML discussion group
PAML理論基礎原則 蛋白編碼基因（protein coding sequence）的自然選擇壓力水平可以通過dN/dS(ω)值的大小來衡量，其中，dS代表同義替換率（synonymous rate），dN代表非同義替換率（non-synonymous rate）。在沒有受到選擇壓力時，同義替換率和非同義替換率相等，此時dN/dS = 1；當受到負選擇或凈化選擇壓力時，自然選擇會阻止氨基酸發生改變，同義替換率會大於非同義替換率，即dN/dS < 1;當受到正選擇壓力時，氨基酸的置換率蔽茄會受自然選擇的青睞，即蛋白功能可能會發生適應性改變，此時dN/dS > 1。
PAML簡易流程
要正常運行paml軟體，需要四個標准文件：
樹文件設置:
需要注意的是，Newick格式的樹尾部一定要有分號，沒有的話程序無法正常運行。
在paml分析裡面，主要涉及四種模型：位點模型（site-model）、支模型（branch model）、支位點模型（branch site model）以及進化支模型。位點模型通常適用於檢測某一支系普遍性、廣泛性的正選擇，這種正選擇是由於位點持續改變所引起的，例如，適應多種病原體；支模型主要是檢測某一支系是否存在快速進化、選擇壓力約束以及正選擇，但卻無法檢測到正選擇位點；支位點模型比較准確與穩定，適用於檢測某一支系斷點性的正選擇事件，此結果是由於適應某一時期環境改變所引起的，通常會保留在後代中；進化支模型（clade model），則主要是判別不同物種之間是否存受分化選擇壓力作用，不局限於正選擇，它一次可以標記多個分支進行比較。
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且呈現beta分布；
M8 ：在M7模型的基礎上，增加一類正選擇位點（ω >1 ）;
M8a :與M8類似，只不過是將新增的ω固定為1；
位點模型的Codeml.ctl參數設置：
one rario ：假設所有的進化譜系都具有相同的ω值；
free ratio ：假設所有的支系都具有獨立的ω值；
two ratio ：假設前景支與背景支的ω不同
分支模型的Codeml.ctl參數設置：
假定位點間的ω值是變化的，同時也假定支系間的ω值是變化的。該模型主要用於檢測前景支中正選擇作用對部分位點的影響。
modelA null（零假設） ：ω值設定則並侍為固定值1
modelA（備擇假設） ：估算其ω值是否大於1
背景支與前景支具有相同的位點ω值：
K0：前景支與背景支中的位點受到純化選擇0 < ω < 1；
K1：前景支與背景支中的位點處於中性進化0 < ω = 1;
背景支與前景支具有不同的位點ω值：
K2a：前景支處於中性進化，而背景支處於純化選擇；
K2b：前景支受到正選擇壓力（ ω>1 ），而背景支處於中性進化；
支位點模型的的孫吵Codeml.ctl參數設置：
與枝位點模型類型，能同時檢測多個進化枝（Clade），但是該模型並沒有將背景支的dN/dS值約束在（0,1）。
簡單的區別方法就是，有根樹的祖先節點是二叉樹，而無根樹為三叉樹，舉個例子：
在使用codeml時，如果沒有指定有根樹參數卻使用了有根樹作為輸入，那麼在輸出結果中將會得到這樣的報錯信息： "This is a rooted tree. Please check!" 。對於大多數模型，即使使用有根樹，其該模型似然值仍然是正確的，但是root周圍兩個分支的長度不穩定，因為它們的和是估計值。對於其他模型，似然估計和參數估計都是不正確的。 因此，分析時確實應該注意到這一信息，並盡可能使用一棵無根樹。 我們可以使用R包 ape 將有根樹轉換為無根樹：
在多序列比對過程中，對齊gap是極其困難的，paml軟體包沒有辦法無法處理gap。因此，我們可以通過設置 cleandata = 1 來去除gap；此外，還可以將gap當作為模糊字元進行處理。但是，這都不是最好的解決辦法，這兩種策略都低估了序列差異。 就個人而言，我認為除了一個或兩個序列之外，大多數序列都有序列信息的位點也許應該保留，而除了一個或兩個序列之外，所有序列都有對齊間隙的位點最好被移除。因此，選擇合適的多序列比對軟體以及過濾軟體，尤其重要。
如果檢測前景支時，其dN/dS > 1時，我們可以認為它受到正選擇作用。但是，如果其dN/dS <1但大於背景支時，就不能認為它是受正選擇作用的，選擇壓力約束放鬆可能是較為合理的解釋。此外，在 Ohta's 的微有害突變假說下，凈化選擇在大種群中比在小種群中更有效，因此不同譜系的種群規模的差異提供了另一個相容的假設。如果氨基酸的變化是稍微有害的，我們預計在大群體中它們從群體中移除的速度會比在小群體中更高。因此，即使兩系在選擇壓力或基因功能上沒有差異，我們也期望在一個大群體中看到一個較小的dN/dS比值，例如，許多核基因的dN/dS比值在嚙齒類動物中低於靈長類或偶蹄類。
通常，如果一個位點在一種模型下出現在列表中，那麼在另一種模型下也會有相當大的概率。如果你以這種方式看待它們，結果可能沒有太大不同。確定位點的問題很困難，而且容易出錯。這種情況類似於從一個班級中得到排名前幾的學生。列表包含的內容越多，質量就越差。因此，我們通常會認為後驗概率大於95%或者99%的位點，是比較可信的。
如果你是關注某一個節點，則可以使用"#"；如果是關注某一類群，則可以使用"$"。對於位點模型以及free ratio分析，無需標注前景支；對於支模型以及進化支模型，單次分析可以標注多個前景支；但支位點模型單次分析僅能標注一個前景支。
上面的兩個例子其實是等同的，因此 $1 可以標記大的分支，包括其祖先節點以及現生節點。而 #1 僅僅是表示末端枝或者祖先節點。
這里有一些關於嵌套進化枝標簽的規律。符號#的優先順序高於$，樹頂端的進化枝標簽優先於接近根處的祖先節點的進化枝標簽。下面的兩棵樹也是一樣。第一個樹中$1適用於整個胎盤哺乳動物的進化枝（除人類世系外），$2適用於兔與鼠的進化枝。
使用TreeView軟體可以很方便的創建樹文件並且檢查樹和標簽是否正確。上面作為例子的樹都可以被TreeView識別。在TreeView X中，你必須使用單引號來標記。如下：
此外，還可以同時標注多個支為統一前景支，例如
鏈接1
鏈接2
[圖片上傳失敗...(image-4ab1a5-1584449462345)]
某一基因，在動物祖先發生基因復制，分化為A，B兩個不同支系。
較高的dN/dS值可以解釋為正選擇或者快速進化。盡管突變本身也處於選擇壓力之中（大多數為凈化選擇），但不可以解釋為「基因A增加了突變的選擇壓力」，因為突變是隨機發生的。原則上講，突變率會同時影響dN、dS，但通常dN/dS不受突變率影響。
dN/dS是一種進化速率，但它不是突變速率，因為同義和非同義替代率具有不同的選擇約束水平。 選擇壓力測試的基本原理是：它假設同義替換是中性進行，也就是說，它們大多是在遺傳漂移下進化的。如果這是真實情況，那麼dS可以用作（中性）突變率的替代。然而，非同義替代率總是處於凈化選擇壓力，在正選擇下程度較小。因此，dN/dS是中性偏離的度量。所以，dN > dS，既dN/dS > 1,則受到正選擇；如果dN小於dS，則dN/dS < 1，則是凈化選擇。選擇壓力測試的關鍵就是它通過特定基因的同義替換的「中性」進化速率歸一化非同義替換的速率。
在任何情況下，使用代表真實演化歷史的基因樹都是最好的。但是，有時可能無法輕易判斷是否符合真實演化歷史，那麼你可以選擇物種樹作為代替。基因組水平分析，那麼推薦使用物種樹。你可以可以使用基因樹和物種樹進行數據穩健性測試。
如果以哺乳動物祖先為前景，假設該基因在共同祖先中存在適應性，可能是由於獲得了新的功能，但隨後該基因在凈化選擇下得到了保守進化。如果你把整個clade作為前景，那麼假設該基因在整個哺乳動物中的所有分支都承受著持續的壓力來改變或多樣化，如果該基因涉及到防禦或免疫，則可能是這種情況。
對於到底是檢測祖先，還是整個枝系，取決於生物學問題。例如，溶菌酶在所有colobine猴子中都應該具有相同的功能，因此蛋白質在clade內被期望受到選擇性約束，但在colobine clade的分支上，該酶顯然獲得了一個新的功能，其正選擇驅動了氨基酸的變化。有了這個假設，你就應該把分支的祖先貼上clade的標簽，而不是那些在clade中的分支。
經常遇到審稿人的審稿意見是這樣的：前景分支上正選擇的基因的重要支持並不意味著在背景分支上沒有正選擇，這些基因可能仍在許多（即使不是全部）背景分支上處於正選擇狀態。為了檢驗原假設（基因僅在前景支受到正選擇）是否正確，可以進一步通過Clade模型檢驗，因為進化枝模型允許在不限制背景dN/dS小於1的約束下估計前景支與背景支dS/dN的比率。
一種可能性是，對整個基因有積極選擇的證據，但每個單獨位點的信息或證據太弱。 可以查看rst文件，該文件具有所有站點的後驗概率，以查看是否是這種情況，mlc文件只列出後驗概率高於0.5的文件。
可能是因為codeml在刪除gap或模糊字元的列，之後重新編號位點了（cleandata =1）。
遇到這種極大值的dN/dS時，比如ω = 999，首先要確保你的序列是否正確；其次，該位置的dn，ds是否遠小於0.0001，枝長是否過小。顯然，高度相似的序列和非常發散的序列都不是信息豐富的，很難指定確切的值。為了避免此類問題的出現，可以先通過M0模型獲得枝長，再將具有分支長度的進化樹應用到codeml，並在ctl中設置FIX_blength=2。
如圖所示，紅色分支代表表型趨同的進化分支，如果你想通過分支位點模型來檢測適應性趨同進化，應該將所有紅色分支統一設置為前景分支。當然，前提是假設所有前景支均具有相同的位點受到正選擇。 對於背景支是否同樣存在類似的適應性趨同，可以通過clade進化支模型進行檢測。
p0/ω0,p1/ω1,p2=(1-p0-p1)/ω2 ：在正選擇分析的備擇假設結果文件中，通常會獲得三個p值，其中，p0代表處於凈化選擇下的位點概率；p1代表處於中性進化下的位點概率；p2則代表處於正選擇下的位點概率。
可以使用two ratio備擇假設（fix_omega = 0 omega = 1）與two ratio零假設（fix_omega = 1 omega = 1）進行統計假設檢驗。
codeml檢測選擇約束放鬆可以通過2步完成：首先，確定dN/dS顯著增加的情況（由於正選擇或者選擇約束放鬆）；然後，過濾掉顯著正選擇的情況。
針對於CladeC和CladeD模型，一般要設置幾個不同的起始ω，測試其lnL值是否穩定（ω=0.001，ω=0.01，ω=0.1，ω=0.5，ω=1，ω=1.5），最終一般會選取lnL較大的值作為最終值。
正常分析，應該先用M0估計樹的枝長，Kappa值，然後用跑出來的樹作為初始樹，並設置 fix_blength = 2 。
鏈接1
鏈接2
CladeC經常用於檢測，不同分支的分化選擇壓力，但有時候檢測到前景枝dN/dS>1，此時我們需要進一步利用CladeC的零假設進一步測試（fix_omega=1,omega=1），或者利用branch-site進一步驗證一致性。
free ratio模型估計通常會造成較大的抽樣誤差，例如，較短的分支通常會具有較大的dS/dN。所以，一般對於free ratio出現dN/dS > 999或者dN、dS < 0.0005結果不建議採用。
在dataset中，有三個clade：A, B, C，那麼在clade 1和clade 2之間是否存在顯著差異？
假設，前景枝的標注如下所示：
首先，為了測試clades之間的顯著差異，您可以比較CladeC與M2a_rel。M2a_rel假設 2、$0等都是在相同的選擇壓力下進化的，這個測試應該有兩個自由度。
其次，為了測試clade A和B之間的顯著差異，同時允許clade C是不同的，您可以比較使用上面提供的樹運行CMC與使用更簡單的樹運行CMC的契合程度。在這種情況下，更簡單的樹會將clade A和B分配給同一個組，這個測試應該有一個自由度。如下所示：
當數據集中有多個前景支時：1）進行多個測試，然後在每個測試中設置一個感興趣的分支作為前景分支；2）只進行一次測試，在此測試中中將所有感興趣的分支設置為前景分支。那麼，此時又會出現另一個問題就是進行多個測試時其他感興趣分支是否應該被去除？這也許應該取決於具體的生物學問題。
https://groups.google.com/g/pamlsoftware/c/aVj2opOg7PA
如果校正之後沒有顯著地，我們可以選擇adjP排序
使用free ratio結果可能會產生較大誤差 https://groups.google.com/g/pamlsoftware/c/2DrYS0ff7_o
正選擇位點的狀態是多序列比對中第一條參考序列的狀態，而不是前景支的序列狀態。此外，還有注意cutdata是否設置為1。
https://groups.google.com/g/pamlsoftware/c/ZnPaysiZKbI
『叄』 分子進化樹構建及數據分析方法介紹【轉】
首先是方法的選擇。
基於距離的方法有UPGMA、ME（困虧Minimum Evolution，最小 進化 法）和NJ（Neighbor-Joining，鄰接法）等。其他的幾種方法包括MP（Maximum parsimony，最大簡約法）、ML（Maximum likelihood，最大似然法）以及貝葉斯（Bayesian）推斷等方段桐法。其中UPGMA法已經較少使用。
一般來講，如果模型合適，ML的效果較好。對近緣序列，有人喜歡MP，因為用的假設最少。MP一般不用在遠緣序列上，這時一般用NJ或ML。對相似度很低的序列，NJ往往出現Long-branch attraction（LBA，長枝吸引現象），有時嚴重干擾 進化樹 的構建。貝葉斯的方法則太慢。對於各種方法構建分子 進化樹 的准確性，一篇綜述（Hall BG. Mol Biol Evol 2005, 22(3):792-802）認為貝葉斯的方法最好，其次是ML，然後是MP。其實如果序列的相似性較高，各種方法都會得到不錯的結果，模型間的差別也不大。
對於NJ和ML，是需要選擇模型的。對於各種模型之間的理論上的區別，這里不作深入的探討，可以參看Nei的書。對於蛋白質序列以及DNA序列，兩者模型的選擇是不同的。以作者的經驗來說，對於蛋白質的序列，一般選擇Poisson Correction（泊松修正）這一模型。而對於核酸序列，一般選擇Kimura 2-parameter（Kimura-2參數）模型。如果對各種模型的理解並不深入，作者並不推薦初學者使用其他復雜的模型。
Bootstrap幾乎是一個必須的選項。一般Bootstrap的值>70，則認為構建的 進化樹 較為可靠。如果Bootstrap的值太低，則有可能 進化樹 的拓撲結構有錯誤， 進化樹 是不可靠的。
對於 進化樹 的構建，如果對理論的了解並不深入，作者推薦使用預設的參數。需要選擇模型的時候（例如用NJ或者ML建樹），對於蛋白序列使用Poisson Correction模型，對於核酸序列使用Kimura-2參數模型。另外需要做Bootstrap檢驗，當Bootstrap值過低時，所構建的 進化樹 其拓撲結構可能存在問題。並且，一般推薦用兩種不同的方法構建 進化 樹，如果所得到的 進化 樹類似，則結果較為可靠。
軟體的選擇 表1中列出了一些與構建分子 進化 樹相關的軟體。
構建NJ樹，可以用PHYLIP（寫得有點問題，例如比較慢，並且Bootstrap檢驗不方便）或者MEGA。MEGA是Nei開發的方法並設計的圖形化的軟體，使用非常方便。作者推薦MEGA軟體為初學者的首選。雖然多雪列比對工具ClustalW/X自帶了一個NJ的建樹程序，但是該程序只有p-distance模型，而且構建的樹不夠准確，一般不用來構建 進化 樹。
構建MP樹，最好的工具是PAUP，但該程序屬於商業軟體，並不對學術免費。因此，作者並不建議使用PAUP。而MEGA和PHYLIP也可以用來構建 進化 樹。這里，作者推薦使用MEGA來構建MP樹。理由是，MEGA是圖形化的軟體，使用方便，而PHYLIP則是命令行格式的軟體，使用較為繁瑣。對於近緣序列的進化樹構建，MP方法幾乎是最好的。構建ML樹可以使用PHYML，速度最快。或者使用Tree-puzzle，速度也較快，並且該程序做蛋白質序列的進化樹效果比較好。而PAML則並不適合構建進化樹。
ML的模型選擇是看構出的樹的likelihood值，從參數少，簡單的模型試起，到likelihood值最大為止。ML也可以使用PAUP或者PHYLIP來構建。這里作者推薦的工具是BioEdit。BioEdit集成了一些握尺坦PHYLIP的程序，用來構建進化樹。Tree-puzzle是另外一個不錯的選擇，不過該程序是命令行格式的，需要學習DOS命令。PHYML的不足之處是沒有win32的版本，只有適用於64位的版本，因此不推薦使用。值得注意的是，構建ML樹，不需要事先的多序列比對，而直接使用FASTA格式的序列即可。
貝葉斯的演算法以MrBayes為代表，不過速度較慢。一般的進化樹分析中較少應用。由於該方法需要很多背景的知識，這里不作介紹。
表1 構建分子進化樹相關的軟體
軟體
網址
說明
ClustalX
http://bips.u-strasbg.fr/fr/Documentation/ClustalX/
圖形化的多序列比對工具
ClustalW
http://www.cf.ac.uk/biosi/research/biosoft/Downloads/clustalw.html
命令行格式的多序列比對工具
GeneDoc
http://www.psc.e/biomed/genedoc/
多序列比對結果的美化工具
BioEdit
http://www.mbio.ncsu.e/BioEdit/bioedit.html
序列分析的綜合工具
MEGA
http://www.megasoftware.net/
圖形化、集成的進化分析工具，不包括ML
PAUP
http://paup.csit.fsu.e/
商業軟體，集成的進化分析工具
PHYLIP
http://evolution.genetics.washington.e/phylip.html
免費的、集成的進化分析工具
PHYML
http://atgc.lirmm.fr/phyml/
最快的ML建樹工具
PAML
http://abacus.gene.ucl.ac.uk/software/paml.html
ML建樹工具
Tree-puzzle
http://www.tree-puzzle.de/
較快的ML建樹工具
MrBayes
http://mrbayes.csit.fsu.e/
基於貝葉斯方法的建樹工具
MAC5
http://www.agapow.net/software/mac5/
基於貝葉斯方法的建樹工具
TreeView
http://taxonomy.zoology.gla.ac.uk/rod/treeview.html
進化樹顯示工具
需要注意的幾個問題是:
其一，如果對核酸序列進行分析，並且是CDS編碼區的核酸序列，一般需要將核酸序列分別先翻譯成氨基酸序列，進行比對，然後再對應到核酸序列上。這一流程可以通過MEGA 3.0以後的版本實現。MEGA3現在允許兩條核苷酸，先翻成蛋白序列比對之後再倒回去，做後續計算。
其二，無論是核酸序列還是蛋白序列，一般應當先做成FASTA格式。FASTA格式的序列，第一行由符號「>」開頭，後面跟著序列的名稱，可以自定義，例如user1，protein1等等。將所有的FASTA格式的序列存放在同一個文件中。文件的編輯可用Windows自帶的記事本工具，或者EditPlus（google搜索可得）來操作。
文件格式如圖1所示：
圖1 FASTA格式的序列
NCBI的COG介紹：
什麼是 COG ？
「 COG 」是Cluster of Orthologous Groups of proteins（蛋白相鄰類的聚簇）的縮寫。構成每個 COG 的蛋白都是被假定為來自於一個祖先蛋白，並且因此或者是 orthologs 或者是paralogs。Orthologs是指來自於不同物種的由垂直家系（物種形成）進化而來的蛋白，並且典型的保留與原始蛋白有相同的功能。Paralogs是那些在一定物種中的來源於基因復制的蛋白，可能會進化出新的與原來有關的功能。請參考文獻獲得更多的信息。
COG 分類是如何構建的？
COG 是通過把所有完整測序的基因組的編碼蛋白一個一個的互相比較確定的。在考慮來自一個給定基因組的蛋白時，這種比較將給出每個其他基因組的一個最相似的蛋白（因此需要用完整的基因組來定義 COG 。注1）這些基因的每一個都輪番的被考慮。如果在這些蛋白（或子集）之間一個相互的最佳匹配關系被發現，那麼那些相互的最佳匹配將形成一個 COG （注2）。這樣，一個 COG 中的成員將與這個 COG 中的其他成員比起被比較的基因組中的其他蛋白更相像，盡管如果絕對相似性比較的。最佳匹配原則的使用，沒有了人為選擇的統計切除的限制，這就兼顧了進化慢和進化快的蛋白。然而，還有一個加的限制就是一個COG必須包含來自於3個種系發生上遠的基因組的一個蛋白。
注1：僅僅應用在形成COG時，不包含新蛋白的信息。
注2：為了簡化，許多步驟都省略的，請參考文獻。
使用COG可以得到什麼樣的信息？
簡單的說，有三方面的信息：
1，蛋白的註解。COG的一個蛋白成員的已知功能（以及二維或三維結構）可以直接應用到COG的其他成員上去。然而，這里也要警告，因為有些COG含有paralogs，它們的功能並非對應與那些已知蛋白。
2，種系發生圖譜。這給出在一個特定的COG中一個給定物種是否存在某些蛋白。系統使用，這些圖譜可以用來確定在一個物種中是否一個特定的代謝途徑。
3，多重對齊。每一個COG頁麵包括了一個鏈接到COG成員的一個多重對齊，那可以被用來確定保守序列殘基和分析成員蛋白的進化關系。
COG分類有哪些？
目前COG分類中每個字母代表的功能分類含義：
INFORMATION STORAGE AND PROCESSING
[J] Translation, ribosomal structure and biogenesis
[A] RNA processing and modification
[K] Transcription
[L] Replication, recombination and repair
[B] Chromatin structure and dynamics
CELLULAR PROCESSES AND SIGNALING
[D] Cell cycle control, cell division, chromosome partitioning
[Y] Nuclear structure
[V] Defense mechanisms
[T] Signal transction mechanisms
[M] Cell wall/membrane/envelope biogenesis
[N] Cell motility
[Z] Cytoskeleton
[W] Extracellular structures
[U] Intracellular trafficking, secretion, and vesicular transport
[O] Posttranslational modification, protein turnover, chaperones
METABOLISM
[C] Energy proction and conversion
[G] Carbohydrate transport and metabolism
[E] Amino acid transport and metabolism
[F] Nucleotide transport and metabolism
[H] Coenzyme transport and metabolism
[I] Lipid transport and metabolism
[P] Inorganic ion transport and metabolism
[Q] Secondary metabolites biosynthesis, transport and catabolism
POORLY CHARACTERIZED
[R] General function prediction only
[S] Function unknown
遺傳密碼的新排列和起源探討
肖景發, 於軍 中國科學院北京基因組研究所, 中國科學院「基因組科學及信息」重點實驗室
摘要根據DNA核苷酸組分的動態變化規律將遺傳密碼的傳統排列按 密碼子 對GC和嘌呤含量的敏感性進行了重排. 新密碼表可劃分為兩個半區(或1/2區)和四個四分區(或1/4區). 就原核生物基因組而言, 當 GC含量 增加時, 物種蛋白質組所含的氨基酸傾向於使用GC富集區和嘌呤不敏感半區所編碼的氨基酸, 它們均使用四重簡並密碼, 對DNA序列的突變具有相對魯棒性(Robustness). 當 GC含量 降低時, 大多數 密碼子 處於AU富集區和嘌呤敏感半區, 這個區域編碼的氨基酸具有物理化學性質的多樣性. 因為當 密碼子 第三位核苷酸(CP3)在嘌呤和嘧啶之間發生轉換時, 密碼子 所編碼的氨基酸也傾向於發生變化.
關於遺傳密碼的 進化 存在多種假說, 包括凝固事件假說、共 進化 假說和立體化學假說等, 每種假說均試圖解釋遺傳密碼所表現出來的某些化學和生物學規律. 基於遺傳密碼的物理化學性質、基因組變異的規律和相關的生物學假說, 我們提出了遺傳密碼 分步進化假說 (The Stepwise Evolution Hypothesis for the Genetic Code). 在人們推斷的最原始的RNA世界裡, 原初(Primordial)遺傳密碼從只能識別嘌呤和嘧啶開始, 編碼一個或兩個簡單而功能明確的氨基酸. 由於胞嘧啶C的化學不穩定性, 最初形成的遺傳密碼應該僅僅由腺嘌呤A和尿嘧啶U來編碼, 卻可得到一組七個多元化的氨基酸. 隨著生命復雜性的增加, 鳥嘌呤G從主載操作信號的功能中釋放出來, 再伴隨著C的引入, 使遺傳密碼逐步擴展到12、15和20個氨基酸, 最終完成全部 進化 步驟.
遺傳密碼的 進化 過程同時也伴隨以蛋白質為主體的分子機制和細胞過程的 進化 , 包括氨醯tRNA合成酶(AARS)從初始翻譯機器上的脫離、DNA作為信息載體而取代RNA以及AARS和tRNA共 進化 等基本過程. 分子機制和細胞過程是生命的基本組成元件, 它們不但自己不斷地趨於完善, 也促使生命體走著不盡相同的道路, 要麼維持魯棒性(Robustness, 如細菌), 要麼尋覓多元化(Diversity, 如節肢動物和植物), 要麼追求綜合性(Complexity, 如脊椎動物).
自從 密碼子 被全部發現以來, Crick[1]
就將遺傳密碼表排列成化學家所認可的形式. 盡管後來有些特殊表現形式的列方式(如同心圓、八卦式和二元密碼等), 但其基本排布一直延續至今[1~3]
. 遺傳密碼以4個脫氧核糖核苷酸作為基本符號來組成遺傳信息, 並以20個氨基酸作為基本結構單元來構建蛋白質. 遺傳密碼是使用4個鹼基(兩個嘌呤: 腺嘌呤A和鳥嘌呤G; 2個嘧啶: 尿嘧啶U和胞嘧啶C)構成的三聯體 密碼子 , 共64個, 分別對應20個氨基酸或翻譯起始和終止信號. 生物體要將DNA分子中儲存的信息內涵轉變成功能內涵, 就要利用信使mRNA、解碼分子tRNA和完整翻譯機器等多重功能. 各種復雜分子機制和細胞過程的誕生和成熟一定會反映生命從RNA世界到RNA-蛋白質世界, 再到RNA-蛋白質-DNA世界逐漸轉變的過程, 遺傳密碼作為一個獨立的生物學機制也一定是漫長生命 進化 過程中的一個必然產物.
20世紀60年代初, 實驗分子生物學最大的進展就是解碼遺傳密碼, 發現它在生命有機體中, 基本是統一的. 自此不同的假設均試圖解釋遺傳密碼的信息和化學特性, 從簡單的凝固事件假說到更復雜的統計學、共 進化 和立體化學理論. 凝固事件假說認為 密碼子 與氨基酸的對應關系是在某個生命發生時段里被固定下來, 並且很難被改變[2]
, 這個假說一直被基於適應性、歷史性和化學性的不同論點所挑戰[4]
. 盡管關於遺傳密碼的 進化 也有人提出過不同的假設, 但是解釋 密碼子 的分配原則、物理化學性質的相關性和DNA組分變化對 密碼子 使用頻率的牽動, 從而揭示遺傳密碼表的生物學本質仍然是一個不小的挑戰[2,5]
.
1 重排遺傳密碼表
重排遺傳密碼表有3個重要原因. 首先, DNA序列有4個最基本的可度量的變化, 即核苷酸序列、序列長度、 GC含量 和嘌呤(R或AG)含量. 假如把核苷酸序列和長度相對於時間的變化暫時不考慮, 那麼只有後面的兩個變數對於傳統的遺傳密碼表具有影響力, 所以重排應該以GC和嘌呤含量的變化為主線. 但以前大家熟知的密碼表排列只是為了簡明和清晰地顯示 密碼子 和氨基酸的一一對應關系, 卻忽略了密碼表本身對氨基酸物理化學性質多樣性的表現和DNA編碼承受突變的魯棒性等明顯信息. 因此, 有必要把傳統的密碼表進行重新排列[6]
, 使其能夠表現信息內涵和功能內涵之間的基本關系. 其次, 當 GC含量 和嘌呤含量變化時, 希望從密碼表中找出相應蛋白質組成變化的線索. 圖1展示了 GC含量 和嘌呤含量在極端狀態下4個微生物基因組的氨基酸組分分布.
(1) RNA世界和早期遺傳密碼. RNA世界的存在首先被RNA分子具有相應催化功能的生物學特性所支持[22~26]
. 在RNA世界裡, RNA具有雙重的功能, 既是信息載體也是功能載體. 因為生命的基本分子機制和細胞過程起源於RNA世界, 所以沒有理由說遺傳密碼不起源於RNA世界. 在RNA世界裡RNA分子可以組成簡單的核苷酸多聚物, 這種多聚物在近億年的成熟期里, 為生命提供了足夠的功能上的復雜性和多樣性. 原始細胞可以通過相互爭斗和吞噬獲得基本的組成成分，因此基於模板的RNA合成可能對於生命的初始不是必需的. 可以想像這些RNA分子可以通過簡單的聚合酶來合成, 通過自身剪接或化學修飾轉變為其他相似的結構, 從而達到結構的可變性和功能的多樣性. 此外, RNA的編輯(RNA Editing)也一定起了非常重要的作用, 這一分子機制一直延續到現在, 在包括人類在內的高級物種中仍然存在.
在現代生物世界裡, 剪接體(Spliceosome)通常是用於RNA分子的剪接, 由蛋白質和RNA分子組成. 可以做兩個假設, 生命可能起源於類真核有機體的原型細胞(在DNA引入之前)而不是類原核有機體的原型細胞. 在RNA組成的翻譯機器(Translational Machinery)沒有形成之前, 初始遺傳密碼可能不是必需的. 一旦這個初始生命進入到RNA-蛋白質組成的世界時, 多肽才逐漸按照密碼子開始有序合成, 遺傳密碼就開始發揮其作用了. 可以認為有序的生命可以在與相對無序生命的爭斗中更容易獲勝和取得繁衍的空間.
現在可以推測初始遺傳密碼在RNA世界存在和 進化 的基本過程和起源時的基本邏輯關系. 初始生命一定比較簡單, 分子間相互作用也比較寬松, 最小的編碼系統可能只要區分嘌呤R和嘧啶Y就夠了. 假定現代密碼在生命的早期階段已經被統一並相對忠實地繼承了RNA密碼的基本關系, 這個可能的原始編碼就至少有7個氨基酸(I和M視為等同; 圖5), 同時也有起始和終止密碼子. 這7個氨基酸的側鏈具有廣泛的物理化學性質(氨基、醯基、苯環、羥基、酚基、烴鏈和甲硫基等), 但是沒有小的和酸性的氨基酸. 可以推測: 體積小的氨基酸在初始蛋白質相互作用中的作用顯然不如大的重要, 而鹼性氨基酸的功能對於酸性DNA則是顯而易見的. 另外的一種可能性是氨基酸與tRNA以及AARS之間的關系不是十分明確, 一個密碼子對應多氨基酸的情況可能在遺傳密碼成熟前是普遍存在的[27]
. 由於7個氨基酸的編碼區處於現代密碼表的AU富集區, 可以確信初始密碼子始於這個區域, 後來擴展到嘌呤敏感區即所謂趨變半區. 這個階段的存在既復合由簡到繁的邏輯, 也迎合了實驗的證據, 那就是C的不穩定性和G在RNA操作功能上的作用[21, 28~29]
.
(3) 遺傳密碼的第二次拓展. 當GU和AG從作為剪接信號功能釋放出來以後(剪接體的結構和功能隨著蛋白質的演變而復雜化和精密化), 遺傳密碼引入了Arg, Ser和Val. 氨基酸的個數變成15個, 這次擴展是對已經存在的氨基酸物理化學性質和二級結構特性的擴展. Arg是Lys的替代體, Ser則對應Tyr, Val是疏水性氨基酸Leu, Ile和Met的補充[32~35]
.
最具吸引力的是六重簡並的3個氨基酸Arg, Leu和Ser. 這些氨基酸在被引入後, 又由於核苷酸C在RNA世界的應用而擴展出各自的新四聯碼, 成為六重簡並. 首先, Leu是在現代基因組中包括所有三界生物在內最豐富的氨基酸, Ser是真核生物第二豐富的氨基酸, Arg也是一個富有的氨基酸, 通常在細菌基因組中位於前10位. 其次, Leu在二重簡並密碼和四重簡並之間最容易轉換, 只需要通過簡單U到C轉換(UUR-CUR)即可, 這也說明Leu對於大多數蛋白質來說是用於當 GC含量 增加時維持蛋白質功能的完整性. 這些觀察引出相應的假設: 這3個氨基酸的附加密碼是為了當 GC含量 或AG含量增加時平衡富有氨基酸, 相應的密碼分布按照平衡遺傳密碼的蛋白質多樣性和蛋白質魯棒性二等分. 這種平衡能力用於當編碼序列突變發生時穩定蛋白質的氨基酸組成, 從而維護蛋白質結構的完整性.
(4) 遺傳密碼的最終拓展. 遺傳密碼的最終拓展是在DNA作為信息載體取代RNA使得信息載體具有更高的准確性和穩定性, 同時也產生了最為關鍵的從RNA到DNA的逆轉錄機制. 基於模板的DNA復制機制開辟了新的DNA-蛋白質-RNA世界. 很多新分子機制的 進化 包括DNA復制和修復、RNA的轉錄等, 使這個生物界里分子機制和細胞過程更趨於多元和完善. 同時當C和其脫氧衍生物分別作為結構模板加入RNA和DNA時, 標准遺傳密碼也就隨之產生並被固定下來. 遺傳密碼本身得到新的補充並且編碼能力有了很大提高. 組氨酸(His)和Glu立刻加入進來, 主要是由於它們具有相應的催化性質以及和原有的兩個鹼性氨基酸的相似性, Thr擴展了Ser的功能, 同時使蛋白質的結構增加了精細度, Ala同Ser相比具有類似的體積和尺度, 但其和Ser比具有很強的疏水性質[32,33]
. 這些新引入的氨基酸在蛋白質結構和功能多樣性上起到非常關鍵的作用. 不容懷疑的是Pro的最後加入, 它具有其他氨基酸所不具備的性質, 即通過特有的方式使蛋白質的骨架結構扭曲達到蛋白質結構的緊密折疊. 相應的擴展模式在AARS同樣得到支持遺傳密碼擴展的假設, 除了3個六重簡並的遺傳密碼外, 這次共有六組遺傳密碼最終被引入, 同時編碼6個氨基酸. 這6個氨基酸的AARS分類按照G和I 的配對原則延伸而來. 例如AARS對於雙重編碼的氨基酸His(CAR)和Gln(CAY)的對應, Glu(GAR)和Asp(GAY)的對應等.
遺傳密碼的 進化 就是密碼子的有序發生和合理分布, 這個分布的合理性一定經過一個復雜選擇過程. 首先, 通過長時間的創造和優化, 使其在基因組核苷酸序列發生突變時對蛋白質的結構起到緩沖的作用; 第二, 密碼子採取這樣一種特殊的排布方式: 當DNA組成從AU富集區到GC富集區改變時, 氨基酸的分布傾向於從具有催化性質的氨基酸轉到具有結構性質的氨基酸; 第三, 充分利用密碼子第三位多變的優勢(通常體現在R和Y之間的轉換), 來改變編碼氨基酸的物理化學性質, 致使在趨變半區里大約有15個氨基酸對第三個位置R和Y之間的轉換呈現敏感.
(5) 分子機制與細胞過程的 進化 . 盡管分子機制與細胞過程的根本界限有時會很模糊, 但還是將它們分開: 前者強調物理性的相互作用、發生的空間和組分的存在, 後者強調化學反應的結果、發生的時間和過程. 從一方面講, DNA的變異顯然是細胞過程的產物, 遺傳密碼的發生和最終形成也是它的產物. 從另一方面講, 密碼子與氨基酸的關系影響到細胞的蛋白質組分的變化, 即分子機制的變化[21]
. 比如, 如果在RNA世界需要產生多個拷貝的RNA分子, 一定需要一個分子機制來實現. 在現代生物世界裡, 通常是由以DNA為模板的轉錄機制來完成, 但在RNA世界裡沒有RNA的復制, 多個RNA分子產生是由多聚酶和編輯體(Editosome)共同來完成的. 也許就是那個最原始的細胞機制. RNA世界的第二個分子機制發明可能是就剪接體, 這個分子機制在現代生物世界裡仍然在發揮其重要的作用. 第三個分子機制也許是翻譯體(Translatosome)的形成, 其用於直接進行蛋白質分子的加工, 這一分子機制是從原始的RNA世界到成熟的RNA世界再到現代生物世界裡轉折的重要標志. 在轉折期里, 分子機制在蛋白質精確度的變化中不斷完善和復雜, 直到DNA通過RNA和蛋白質的復合體引進到生命世界
『肆』 選擇壓力分析，從序列下載到結果繪圖
選擇壓力及基因在進化過程中所受到的壓力，也稱自然選擇。
具體知識我不作贅述，感興趣的可以去看其他文章，本文主要介紹選擇壓力分析具體的冊仔流程晌談及方法。這方面網上的教程很少，我也只是在各種網上的教程里自己摸索整理出的方法，方法不一定對，僅供參考。宴姿碰
然後對序列最長CDS序列進行提取，並將物種名稱加到ID的前面，得到如下格式：
這里使用 OrthoFinder尋找同源基因並建樹 ,得到直系單拷貝同源基因和物種系統發育樹結果。
目前我們有了選擇壓力分析所必須的文件：
選擇壓力分析使用的是Paml中的codeml軟體包。 EasyCodeML 做了很好的可視化界面，但是不能批量運算，只能一次分析一個，並且容易卡死（Mac 和 windows上都遇到過）。
所以我自己寫了一個腳本 callCodeml （使用方法： python3 callCodeml.py 序列文件夾 物種發育樹文件 ）來計算前景枝受到的選擇壓力，調用codeml批量運算 枝位點模型（Branch site model） ，並使用chi2檢驗，只輸出 P < 0.05 的結果，最終將結果自動整理成表格。
至此，選擇壓力分析完成，我們只需要根據OrthoFinder的文件即可找到OG00XXXX同源組所對應的基因名稱。
困了，睡覺，改天有緣再寫。
『伍』 進化樹構建需要哪些工具
序列比對建議用ClustalX
建NJ或MP樹，用MEGA就可以了，非常方便
若要建ML樹推薦用phyML
建Bayes樹推薦用Parallel MrBayes @ BioHPC
如果不是專業建樹的話，MEGA足夠用了，建議參考下面這篇文章：
一、引言
開始動筆寫這篇短文之前，我問自己，為什麼要寫這樣的文章？寫這樣的文章有實際的意義嗎？我希望能夠解決什麼樣的問題？帶著這樣的疑惑，我隨手在丁香園（DXY）上以關鍵字「進化 分析 求助」進行了搜索，居然有289篇相關的帖子（2006年9月12日）。而以關鍵字「進化分析」和「進化」為關鍵字搜索，分別找到2,733和7,724篇相關的帖子。考慮到有些帖子的內容與分子進化無關，這里我保守的估計，大約有 3,000~4,000篇帖子的內容，是關於分子進化的。粗略地歸納一下，我大致將提出的問題分為下述的幾類：
1．涉及基本概念
例如，「分子進化與生物進化是不是一個概念」，「關於微衛星進化模型有沒有什麼新的進展」以及「關於Kruglyak的模型有沒有改進的出現」，等等。
2．關於構建進化樹的嫌跡畢方法的選擇
例如，「用boostrap NJ得到XX圖，請問該怎樣理解？能否應用於文章？用boostrap test中的ME法得到的是XXX樹，請問與上個樹比，哪個更好」，等等。
3．關於軟體的選擇
例如，「想做一個進化樹，不知道什麼軟體能更好的使用且可以說明問題，並且有沒有說明如何做」，「拿到了16sr RNA數據，打算做一個系統進化樹分析，可是原來沒有做過這方面的工作啊，都要什麼軟體」，「請問各位高手用ClustalX做出來的進化樹與 phylip做的有什麼區別」，「請問有做過進化樹分析的朋友，能不能提供一下，做樹的時候參數的設置，以及代表的意思。還有各個分支等數值的意思，說明的問題等」，等等。
4．蛋白家族的分類問題
例如，「搜集所有的關於一個特定domain的序列，共141條，做的進化樹不知具體怎麼分析」，等等。
5．新基因功能的推斷
例如，「根據一個新基因A氨基酸序列構建的系統發生樹，這個進化樹能否說明這個新基因A和B同源，屬於同一基因家族」，等等。
6．計算基因分化的年代
例如，「想在基因組水平比較兩個或三個比較接近物種之間的進化年代的遠近，具體推算出他們之間的分歧時間」，「如何估計病毒進化中變異所需時間」，等等。芹芹
7．進化樹的編輯
例如生成的進化樹圖片，如何進行後續的編輯，比如希望在圖片上標注某些特定的內容，等等。
由於相關的帖子太多，作者在這里對無法閱讀全部的相關內容而致以歉意。同時，作者歸納的這七個問題也並不完全代表所有的提問。對於問題1所涉及到的基本的概念，作者推薦讀者可參考由Masatoshi Nei與Sudhir Kumar所撰寫的《分子進化與系統發育》（Molecular Evolution and Phylogenetics）一書，以及相關的分子進化方面的最新文州簡獻。對於問題7，作者之一lylover一般使用Powerpoint進行編輯，而 Photoshop、Illustrator及Windows自帶的畫圖工具等都可以使用。
這里，作者在這里對問題2-6進行簡要地解釋和討論，並希望能夠初步地解答初學者的一些疑問。
二、方法的選擇
首先是方法的選擇。基於距離的方法有UPGMA、ME（Minimum Evolution，最小進化法）和NJ（Neighbor-Joining，鄰接法）等。其他的幾種方法包括MP（Maximum parsimony，最大簡約法）、ML（Maximum likelihood，最大似然法）以及貝葉斯（Bayesian）推斷等方法。其中UPGMA法已經較少使用。
一般來講，如果模型合適，ML的效果較好。對近緣序列，有人喜歡MP，因為用的假設最少。MP一般不用在遠緣序列上，這時一般用NJ或ML。對相似度很低的序列，NJ往往出現Long-branch attraction（LBA，長枝吸引現象），有時嚴重干擾進化樹的構建。貝葉斯的方法則太慢。對於各種方法構建分子進化樹的准確性，一篇綜述（Hall BG. Mol Biol Evol 2005, 22(3):792-802）認為貝葉斯的方法最好，其次是ML，然後是MP。其實如果序列的相似性較高，各種方法都會得到不錯的結果，模型間的差別也不大。
對於NJ和ML，是需要選擇模型的。對於各種模型之間的理論上的區別，這里不作深入的探討，可以參看Nei的書。對於蛋白質序列以及DNA序列，兩者模型的選擇是不同的。以作者的經驗來說，對於蛋白質的序列，一般選擇Poisson Correction（泊松修正）這一模型。而對於核酸序列，一般選擇Kimura 2-parameter（Kimura-2參數）模型。如果對各種模型的理解並不深入，作者並不推薦初學者使用其他復雜的模型。
Bootstrap幾乎是一個必須的選項。一般Bootstrap的值>70，則認為構建的進化樹較為可靠。如果Bootstrap的值太低，則有可能進化樹的拓撲結構有錯誤，進化樹是不可靠的。
對於進化樹的構建，如果對理論的了解並不深入，作者推薦使用預設的參數。需要選擇模型的時候（例如用NJ或者ML建樹），對於蛋白序列使用Poisson Correction模型，對於核酸序列使用Kimura-2參數模型。另外需要做Bootstrap檢驗，當Bootstrap值過低時，所構建的進化樹其拓撲結構可能存在問題。並且，一般推薦用兩種不同的方法構建進化樹，如果所得到的進化樹類似，則結果較為可靠。
三、軟體的選擇
表1中列出了一些與構建分子進化樹相關的軟體。
構建NJ樹，可以用PHYLIP（寫得有點問題，例如比較慢，並且Bootstrap檢驗不方便）或者MEGA。MEGA是Nei開發的方法並設計的圖形化的軟體，使用非常方便。作者推薦MEGA軟體為初學者的首選。雖然多雪列比對工具ClustalW/X自帶了一個NJ的建樹程序，但是該程序只有p- distance模型，而且構建的樹不夠准確，一般不用來構建進化樹。
構建MP樹，最好的工具是PAUP，但該程序屬於商業軟體，並不對學術免費。因此，作者並不建議使用PAUP。而MEGA和PHYLIP也可以用來構建進化樹。這里，作者推薦使用MEGA來構建MP樹。理由是，MEGA是圖形化的軟體，使用方便，而PHYLIP則是命令行格式的軟體，使用較為繁瑣。對於近緣序列的進化樹構建，MP方法幾乎是最好的。
構建ML樹可以使用PHYML，速度最快。或者使用Tree-puzzle，速度也較快，並且該程序做蛋白質序列的進化樹效果比較好。而PAML則並不適合構建進化樹。ML的模型選擇是看構出的樹的likelihood值，從參數少，簡單的模型試起，到likelihood值最大為止。ML也可以使用 PAUP或者PHYLIP來構建。這里作者推薦的工具是BioEdit。BioEdit集成了一些PHYLIP的程序，用來構建進化樹。Tree- puzzle是另外一個不錯的選擇，不過該程序是命令行格式的，需要學習DOS命令。PHYML的不足之處是沒有win32的版本，只有適用於64位的版本，因此不推薦使用。值得注意的是，構建ML樹，不需要事先的多序列比對，而直接使用FASTA格式的序列即可。
貝葉斯的演算法以MrBayes為代表，不過速度較慢。一般的進化樹分析中較少應用。由於該方法需要很多背景的知識，這里不作介紹。
表1 構建分子進化樹相關的軟體
軟體 網址 說明
ClustalX http://bips.u-strasbg.fr/fr/Documentation/ClustalX/ 圖形化的多序列比對工具
ClustalW http://www.cf.ac.uk/biosi/research/biosoft/Downloads/clustalw.html 命令行格式的多序列比對工具
GeneDoc http://www.psc.e/biomed/genedoc/ 多序列比對結果的美化工具（可以導入fasta格式的文件，出來的圖可用於發表，我用過）
BioEdit http://www.mbio.ncsu.e/BioEdit/bioedit.html 序列分析的綜合工具
MEGA http://www.megasoftware.net/ 圖形化、集成的進化分析工具，不包括ML
PAUP http://paup.csit.fsu.e/ 商業軟體，集成的進化分析工具
PHYLIP http://evolution.genetics.washington.e/phylip.html 免費的、集成的進化分析工具
PHYML http://atgc.lirmm.fr/phyml/ 最快的ML建樹工具
PAML http://abacus.gene.ucl.ac.uk/software/paml.html ML建樹工具
Tree-puzzle http://www.tree-puzzle.de/ 較快的ML建樹工具
MrBayes http://mrbayes.csit.fsu.e/ 基於貝葉斯方法的建樹工具
MAC5 http://www.agapow.net/software/mac5/ 基於貝葉斯方法的建樹工具
TreeView http://taxonomy.zoology.gla.ac.uk/rod/treeview.html 進化樹顯示工具
（加紅色標注的為最通用的分析軟體）
需要注意的幾個問題是，其一，如果對核酸序列進行分析，並且是CDS編碼區的核酸序列，一般需要將核酸序列分別先翻譯成氨基酸序列，進行比對，然後再對應到核酸序列上。這一流程可以通過MEGA 3.0以後的版本實現。MEGA3現在允許兩條核苷酸，先翻成蛋白序列比對之後再倒回去，做後續計算。
其二，無論是核酸序列還是蛋白序列，一般應當先做成 FASTA格式。FASTA格式的序列，第一行由符號「>」開頭，後面跟著序列的名稱，可以自定義，例如user1，protein1等等。將所有的FASTA格式的序列存放在同一個文件中。文件的編輯可用Windows自帶的記事本工具，或者EditPlus（google搜索可得）來操作。
另外，構建NJ或者MP樹需要先將序列做多序列比對的處理。作者推薦使用ClustalX進行多序列比對的分析。多序列比對的結果有時需要後續處理並應用於文章中，這里作者推薦使用GeneDoc工具。而構建ML樹則不需要預先的多序列比對。
因此，作者推薦的軟體組合為：MEGA + ClustalX + GeneDoc + BioEdit。
四、數據分析及結果推斷
一般碰到的幾類問題是，（1）推斷基因/蛋白的功能；（2）基因/蛋白家族分類；（3）計算基因分化的年代。關於這方面的文獻非常多，這里作者僅做簡要的介紹。
推斷基因/蛋白的功能，一般先用Blast工具搜索同一物種中與不同物種的同源序列，這包括直向同源物（ortholog）和旁系同源物（paralog）。如何界定這兩種同源物，網上有很多詳細的介紹，這里不作討論。然後得到這些同源物的序列，做成FASTA格式的文件。一般通過NJ構建進化樹，並且進行Bootstrap分析所得到的結果已足夠。如果序列近緣，可以再使用MP構建進化樹，進行比較。如果序列較遠源，則可以做ML樹比較。使用兩種方法得到的樹，如果差別不大，並且Bootstrap總體較高，則得到的進化樹較為可靠。
基因/蛋白家族分類。這方面可以細分為兩個問題。一是對一個大的家族進行分類，另一個就是將特定的一個或多個基因/蛋白定位到已知的大的家族上，看看屬於哪個亞家族。例如，對驅動蛋白（kinesin）超家族進行分類，屬於第一個問題。而假如得到一個新的驅動蛋白的序列，想分析該序列究竟屬於驅動蛋白超家族的14個亞家族中的哪一個，則屬於後一個問題。這里，一般不推薦使用MP的方法。大多數的基因/蛋白家族起源較早，序列分化程度較大，相互之間較為遠源。這里一般使用NJ、ME或者ML的方法。
計算基因分化的年代。這個一般需要知道物種的核苷酸替代率。常見物種的核苷酸替代率需要查找相關的文獻。這里不作過多的介紹。一般對於這樣的問題，序列多數是近緣的，選擇NJ或者MP即可。
如果使用MEGA進行分析，選項中有一項是「Gaps/Missing Data」，一般選擇「Pairwise Deletion」。其他多數的選項保持預設的參數。
五、總結
在實用中，只要方法、模型合理，建出的樹都有意義，可以任意選擇自己認為好一個。最重要的問題是：你需要解決什麼樣的問題？如果分析的結果能夠解決你現有的問題，那麼，這樣的分析足夠了。因此，在做進化分析前，可能需要很好的考慮一下自己的問題所在，這樣所作的分析才有針對性。
六、致謝
本文由mediocrebeing在2005年9月8日所發起的討論《關於建樹的經驗》擴充、修改而來。文章的作者按原貼ID出現先後排名，由 lylover執筆。作者同時感謝所有參與討論的戰友。作者lylover感謝中國科大細胞動力學實驗室的金長江博士所給的一些有益的建議。
來源:丁香園(mediocrebeing, rodger, lylover , klaus, oldfish, yzwpf)
『陸』 使用R語言將fasta轉phylip格式
使用paml軟體里的mcmctree功能需要phylip格式的比對結果,找了以下,發敗漏耐現用R包phylotools轉化最方便，搜銀另外提醒一下要用mcmctree的朋友，我發現mcmctree的序列名字長度不能超過50個字元，如果超出的話，建議先改名字再轉化。
先安裝依賴的包（要先裝BioConctor）：
假設fasta文件名為: aligned_fasta.fasta 讀取fasta文件，轉化：
結果文件為out.phy
注意:生成out.phy里,第一列序列名和第二列序列只有一個空格,而mcmctree要求兩個以上,所以需要察春人工加多一個空格.
也可以用shell命令生成:
『柒』 重新認識TBtools，減少你的生信分析煩惱（20190422）
或許，沒有人知道TBtools到底是什麼？能幹啥。
但是看完這個推文，或許你就知道了其中的一部分。
TBtools對外開放兩年多，不時會有熟悉的不熟悉的人與我聊到TBtools。TBtools在每個人的認知上，或許都不一樣。
有的人覺得TBtools就是 耽誤了他們所謂的生物信息學習 。
有的人說，TBtools 保了他們畢業 。
也有的人說，TBtools 幫助了他發了文章
更或者，TBtools。。。
三年前開始寫TBtools。功能很簡單，主要是做做序列提取，也做了BlastWrapper（當時並不穩健）。目的很純粹，課題組的人不會找我提取序列，也可以直接Blast到轉錄組找序列。那時候接觸TBtools的朋友，或許都這么認為。
當然，後來我對這方面進行了各種增強，也保證了其現在的穩健性。 一個輸入窗口，支持不同的輸入 ，無論是提取序列全長，還是提取序列區段；不僅支持ID提取，還支持ID子串匹配...
此外，也增加了基於gff3進行序列提取的功能，比如提取所有序列的全長CDS，全長EXON，甚至是 可以批量一次性提取一個物閉模種的所有啟動子序列。
除了Blast Wrapper，可能還需要調用外部程序的，那麼是muscle（主要是我實現的NW演算法運行效率一般，用到多序列比對，就不提了）
可以從一個子菜單看到
功能相對豐富，包括
無論是兩條序列的直接比較，還是兩個序列文件，甚至是兩個基因組的兩個指定區間的Blast，TBtools中都已經提供了GUI；不僅於此，四種可視化方式，常常能滿足大多數人的需求。
比如
後來，由於Blast2GO太慢了。基於IDmapping的邏輯，我大體寫了一個GO注釋的功能。當然，更重要的或許是直接寫了GO和KEGG富集分析的功能。所以後來，也有不少做非模式生物的朋友認為，TBtools事實就做這個事情。其中也包括一些可視化，比如GO Level2的可視化。後來，我也寫了一些富集結果的可視化。
慢慢地，我發現，網頁版工具，如Venny，明明是很小的韋恩圖繪制功能，轎察緩網路太差，等待緩沖總是佔用了我太多的時間。應該本地化。所以我索性寫了一個 最高支持六組的Venn圖 工具。當然，也沒談有後來的UpsetPlot工具。
基因展示在染色體上的，類似MapChart的工具等。
此外也由於一些工具，如熱圖繪制上，我覺得用起來真的不順手。或者參數太少，或者不容易調整各種細節。所以我也寫了熱圖工具。
所以，或許確實有的朋友就覺得，你這工具，就是一個畫圖工具包。
Venn圖與UpSetPlot
當然，你還可以直接繪制SeqLogo
基於前面我發送過的推文，總的來說，有了TBtools， 所有人無需任何一行命令，也不需要Linux或者虛擬機操作，可以完成常見的基因家族分析。
大體包括的工具有：
似乎， 有一些培訓機構提供的線上線下的基因家族分析培訓，需要使用的各種虛擬機，Linux，命令，腳本 ，統統都可以扔掉。常見的基因家族分析項目，可能只需要TBtools就完全足夠了。
為此，慢慢地開始有朋友給TBtools下了定義： TBtools是一個基因家族分析工具包 。在我看來，事實上，這些朋友對TBtools有很大的誤解。我從來就沒想過寫一個基因家族分析工具。不是因為TBtools要做一個基因家族分析工具，而是因為基因家族分析本身就是 所有人都需要，都懂的生物數據分析的基本技能 。我只是簡化這些技能的實現。
正如下面，我所寫的TBtools中，或許很多人想想不到功能一樣。
與其說TBtools是你以為的基因家族分析工具，你不如說他是比較基因組分析工具，那麼還顯得高大上一些。
這半個月以來，課題組的安排下，我參與了一些基因組分析相關的工作；
基因組分析，本身確實是一個 耗費智力和體力 的活。分析過程中，也發現了一些或許可以 讓所有人都從中獲取生物信息 的分析手段。
於是我 花了兩個晚上 的時間，寫了幾個工具。
加上一些以前TBtools中就有的工具，相信對做 基因組分析 的朋友會有一定的幫助。
但是，請注意，除非是贊助我們課題組的戶外拓展活動或者合作單位，否則我並不保證這些工具的使用。
MCScanX是比較基因組分析中常用的工具。我已經將其打包到TBtools中，所以即使是windows用戶，也可以輕松進行分析。此外，也不要求用戶保證gff文件和blast文件的名字一致。
關於這個工具... 詳細見公眾號以前的推文。
Ka/Ks的計算，常常會被人問題。事實上，如果只是簡單的進行 NG 演算法的計算，是非常容易實現的。目前用的廣泛的，或許是KaksCalculator2和PAML。這兩個軟體都是大牛級的軟體。在TBtools中，我終於還是開放了去年還是前年實現的NG計算邏輯，並打了非常方便的GUI。用戶幾乎只有有CDS序列和基因對信息就可以直接進行計算，而完全不用浪費時間在文件格式整理上。
工具是不斷地優化和發展的。
分析門檻也是會被不斷打破的。
或許讓所有人，都能有開展一些分析的能力，也是推進一些事物發展的方式。
歡迎尚未加到TBtools使用交流QQ群的朋友，加入
課題組每年暑期有 內部生信入門培訓 ，主要是對實驗室新生開展（以及濕實驗為主的成員）培訓。一直有收到其他課題組想要了解 我們課題組生信數據分析 的想法。故，在博導的提議和課題組的討論後，我們近期計劃，在本年度暑期（7~8月份之間） 對外 增設 生信基礎培訓 名額10枚（前面每年只是課題組內培訓，而不對外）。具體請見 https://mp.weixin.qq.com/s/OtmeTErd9f9rvjJPtBKjMw
園藝植物小分子RNA與基因組研究 -夏瑞課題組
課題組主頁： http://xialab.scau.e.cn/
『捌』 Hyphy，不亞於Paml的選擇壓力分析的優秀軟體，使用指北
近幾年來，Hyphy的使用人數越來越來多，雖然不及paml，但這款軟體的一些優秀特性使得它值得受到使用和關注。
首先相比paml，hyphy有以下幾大優點：
接下來介紹的一系列東西，實際上是對Hyphy官方網站的一系列教程的總結，很多東西官網都寫得很清楚，官網地址為 https://www.hyphy.org/ 。
如果你不想看長篇大論，直接跳到最後的總結部分，那裡有最簡練的總結
關於Hyphy的不同版本 ,hyphy的網頁版即是datamonkey，並且還有GUI版本，這里介紹的主要是命令行版本，並且命令行版本也可以分為互動式運行和一行命令運行，這里不介紹互動式方法的使用。
關於hyphy的安裝 ,只需要用conda就可以安裝了
關於hyphy的輸老畝芹入文件 ，要求一顆newick格式（只能是此格式）系統發育樹以及相對應的fasta序列比耐頃對文件（可以是FASTA, phylip, 等等），標注foreground branch，即前景支的方法和paml略微不同，即在newick文件中在分支名和支長（如果有的話）之間加上{Foreground}來標注，或者你可以去hyphy官網的phylotree來在線標注，地址為 http://phylotree.hyphy.org/ 。
關於多線程支持方面 ，2.4.0版本當中，軟體中的命令 hyphy 和 HYPHYMPI 已經等同，都是調用多線程，在這個版本之前， hyphy 調用的是單核，而 HYPHYMPI 則對應的是多線程版本命令。
關於具體使用方法 ，hyphy的使用非常簡單，以2.2.0版本為例，如果你要使用多線程命令，則如下有兩種方法，分別對應位點模型slac以及指定了前景支的支位點模型absrel,它們都需要openmp支持
不同的模型，只需要改相應的模型名稱就可以調用了(替換上面命令的slac或者absrel)，用法非常簡單,如果不特別用branches指定Foreground，那麼則會默認對整個系統發育應用模型。
關於輸出結果及結果可視化 ，Hyphy運行的時候，默認列印到屏幕上的結果是以markdown格式輸出的，這個結果還是很直觀的，而保存到本地文件的結果是以json格式輸出的，並不是很直觀( 但json格式可以很方便的用python的json模塊提取各種信息，例如pvalue和正選擇位點，在多個任務批量操作的時候，非常的方便，這種保存的格式非常具有通用性，其實是件好事 )，默認是輸出到和多序列比對文件相同的文件夾，可以用 --output 來改變輸出位置，可以去官網 http://vision.hyphy.org/ 來可視化輸出結果，具體的格式介紹，詳見 https://www.hyphy.org/resources/json-fields.pdf 。
關於Hyphy的各種模型，基本上都可以分為不指定foreground和指定foreground運行的兩種方式，前者對應的是檢測 pervasive (across the whole phylogeny) positive or purifying selection ，即整個系統發育中的普遍的正選擇/純化選擇，而後者對應的是檢測 episodic (at a subset of branches) positive or purifying selection ，即侍畢檢測一部分branches的獨立正選擇/純化選擇。
①FEL 固定效應似然法(FEL, Fixed Effects Likelihood)
使用最大似然(ML)方法來推斷每個位點上的非同義(dN)和同義(dS)替換率，用於給定的編碼比對和相應的系統發育。該方法假設在整個系統發育過程中，每個位點的選擇壓力是恆定的。 注意，FEL適合小到中型數據
②SLAC (Single-Likelihood Ancestor Counting)
對於給定的編碼比對和相應的系統發育使用最大似然(ML)和計數方法的結合來推斷每個位點上的非同義(dN)和同義(dS)替換率。像FEL一樣，該方法假設在整個系統發育過程中，每個位點的選擇壓力是恆定的。 SLAC和FEL精準度相似，但適合更大的數據，並且不適合高度分歧的序列
③⭐FUBAR(Fast, Unconstrained Bayesian AppRoximation)⭐
使用貝葉斯方法來推斷給定編碼比對和相應系統發育的每個位點上的非同義(dN)和同義(dS)替換率。該方法假設在整個系統發育過程中，每個位點的選擇壓力是恆定的。 FUBAR適用於中到大數據集，預計在檢測位點的普遍選擇方面比FEL更有效。FUBAR是推斷pervasive selection的首選方法。
MEME(Mixed Effects Model of Evolution)
MEME(混合效應進化模型)採用混合效應最大似然方法來檢驗個別位點是否受到episodic positive或多樣化選擇的影響的假設。換句話說，MEME的目的是檢測在一定比例的分支下正選擇下進化的位點。
對於每個位點，MEME推測兩種ω值，以及在給定的分支下，以此ω進化的概率。為了推斷ω，MEME會推斷α(dS) 和兩個不同的β(dN),β−和β+。在空模型和備擇模型中，MEME強制β−≤α。因此β+是空模型和備擇模型不同的關鍵：在空模型中，β+被限制為≤α，但在備擇模型中不受限制。最終，當β+>α時，位點被推斷為正選擇，並使用似然比檢驗顯示顯著。
FADE(FUBAR Aproach to Directional Evolution)
是一種基於FUBAR引入的貝葉斯框架(Bayesian framework)的方法，用來測試蛋白質比對中的位點是否受定向選擇的影響。具體地說，FADE將系統地測試，對於比對中的每個位點，與背景分支相比，一組指定的前景分支是否顯示對特定氨基酸的替代偏向。該偏差參數的高值表明該位點對特定氨基酸的取代作用大大超過預期。使用貝葉斯因子(BF)評估FADE的統計顯著性，其中BF>=100提供了強有力的證據，表明該位點正在定向選擇下進化。
重要的是，與HyPhy中的大多數方法不同，FADE不使用可逆的馬爾可夫模型，因為它的目標是檢測定向選擇。因此， FADE分析需要一個有根的系統發育。在使用FADE進行分析之前，可以使用基於瀏覽器的交互工具「Phylotree.js」來幫助建立樹的根。
aBSREL (adaptive Branch-Site Random Effects Likelihood)
是常見的「Branch-Site」類模型的改進版本。aBSREL既允許分支的先驗指定（即指定foreground branches)來測試選擇，也可以以探索性的方式測試每個譜系以進行選擇( p-value將自動進行BH校正，為什麼叫探索性的方法呢，因為你可以先不指定foreground，來看看哪個支的pvalue更低，然後來針對那一支進行進一步的選擇壓力分析 ）。請注意，探索性的方法將犧牲功效。 ⭐aBSREL是在各個分支檢測正選擇的首選方法⭐，需要注意一點的是，aBSREL是多次獨立對指定的每一支進行檢驗的，也就是說，你指定了許多的branches，實質上和多次指定不同一個branch來多次運行，效果是一樣的，而並非將這些branches視為一個整體去做檢測
BUSTED(Branch-Site Unrestricted Statistical Test for Episodic Diversification)
通過測試一個基因是否在至少一個分支的至少一個位點上經歷了正選擇，BUSTED(分支位點無限制統計檢驗)提供了一個全基因(非位點特異性)正選擇的測試。當運行BUSTED時，用戶可以指定一組前景支來測試正選擇(其餘分支被指定為「背景」)，或者用戶可以測試整個系統發生的正選擇。在後一種情況下，整個樹被有效地視為前景，正選擇的檢驗考慮整個系統發育。這種方法對於相對較小的數據集(少於10個分類單元)特別有用，在這些數據集中，其他方法可能沒有足夠的功效來檢測選擇。 這種方法不適用於確定有正選擇的特定位點。
對於每個系統發育分區（前景和背景分支位點），BUSTED擬合了一個具有三個速率類的密碼子模型，約束為ω1≤ω2≤1≤ω3。與其他方法一樣，BUSTED同時估計每個分區屬於每個ω類的位點的比例。這種模型作為選擇檢驗中的替代模型，被稱為無約束模型。然後，BUSTED通過比較這個模型與前景分支上ω3=1（即不允許正選擇）的空模型的擬合度來測試正選擇。這個零模型也被稱為約束模型。如果零假設被拒絕，那麼就有證據表明，至少有一個位點在前景枝上至少有一部分時間經歷了正選擇。重要的是，一個顯著的結果並不意味著該基因是在整個前景的正選擇下進化的。
RELAX
RELAX是一種假設檢驗框架，它檢測自然選擇的強度是否沿著一組指定的測試分支被放鬆或加強。因此，RELAX不是明確檢測正選擇的合適方法。相反，RELAX在識別特定基因上自然選擇嚴格程度的趨勢和/或變化方面最有用。K>1表示選擇強度加強，K<1表示選擇壓力放鬆。
RELAX需要一組指定的 "測試 "分支與第二組 "參考 "分支進行比較（注意，不必分配所有的分支，但測試集和參考集各需要一個分支）。RELAX首先對整個系統發育過程擬合一個具有三個ω類的密碼子模型（空模型）。然後，RELAX通過引入作為選擇強度參數的參數k（其中k≥0）作為推斷ω值的指數：ωk來測試放鬆/強化選擇。具體來說，RELAX固定推斷的ω值（都是ωk<1,2,3>），並對測試分支推斷出一個將比率修改為ωk<1,2,3>的k值（替代模型）。然後，RELAX進行似然比檢驗，比較替代模型和空模型。
用法來說，以我用的2.2.0版本為例子，（2.4.0直接用hyphy命令即可）
模型上來說：
如果你要檢測類似paml中的M8位點模型，最好用FUBAR，如果是小數據，則用FEL，大數據並且分歧度不是很高用SLAC。
如果你要檢測某個前景支當中正選擇位點，最好用MEME。
如果你要檢測單獨的某個branch是否存在正選擇，最好用aBSREL。
如果你要檢測一系列的branches的正選擇，即檢驗你的這個基因，在指定的branches的任意一個位點是否在某段時間經歷過正選擇，則用BUSTED，BUSTED是不適合檢測單獨位點的正選擇的。
如果你要檢測選擇壓力的放鬆/加強，用RELAX。
如果你要用蛋白序列來檢測氨基酸位點正選擇/定向選擇，用FADE。
最後再提一下，幾乎所有模型(還有一些沒常用的模型沒有提到)都可以分為指定前景和不指定前景的模式運行，但不是都適合，就像官方說的那樣，根據你的目的不同，會有最優選擇，當然你也可以把某種模型都跑一遍，比如各種位點模型都走個流程，並且你也可以結合paml的模型，例如，對於檢測Pervasive selection的位點模型，你可以結合paml的M8、M2a來分析。對於檢測episodic selection的branch-site，你可以結合paml的branch-site modelA和BUSTED/aBSREL來比較分析。
以上的所有文字，都是筆者根據官方以及一些文獻當中對於hyphy的使用總結、翻譯，如有錯誤使用之處，還請各位多多指正。