‘壹’ 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的使用总结、翻译，如有错误使用之处，还请各位多多指正。