【系统发育】摘PAMl学习笔记-蜗牛

PAML与Branchmodel——以灵长目动物溶菌酶编码基因适应性进化分析为例解读BranchModel1.什么是Branchmodel？Branchmodel是PAML软件CODEML程序中通过likelihoodratiotest（LRT）进行不同支系间（lineages）适应性进化检测的一种模型。该模型通过限制（constraint）系统发育树中不同分支上的omega(dN/dS)值的异同，并对不同的限制进行显著性分析（PAML软件中的Chi2程序），进而得到较为可靠地分析结果。在该法提出之前，不少学者通过简约法（parsimonymethod）或者似然法（likelihoodmethod）先重建祖先序列（ancestralsequence），然后通过对构建的祖先序列的omega值估算进而预测不同支系的适应性进化特征。诸如Prof.Messier等对于灵长目动物溶菌酶的分析便是如此。Prof.Yang认为，从统计学的角度而言，这种将预测的数据当做真实观测数据的分析理念存在一定的随机误差（randomerrors）和系统误差(systematicerrors),本身并不是一种严谨的统计学方法。Prof.Yang所提出的Branchmodel巧妙地避开了直接利用ancestralsequence进行支系间适应性进化检测的流程，而是通过平均统计每一个节点（eachnode）中可能的ancestralsequence，根据其相对发生似然率（relativelikelihoodsofoccurance）进行加权分析。此外，Branchmodel还考虑到了（takeintoaccount）密码子转移/颠换速率偏差（transition/transversionratebias）和非均匀密码子（nonuniformcondonusage）这些与omega值计算有着显著关系的影响因素。2.Branchmodel中存在哪些假设模型，在CODEML程序的controlfile文本中如何选择？Branchmodel主要是对系统发育树中的不同支系的omega值的异质性进行界定，主要的model有：one-ratiomodel，即系统发育树中所有支系的omega值是相等的；free-ratiomodel，该模型指的是系统发育树中所有支系的omega值是不相等的。这两个假设是不同支系omega取值的两个极限。此外，还可以设定前景枝（foregroundclade）,假定其与其余支系（又称背景枝backgroundclade）的omega值不同。前景枝可以根据需要设置多个。在controlfile中，Model=0表示one-ratiomodel，Model=1表示free-ratiomodel.Model=2表示系统发育树中不同omega值得个数，其中所选择前景枝的个数为(n-1)。值得注意的是，当设置Model=2,3，……，n时，需要在treefile中标记所要设置的前景枝，可以标记一个，也可以标记多个。树标记格式如下所示：((1,2),3)#1,4,5);该treefile表示Clade1，2and3为前景枝，其对应的omega值为ω1（用#1表示），其余Clade为背景枝，对应的omega值为ω0（用#0指定，但在PAMl软件中，#0为默认值，故不需要在树中注明）。在resultfilemlc文件中，我们可以得到两个不同的omega值。3.通过BranchModel可以得到什么样的结论？3.1不同支系间的omega值是否显著不同这主要通过比较one-ratioModel和free-ratioModel对应的likelihoodvalues的差异进行说明。3.2前景枝和背景枝的omega是否显著不同这主要通过比较one-ratioModel和two-ratioModel的likelihoodvalues进行说明。3.3前景枝的omega值是否显著大于1（greaterthanone）这主要通过比较two-ratioModel中存在与不存在ω11这一约束的两种情况下所对应的likelihoodvalues进行说明。4.如何对不同Model的差异性进行比较？该比较主要在PAML软件Chi2程序中进行，首先在mlc文件中查找不同Model所对应的likelihoodvalues并计算不同Modellikelihood差值绝对值的两倍，即2△l=|l1-l2|。打开Chi2程序，界面如图1所示。图1Chi2程序主界面Fig1ThemainofChi2program通过上图查询df=10时，2△l值所对应的显著性水平，小于等于0.05时，被认为是存在显著性差异的,如图1中绿色框所标注。注：该法与linxiao.name（）网站中所述方法有异（将df值与2△l值输入程序中，回车查看显著性水平），但网站中是在Linux平台下操作，而本法是在windows平台下操作。另外，在Chi程序的windows版本中并未发现任何输入的光标。5.Prof.Yang对灵长目动物溶菌酶不同支系的适应性进化分析5.1数据和方法（Dataandmethods）5.1.1数据（Data）本文所涉及的数据主要分为两部分内容，首先是大数据集（largedataset），这主要包括Prof.Messier分析的24条灵长目动物溶菌酶编码基因中有显著不同的19条序列（Distinctsequence），其系统发育树如图2所示。图2Prof.Yang所分析的灵长目溶菌酶系统发育树Fig2.PhylogenyoftheprimatelysozymeanalyzedbyProf.Yang6.langurSen&Sve7.langurTob&Tfr8.DouclangurPne9.probiscisNla5.colobusCgu&Can10.baboonPcy11.mangabeyCat12.rhesusMmu13.AllenAni14.talapoinMta15.patasEpa16.vervetCae1.human2.chimpbonobogorilla3.orangutanPpy4.gibbonGgo17.squirrelm18.tamarinSoe19.MarmosetCja0.02hc颊囊猴人科动物叶猴新大陆猴在图2中，Branchh和Branchc是本文分析是所选取的前景枝（foregroundclade）。文中这两个前景枝的选取是根据Prof.Messier于1997年的研究结果。小数据集包括7条序列，其来源是从图2中四个分支中各挑选出几条具有代表性的序列，重新进行分析，其系统发育树如图3所示。图3从图2挑选出来的四个分支代表序列的系统发育分析Fig3PhylogenyofasubsetofsevenprimatelysozymeselectedfromthoseofFig1.torepresentthefourmajorgroupsofspecies.Prof.Yang认为，对于这种大数据集和代表性的小数据集的分析比较能够之时取样方法对于所得结果的一个敏感度。（原文：Differencesbetweentheanalysesofthetwodatasetswillgiveusanindicationofthesensitivelyoftheresultstospeciessampling.）5.1.2方法（Methods）第一，对小数据集controlfile文件的解读图4溶菌酶小数据集的controlfileFig4Thecontrolfileoflysozymesmalldataset在图4中，seqfile，treefile和outfile分别指代的是分析所用的序列文件，树文件以及主要结果输出文件。Noisy命令控制屏幕中输出文件的数量，一般取值为9；Verbose命令控制文件中输出结果的数量；Runmode命令控制是否使用树文件，常见的值有0，1，2，3，-2等，常见取值为0，表示从树文件中获取树的拓扑结构，-2是pairwisealignment的命令，若选择runmode=-2，则不需要输入树文件。Seqtype命令是指代输入系列文件的类型，可取值是1，2，3，seqtype=1表示序列文件为condonssequence，seqtype=2表示序列文件为aminoacidsequence，seqtype=3表示序列文件为通过condonssequence翻译得到的aminoacidsequence；CodonFreq命令是指定密码子使用频率模型的类型，其所取数值0,1,2,3分别代表不同的密码子使用频率模型。一般取值为CondonFreq=2,有人对不同CondonFreq的取值进行了比较，发现CondonFreq=1和CondonFreq=2并未有显著性的差别。Clock命令用来指代不同支系间的速率（指的是什么速率？）是恒定的还是变化的（原文：clockspecifiesmodelsconcerningrateconstancyorvariationamonglineages）。可取的值有0，1，2和3，其中0表示支系与支系之间的速率是不断变化的；1表示所有支系的速率相同；2表示系统发育树中大多树支系的速率是相同的，只有少部分特意指定的支系速率不同；3用于对多个基因或者多分区数据进行联合分析时所选取的参数。Runmode=1or2时，treefile必须是rootedtree。在常见的分析中，clock=0的命令居多。Model命令主要是对系统发育树中不同支系的omega值进行限定。可取值有0,1,2。0表示系统发育树所有支系均有相同的omega值，1表示omega值是自由的，不同支系间各不相同。2表示指定的前景枝和背景枝之间的omega值不同，一般而言，有n个前景枝，便有n+1个omega值。对于BranchModel，一般选择Model=2，但在该模式下，需要对树文件进行前景支的定义。另外，对于小数据集的BranchModel分析还可以选择Model=1，但由于Model=1所涉及的参数数量较多，不建议对大数据集进行Model=1的分析。NSsites是在SiteModel和Branch-siteModel中进行设置的参数，在BranchModel中一般设置为0。icode用来指定遗传密码子的类型,可取值有0——10共11个数字，0表示通用密码子，也是最为常见的一个设置。该命令主要根据所分析的序列来源进行具体的设置。fix_kappa命令用来指代K80，F84和HKY85中的kappa值是一个固定值还是根据数据进行迭代分析得到的。Fix_kappa=1表示kappa值是一个固定值。Fix_kappa=0表示kappa值是通过数据迭代运算得到的。而kappa是定义一个kappa的初始值，这个需要使用者自行设置。Fix_omegaandomega，fix-alphaandalpha以及fix-rhoandrho等的设置参照fix_kappa进行。注：在本文中，对于图5所涉及的一些模型，需要改变fix_omegaandomega值进行运算。图5本文中的一些模型Fig5somemodelsofthispaper若要限定omega值等于1，则需要在controlfile中设置fix_omega=1andomega=1.getSE命令是来设置是否需要计算估算参数的标准误（stranderrors）,若需要则getSE=1，不需要则getSE=0.Rateancestor=0or1,该值常设为0。Method命令主要控制在noclock模型下计算枝长的迭代算法。Method=0表示该迭代算法为老式算法，Method=1表示该迭代算法为新式算法，注意当clock=1,2,3时，Method只能取0.PAML与sitesmodel——以人类HIV病毒适应性进化分析为例解读SitesModel1