使用PAML软件利用核酸序列进行系统发育参数分析

PAML软件中的baseml命令可以利用Maxmum Likelihood方法计算系统发育中的相关参数。basml命令能以分支树拓扑结构信息和多序列比对信息作为输入,在参数配置文件中设置一些控制程序运行的参数,来对一些进化过程中的参数进行分析。这些参数主要有:

1. 序列中各碱基频率。
2. 两两物种之间的距离。
3. 系统发育树的枝长信息。
4. 碱基替换矩阵。

PAML软件的baseml示例数据位于软件安装目录下。

1. 输入文件(1/3):含有拓扑结构的树文件 input.trees

输入一个树文件input.trees。该文件可以是有根树,也可以是无根树,也可以带有枝长信息。例如:

5  1
(((Human:0.1, Chimpanzee:0.2):0.8, Gorilla:0.3):0.7, Orangutan:0.4, Gibbon:0.5);

文件内容分两行:第一行表述树中有5个物种,共计1个树,两个数值之间用任意长度的空分割。此外,Newick格式的树尾部一定要有分号,没有的话程序可能不能正常运行。

2. 输入文件(2/3):核酸序列比对结果 input.txt

5   895

Human AAGCTTCACCGGCGCAGTCATTCTCATAATCGCCCACGGACTT...
Chimpanzee AAGCTTCACCGGCGCAATTATCCTCATAATCGCCCACGGACTT...
Gorilla AAGCTTCACCGGCGCAGTTGTTCTTATAATTGCCCACGGACTT...
Orangutan AAGCTTCACCGGCGCAACCACCCTCATGATTGCCCATGGACTC...
Gibbon AAGCTTTACAGGTGCAACCGTCCTCATAATCGCCCACGGACTA...

PAML要求输入的Phylip格式,其物种名和后面的序列之间至少间隔两个空格(是为了允许物种名的属名和种名之间有一个空格)。

3. 输入文件(3/3):baseml命令的配置文件baseml.ctl

软件难点在于配置文件的理解,常用示例如下:

*输入输出参数:
      seqfile = input.txt   *输入核酸多序列比对结果文件路径
     treefile = input.trees *输入系统发育树文件路径
      outfile = mlb         *设置输出文件路径
        noisy = 2           *设置屏幕上输出信息模式,值越大则输出越多信息:0,1,2,3。
      verbose = 0           *结果文件中的输出信息方式:1,输出详细信息,包含碱基序列;0,精简模式结果(推荐)。
        getSE = 0           *设置是否计算并获得各参数的标准误:0,不需要;1,需要。
 RateAncestor = 0           *设置是否计算序列中每个位点的碱基替换率:0,不需要;1,需要。当设置成1时,会在结果文件rates中给出各个位点的碱基替换率;同时也会进行祖先序列的重构,在结果文件rst中有体现。

*数据使用说明参数:
      runmode = 0     *设置程序获取进化树拓扑结构的方法:0,直接从输入文件中获取进化树拓扑结构(推荐);1,程序以输入文件中的多分叉树作为起始树,并使用星状分解法搜索最佳树;2,程序直接以星状树作为起始树,并使用星状分解法搜索最佳树;3,使用逐步添加法搜索最佳树;4,使用简约法获取起始树,再进行邻近分枝交换寻找最优树;5,以输入文件中的树作为起始树,再进行邻近分枝交换寻找最优树。 
* fix_blength = 1     *设置程序处理输入树文件中枝长数据的方法:0,忽略输入树文件中的枝长信息;-1,使用一个随机起始进行进行计算;1,以输入的枝长信息作为初始值进行ML迭代分析,此时需要注意输入的枝长信息是碱基替换率,而不是碱基替换数;2,不需要使用ML迭代计算枝长,直接使用输入的枝长信息,需要注意,若枝长信息和序列信息不吻合可能导致程序崩溃;3,让ML计算出的枝长和输入的枝长呈正比。
    cleandata = 1     *设置是否移除不明确字符(N、?、W、R和Y等)或含以后gap的列后再进行数据分析:0,不需要,但在序列两两比较的时候,还是会去除后进行比较;1,需要。
        clock = 0     *设置进化树中各分支上的变异速率是否一致,服从分子钟理论:0,变异速率不一致,不服从分子钟理论(当输入数据中物种相差较远时,各分支变异速率不一致);1,所有分枝具有相同的变异速率,全局上服从分子钟理论;2,进化树局部符和分子钟理论,程序认为除了指定的分支具有不同的进化速率,其它分枝都具有相同的进化速率,这要求输入的进化树信息中使用#来指定分支;3,对多基因数据进行联合分析。
        Mgene = 0     *设置是否有多个基因的多序列比对信息输入,以及多各基因之间的参数是否一致:0,输入的多序列比对文件中仅包含一个基因时,或多个基因具有相同的Kappa和Pi参数;1,输入文件包含多个基因,这些基因之间是相互独立的(这些基因之间具有不同的Kappa和Pi值,且其进化树的枝长也不相关);2,输入文件包含多个基因,这些基因具有相同的Kappa值,不同的Pi值;3,输入文件包含多个基因,这些基因具有相同的Pi值,不同的Kappa值;4,输入文件包含多个基因,这些基因具有不同的Kappa值和不同的Pi值。当值是2、3或4时,多个基因的进化树枝长虽然长度不一样,但是呈正比关系。这点和参数值等于1是不一样的。
*       ndata = 1     *设置输入的多序列比对文件中的数据个数。
*       icode = 0     *设置遗传密码。当设置参数RateAncestor = 1时生效,有利于根据密码子重构编码蛋白的祖先DNA序列。
   Small_Diff = 7e-6  *设置一个很小的值,一般位于1e-8到1e-5之间。推荐检测设置不同的值在比较其结果,需要该参数值对结果没什么影响。
	  
*位点替换模型参数:
        nhomo = 1     *设置碱基替换模型中的频率参数计算方法:0,使用输入数据中的观测到的所有序列的平均碱基频率;1,使用ML迭代方法计算碱基频率Pi(T)、Pi(C)、Pi(A)和Pi(G),适合model参数值为F81、F84、HKY85、T92、T93或REV模型;2,进化树所有的分枝具有相同的Kappa值,适合K80、F84和HKY85这三个含有Kappa参数的模型;3,外部节点所在的分枝使用一组碱基频率,所有内部节点所在分枝使用一组碱基频率,root节点所在分枝使用一组碱基频率,也称为N1模型,需要输入有根树,适合model参数值为F84、HKY85和T95模型;4,root节点所在分枝使用一组碱基频率,所有其它分枝使用一组碱基频率,也称为N2模型,需要输入有根树,适合model参数值为F84、HKY85、T95和GTR/REV模型;5,让用户自己定义各组分枝使用的碱基频率,需要在进化树使用#对各分枝进行编号。
        model = 7     *选取碱基替换模型方法:0,JC69;1,K80;2,F81;3,F84;4,HKY85;5,T92;6,TN93;7,REV;8,UNREST;9,REVu;10,UNRESTu。
                      *一般为了得到较为准确的枝长信息,至少选择4,使用HKY85碱基替换模型,程序会使用ML方法计算Kappa值(转换率/颠换率比值)。
                      *推荐选择7,也称为GTR,程序会使用ML计算碱基替换模型。
    fix_alpha = 0     *序列中不同的位点具有不同的碱基替换率,服从discrete-gamma分布,该模型通过alpha(shape参数)和ncatG参数控制。该参数设置是否给定一个alpha值:0,使用ML方法对alpha值进行计算;1,使用下一个参数设置一个固定的alpha值。
        alpha = 0.5   *设置一个固定的alpha值或初始alpha值(推荐设置为0.5)。该值小于1,表示只有少数热点位置的碱基替换率较高;该值越小,表示碱基替换率在各位点上越不均匀;若设置fix_alpha = 1且alpha = 0则表示所有位点的碱基替换率时恒定一致的。
       Malpha = 0     *当输入数据时多基因时,设置这些基因之间的alpha是否不一致:1,不同基因之间具有不同的alpha值;0,所有基因都具有相同的alpha值。
        ncatG = 5     *序列中不同位点的变异速率服从GAMMA分布的,ncatG是其一个参数,一般设置为5,4,8或10,且序列条数越多,该值设置越大。
        nparK = 0     *设置rate-class models:1,不同位点的变异速率是独立的,限定ncatG各类的概率都相等;2,不同位点的变异速率是独立的,不限定ncatG各类的概率都相等;3,邻近位点的变异速率是Markov相关的,限定ncatG各类的概率都相等;4,邻近位点的变异速率是Markov相关的,不限定ncatG各类的概率都相等。此外,设置该参数值不为0,比较消耗计算并且要求nhomo = 0,否则程序报错。
    fix_kappa = 0     *设置是否给定一个Kappa值(该参数在model设置为K80、F84或HKY85时有效):0,通过ML迭代来估算Kappa值;1,使用下一个参数设置一个固定的Kappa值。对于JC69和F81,该参数没有效果;对于TN93和REV,则分别对应2个和5个rate参数。
        kappa = 5     *设置一个固定的Kappa值,或一个初始的Kappa值。
       method = 0     *设置评估枝长的ML迭代算法:0,使用PAML的老算法同时计算所有枝长,在clock = 0下有效;1,PAML新加入的算法,一次对一个枝长进行计算,该算法仅在clock参数值为1,2或3下工作。

4. 运行baseml命令对核酸序列进行进化参数计算

程序运行的命令:

baseml baseml.ctl

5. 结果文件解释

程序的主要结果文件是mlb,其主要内容有:

1. 四种碱基观测到的频率及其平均值。
2. 各序列之间的距离矩阵。
3. 各分枝的长度、总的枝长、含有枝长信息的系统发育树。
4. ML计算出来的碱基频率。
5. 碱基替换矩阵。
6. Gamma分布的alpha值。
7. 若使用K81、F84和KEY85模型,则给出Kappa参数值;若选择TN93或REV模型,则分别给出2个和5个rate参数值。

6. baseml命令的意义

PAML软件不擅长搜索进化树的最优拓扑结果,它仅适合于对小数据(序列条数少于10)进行最优树的搜索。一般使用其它ML软件进行系统发育树分析,就可以得到包含枝长信息的最优ML树。再将该树信息输入给basml命令进行分析,又得到一个包含枝长信息的进化树,此外还得到一些没有什么意义的参数值。很多人会疑惑basml命令又有什么用呢?

在我看来,baseml命令的用处:(1) 通过这个软件的使用和练习,可以更加深入了解进化中各参数的意义和系统发育树枝长计算的原理。(2) 再进行物种分歧时间计算时,是需要先通过baseml命令计算相关参数和枝长的。使用mcmctree命令进行物种分歧时间计算时,第一步会调用baseml进行计算。(3) 若使用数百上千个基因构建物种树,可以使用baseml命令对各个基因分别计算进化树枝长,然后通过比较枝长来去除掉和物种树枝长差异较大的基因,再进一步得到更准确的物种树。