系统发育摘PAMl学习笔记蜗牛.docx

系统发育摘PAMl学习笔记蜗牛.docx

《系统发育摘PAMl学习笔记蜗牛.docx》由会员分享，可在线阅读，更多相关《系统发育摘PAMl学习笔记蜗牛.docx（15页珍藏版）》请在冰豆网上搜索。

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

PAML与Branchmodel

——以灵长目动物溶菌酶编码基因适应性进化分析为例解读BranchModel

1。

什么是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中存在与不存在ω1〈1这一约束的两种情况下所对应的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（http：

//linxiao.name/archives/172）网站中所述方法有异（将df值与2△l值输入程序中，回车查看显著性水平）,但网站中是在Linux平台下操作,而本法是在windows平台下操作。

另外，在Chi程序的windows版本中并未发现任何输入的光标。

5.Prof.Yang对灵长目动物溶菌酶不同支系的适应性进化分析

5。

1数据和方法（Dataandmethods）

5.1.1数据（Data）

本文所涉及的数据主要分为两部分内容，首先是大数据集（largedataset）,这主要包括Prof.Messier分析的24条灵长目动物溶菌酶编码基因中有显著不同的19条序列（Distinctsequence），其系统发育树如图2所示。

新大陆猴

叶猴

6.langurSen&Sve

7.langurTob&Tfr

8.DouclangurPne

9.probiscisNla

5.colobusCgu&Can

10.baboonPcy

11.mangabeyCat

12.rhesusMmu

13.AllenAni

14.talapoinMta

15.patasEpa

16.vervetCae

1.human

2.chimpbonobogorilla

3.orangutanPpy

4.gibbonGgo

17.squirrelm

18.tamarinSoe

19.MarmosetCja

0.02

h

c

颊囊猴

人科动物

图2Prof。

Yang所分析的灵长目溶菌酶系统发育树

Fig2。

PhylogenyoftheprimatelysozymeanalyzedbyProf。

Yang

在图2中,Branchh和Branchc是本文分析是所选取的前景枝（foregroundclade）。

文中这两个前景枝的选取是根据Prof。

Messier于1997年的研究结果.

小数据集包括7条序列，其来源是从图2中四个分支中各挑选出几条具有代表性的序列,重新进行分析,其系统发育树如图3所示。

图3从图2挑选出来的四个分支代表序列的系统发育分析

Fig3PhylogenyofasubsetofsevenprimatelysozymeselectedfromthoseofFig1.torepresentthefourmajorgroupsofspecies.

Prof.Yang认为，对于这种大数据集和代表性的小数据集的分析比较能够之时取样方法对于所得结果的一个敏感度。

（原文:

Differencesbetweentheanalysesofthetwodatasetswillgiveusanindicationofthesensitivelyoftheresultstospeciessampling。

）

5。

1。

2方法（Methods）

第一，对小数据集controlfile文件的解读

图4溶菌酶小数据集的controlfile

Fig4Thecontrolfileoflysozymesmalldataset

在图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病毒适应性进化分析为例解读SitesModel

1.什么是SitesModel？

SitesModel是PAML软件CODEML程序的一个正选择作用分析模型，其主要观点是同一序列不同位点的omega值不同。

在进行sitesModel分析时，需要设置controlfile中的Model=0,Nssites命令在此是一个变量，根据不同Model的选择设置不同的值。

值得注意的是，以此可以选择多个sitesModel。

如Nssites=01378.

2。

不同的sitesModel表示什么意思？

M0即one-ratioModel,值得是所有位点的omega值是恒定的；

M1表示加假定有一部分位点的omega值为0，其他位点的omega值为1；

M2是在M1的基础上增加了第三类omega值，该类omega是通过数据计算得到的，有可能大于1;

M3假定所有位点的omega值呈简单的离散分布趋势；

M5假定所有位点的omega值呈简单的gamma分布趋势；

M7假定所有位点的omega属于矩阵（0,1）且呈beta分布；

M8是在M7的基础上增加另一类omega值，该值可通过计算得到，可以大于1.

M8a与M8类似，只不过新增加的omega值等于1。

3。

不同Model的比较可以得到什么样的结果？

首先是M0与M3的比较，该比较与BranchModel中的Model之间的比较是一样的，首先计算2△l值,然后在一定df值下进行显著性水平计算。

这里需要注意的是,在参阅了Prof.Yang关于PAML的一些参考材料之后,我们发现sitesModel比较的df值一般取2.

在sitesModel中,M0表示oneratioforallsites，M3表示所有位点的omega值呈简单的离散分布。

对于这两个模型的比较并非用于正选择作用的检测,而是用于位点间omega值是否一致的检测.

M1andM2以及M7andM8是用于正选择作用的检测，但Prof.Yang认为，TheM1—M2comparison与theM7-M8comparison相比，更加的稳定。

（原文：

TheM1—M2comparisonappearstobemorerobust（orlesspowerful）thantheM7-M8comparsion）

此外,还有一类比较是M8toM8a,其中M8和M8a是两个极为类似的Model.在涉及到positiveselection的文章中，对于这两个model的比较并不常见，而且说明书中也并未给出明确的比对结果意义。

4。

如何检测positivesites?

在CODEML中，positivesites的检测流程主要如图1所示.

PPvaluecomputation

Likelihoodratiotest

CODEMLcomputation

图1positivesites的检测流程

Fig1Theprocessofpositivesites

其中CODEMLcomputation主要是对controlfile中的命令值进行设定之后，运行CODEML程序，并在resultfile中查看运算结果。

Likelihoodratiotest如question3所示,即对两个模型进行显著性水平比较。

PPvaluecomputation主要是指位点后验概率的计算，该结果是显示在mainresultfile—mlc文件中。

CODEML程序中常见的计算后验概率的方法有BEB和NEB.与BEB相比，NEB在计算的过程中往往会忽略抽样误差.因此，Prof。

Yang建议在读取运算结果时，可以直接将NEBresult忽略，但值得注意的是，BEB只能在M2a和M8model下运行.

5。

以example中controlfile文件为参考，解读sitemodel下的controlfile.

图2sitemodel下的controlfile截图

Fig2Thefigofcontrolfileundersitemodel

Sitemodel计算的controlfile与Branchmodel中大致相似，但在sitemodel中应当注意，model=0是一个恒定值,Nssites命令可以设置不同的模型参数。

PAML与Branchsitesmodel

1.什么是Branch—sitemodel？

Branch—sitemodel其实是Branchmodelandsitemodel的合集，在该model下，不仅假定位点间的omega值是变化的，同时也假定支系间的omega值是变化的。

该model主要用于检测前景枝中正选择作用对部分位点的影响.

2。

Branch-sitemodel中都有哪些模型？

ModelA=（model=2NSsites=2）；ModelB=（Model=2Nssites=3）;ModelC=（Model=3Nssites=2）；ModelD=（Model=3Nssites=3）.注意ModelD需要ncatG参数设置位点的类型（以omega值进行分类）,可以使用的ncatG值有2或者3，在ModelA，ModelBandModelC中，ncatG的取值是被自动忽略的。

3.在Branch—sitemodel中如何进行模型比较？

ModelA与ω2=1的nullModel进行比较，该nullModel可以通过设定fix_omega=1andomega=1进行确定。

ThecomparisonofModelAandnullModel对应的df值为1；（注：

在nullModel中omega值的设定是ω2=1，这便限定了在nullModel中，background中的sites均处于negativeselection，foregroundBranch中的sites均处于neutralselection。

而在alternativeModel中,foregroundBranch中的site对应的omega值是大于1的。

因此，alternativeModelandnullModel之间的显著性水平检测可以直接用来检测foregroundBranch中的positiveselection位点）

在Branch-siteModel中，ModelA还可以与siteModel中的M1a进行显著性水平比较。

若这两者之间存在显著性差异，则可以说明要么foregroundBranch受到的选择约束较为宽松,要么foregroundBranch受到明显的正选择作用.

ModelB通常与siteModel中的M3进行显著性比较，对应的df值为2

ModelC与siteModel中的M1a进行比较，对应的df值为3.

ModelD通常与siteModel中的M3进行比较，若树文件为无根树，则df=1，若树文件为有根树，则df=2.

就Model的提出理念而言，ModelAandModelB侧重于寻找在进化过程中受正选择作用的点，而ModelD则与之不同，其不再局限于正选择作用，而是涉及到多种选择作用（divergentselectivepressures）。

4。

在Branch-siteModel中的其他注意事项

在ModelCandModelD中，不同omega值的Branch类型不局限于两种，使用者可以自行设置多种Branchtypes，最多可以设置10-12种。

另外,在ModelC中，其ω1是一个定值,而在ModelD中，ω1则是一个自由变化的参数。

Maximumlikelihoodmethodsfordetectingadaptiveevolutionaftergeneduplication

导读：

基因组数据的不断报道使得基于大数据分析的基因家族进化成为可能，在此基础上，本文提出了一种基于Maximumlikelihoodmethods的方法，对基因家族形成过程中的适应性进化进行检测，并以灵长目的ECP—EDN基因家族为例对此方法进行了说明。

1.选择压力支系特异性检验（Detectinglineage-specificchangesinselectivepressure）

若基因家族的分化是通过正选择作用推动的，则在复制事件之后,会立刻出现非同义替换速率大于同义替换速率的现象.但是,若这一基因家族在其余的时间都是受到purifyingselection的作用,那么仅仅是两个序列之间的比对是很难发现omega值大于1的位点.

鉴于这一问题，本文提出了用于检测基因复制后适应性进化的最大似然法，旨在完成以下目标：

（1）同一个基因进化历史中不时间点的选择压力；

（2）这些选择压力是否不同。

1。

1支系间选择压力可变模型（Modelsofvariableselectivepressuresamongbranches）

1.1.1.模型介绍

（1）one—ratiomodel,即Phylogenetictree中所有sites的omega值是一个恒定值。

ω0=ω1=ω2=ω3；

（2）Two—ratiomodel,复制事件之前的omega值和复制事件之后的omega值不相等.ω0≠ω1=ω2=ω3;

90

图1假设基因家族的进化树

Fig1Phylogenyofahypotheticalgenefamily

（3）Three-ratiomodel，该模型共假设三个omega值,即复制事件之前的omega值,复制事件之后的omega值，以及后期分支中的omega值。

ω0≠ω1≠ω2=ω3。

（4）Four-ratiomodel，该模型共假设四个omega，不仅复制事件前后的omega值不同以外,随后的分支上的两个omega值也不相同。

即ω0≠ω1≠ω2≠ω3。

1。

1.2模型比较可以揭示的问题（likelihoodratiotest）

（1）oneratiomodelwithtworatiomodel：

揭示基因家族的复制前后所受的选择压力是否相同；

（2）tworatiomodelwiththreeratiomodel:

揭示基因家族复制后到分化前以及分化后这段时间的选择压力是否发生变化。

（3）threeratiomodelwith