利用单拷贝直系同源基因蛋白序列做分歧时间树mcmctree

利用单拷贝直系同源基因蛋白序列做分歧时间树

if [ -f "./env.sh" ]; then

 source ./env.sh

else

 echo "env.sh file not in current directory, please check"

 exit

fi



mkdir 04.time_tree_pep

cd 04.time_tree_pep

#选择一个满意的树

ln -s ../03.Phylo_Tree/raxml.cds.raxml.support species.tre


#准备数据


ln -s ../03.Phylo_Tree/single_copy.cds_msa.4d.fasta .

ln -s ../03.Phylo_Tree/single_copy.cds_msa.fasta .

ln -s ../03.Phylo_Tree/single_copy.pep_msa.fasta .


#格式转换

trimal  -in single_copy.cds_msa.4d.fasta -out input.phy -phylip_paml


#进化树去掉不必要枝长和boot值

cat species.tre |sed -r 's/:[0-9\.]+//g' |sed -r 's/\)[0-9\.]+/)/g'>species_formated.tre


#######手动编写进化树,

#添加化石时间点:timetree http://timetree.org/

#设置进化树的root : https://itol.embl.de/

#添加化石校准点时间信息(格式是时间范围’>0.23<0.26’或者时间点‘@0.245’),单位时百万年前100Ma;

#再在首行添加两个数字(物种数量和树的数量),空格隔开,可得到input.tre文件。


# Tip: 可以利用figtree 标记一下,save as 编辑进化树;

# Monocots versus Dicots (173–148 Mya)

# A. thaliana versus P. trichocarpa (109–97 Mya)

# A. thaliana versus V. vinifera (115–105 Mya)

# H. annuus L. versus S. lycopersicum (107–93 Mya).



echo "

15 1

 (((Helianthus_annuus,((Sesamum_indicum,Olea_europaea),Solanum_lycopersicum))'>0.93<1.07',(((Populus_trichocarpa,(((Acer_yangbiense,Acer_truncatum),Citrus_sinensis),(Arabidopsis_thaliana,Gossypium_raimondii)))'>0.97<1.09',(Juglans_regia,Glycine_max)),(Vitis_vinifera,Malania_oleifera))'<1.15>1.05'),Oryza_sativa)'<1.73>1.48';

" >input.tre


##保存新的精进化树:  input.tre


############################################################################

## step1 估算位点替换率

### 将mcmctree.ctl  配置文件中seqtype 设置为2; 将usedata 设置为3 运行mcmctree



echo "

          seed = -1        *设置随机数作为seed,-1代表使用系统当前时间作为随机数

       seqfile = input.phy *输入多序列比对文件

      treefile = input.tre *带校准点(化石时间)的有根树文件

       outfile = out.txt   *输出文件

      mcmcfile = mcmc.txt  *输出的mcmc信息文件,可用Tracer软件查看


         ndata = 1   * 输入的多序列比对的数据区域的数量;多个数据phy格式合并

       seqtype = 2   * 设置多序列比对数据类型;0:核酸数据;1:密码子比对数据;2:氨基酸数据;

       usedata = 3   * 0: no data; 1:seq like; 2:use in.BV; 3: out.BV

                     * 是否利用多序列比对数据;

                     * 0: no data,不使用,不会进行likelihood估算,会快速得到mcmc树,但分歧时间不可用;

                     * 1:seq like,使用多序列比对数据进行likelihood估算,正常进行mcmc; usedata=1时model无法选择;

                     * 2:use in.BV, 进行正常的approximation likelihood分析,不读取多序列比对数据,直接读取当前目录的in.BV文件,in.BV是由usedata = 3时生成的out.BV重命名得来;此外,由于程序BUG,当设置usedata = 2时,一定要在改行参数后加 *,否则程序报错 Error: file name empty..;

                     * 3:out.BV,程序利用多序列比对数据调用baseml/codeml命令对数据进行分析,生成out.BV文件。由于mcmctree调用baseml/codeml进行估算的参数设置可能不太好(特别时对蛋白序列进行估算时),推荐自己修改软件自动生成的baseml/codeml配置文件,然后再手动运行baseml/codeml命令,再整合其结果文件为out.BV文件。


         clock = 2       * 设置分子钟算法,1: global clock,表示所有分支进化速率一致; 2: independent rates,各分支的进化速率独立且进化速率的对数log(r)符合正态分布; 3,correlated rates方法,和方法2类似,但是log(r)的方差和时间t相关。

                         *       TipDate = 1 100  *当外部节点由取样时间时使用该参数进行设置,同时该参数也设置了时间单位。具体数据示例请见examples/TipData文件夹。

        RootAge = '<2'  * constraint on root age, used if no fossil for root.设置root节点的分歧时间,一般设置一个最大值。


         model = 4      * models for DNA:

                        *  0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85;*设置碱基替换模型;当设置usedata = 1时,model不能使用超过4的模型,所以usedata = 1时用model = 4;usedata不等于1时,用model = 7,即GTR模型;

                        * models for codons:

                        *  0:one 恒定速率模型, 1:b 中性模型, 2:2 or more dN/dS ratios for branches 选择模型。

                        * models for AAs or codon-translated AAs:

                        *  0:poisson, 1:proportional, 2:Empirical, 3:Empirical+F

                        *  6:FromCodon, 7:AAClasses, 8:REVaa_0, 9:REVaa(nr=189)

         alpha = 0.5    * alpha for gamma rates at sites;*核酸序列中不同位点,其进化速率不一致,其变异速率服从GAMMA分布。一般设置GAMMA分布的alpha值为0.5。若该参数值设置为0,则表示所有位点的进化速率一致。此外,当userdata = 2时,alpha、ncatG、alpha_gamma、kappa_gamma这4个GAMMA参数无效。因为userdata = 2时,不会利用到多序列比对的数据。

         ncatG = 5      * No. categories in discrete gamma;设置离散型GAMMA分布的categories值。


     cleandata = 0    * remove sites with ambiguity data (1:yes, 0:no)?

                      * 设置是否移除不明确字符(N、?、W、R和Y等)或含以后gap的列后再进行数据分析: 0,不需要,但在序列两两比较的时候,还是会去除后进行比较;


       BDparas = 1 1 0.1   * birth, death, sampling;*设置出生率、死亡率和取样比例。若输入有根树文件中的时间单位发生改变,则需要相应修改出生率和死亡率的值。例如,时间单位由100Myr变换为1Myr,则要设置成'.01 .01 0.1'。

   kappa_gamma = 6 2       * gamma prior for kappa;设置kappa(转换/颠换比率)的GAMMA分布参数。

   alpha_gamma = 1 1       * gamma prior for alpha;设置GAMMA形状参数alpha的GAMMA分布参数。只对usedata=1时起作用,其他2,3不起作用


   rgene_gamma = 2 20 1 0   * au bu a prior, gamma prior for rate for genes;物种替换速率可查文献;设置序列中所所有位点平均[碱基/密码子/氨基酸]替换率的Dirichlet-GAMMA分布参数:alpha=2、beta=20、初始平均替换率为每100million年(取决于输入有根树文件中的时间单位)1个替换。若时间单位由100Myr变换为1Myr,则要设置成'2 2000 1'。总体上的平均进化速率为:2 / 20 = 0.1 个替换 / 每100Myr,即每个位点每年的替换数为 1e-9。

  sigma2_gamma = 1 10 1    * gamma prior for sigma^2     (for clock=2 or 3);设置所有位点进化速率取对数后方差(sigma的平方)的Dirichlet-GAMMA分布参数:alpha=1、beta=10、初始方差值为1。当clock参数值为1时,表示使用全局的进化速率,各分枝的进化速率没有差异,即方差为0,该参数无效;当clock参数值为2时,若修改了时间单位,该参数值不需要改变;当clock参数值为3时,若修改了时间单位,该参数值需要改变。


      finetune = 1: .1 .1 .1 .1 .01 .1  * times, rates, mixing, paras, RateParas;冒号前的值设置是否自动进行finetune,一般设置成1,然程序自动进行优化分析;冒号后面设置各个参数的步进值:times, musigma2, rates, mixing, paras, FossilErr。由于有了自动设置,该参数不像以前版本那么重要了,可能以后会取消该参数。


         print = 1      *设置打印mcmc的取样信息:0,不打印mcmc结果;1,打印除了分支进化速率的其它信息(各内部节点分歧时间、平均进化速率、sigma2值);2,打印所有信息。

        burnin = 2000   *将前2000次迭代burnin后,再进行取样(即打印出该次迭代估算的结果信息,各内部节点分歧时间、平均进化速率、sigma2值和各分支进化速率等)。

      sampfreq = 10      *每10次迭代则取样一次

       nsample = 100000  *当取样次数达到该次数时,则取样结束,同时结束程序。


*** Note: Make your window wider (100 columns) when running this program.


" >mcmctree1.ctl


mcmctree mcmctree1.ctl


### 将生成的tmp0001.ctl 拷贝为codeml.ctl,添加 clock = 1;修改 getSE 为 0.




#cp tmp0001.ctl codeml.ctl



echo "

  seqfile = tmp0001.txt      * 设置输入的多序列比对文件路径。

  treefile = tmp0001.trees   * 设置输入的树文件路径。

  outfile = tmp0001.out      * 设置输出文件路径。

  noisy = 3                  * 设置输出到屏幕上的信息量等级:0,1,2,3,9。

  seqtype = 2               * 设置输入的多序列比对数据的类型:1,密码子数据;2,氨基酸数据;3,输入数据虽然为密码子序列,但先转换为氨基酸序列后再进行分析。

  model = 2                 * 若输入数据是密码子序列,该参数用于设置branch models,即进化树各分枝的omega值的分布:0,进化树上所有分枝的omega值一致;1,对每个分枝单独进行omega计算;2,设置多类omega值,根据树文件中对分枝的编号信息来确定类别,具有相同编号的分枝具有相同的omega值,没有编号的分枝具有相同的omega值,程序分别计算各编号和没有编号的omega值。

                            * 若输入数据是蛋白序列,或数据是密码子序列且seqtype值是3时,该参数用于设置氨基酸替换模型:0,Poisson;1,氨基酸替换率和氨基酸的观测频率成正比;2,从aaRatefile参数指定的文件路径中读取氨基酸替换率信息,这些信息是根据经验获得,并记录到后缀为.dat的配置文件中。这些经验模型(Empirical Models)文件位于PAML软件安装目录中的dat目录下,例如,Dayhoff(dayhoff.dat)、WAG(wag.dat)、LG(lg.dat)、mtMAN(mtman.dat)和mtREV24(mtREV24.dat)等;


  aaRatefile = wag.dat      * 当对蛋白数据进行分析,且model = 2时,该参数生效,用于设置氨基酸替换模型。


  fix_alpha = 0             * 序列中不同的位点具有不同的碱基替换率,服从discrete-gamma分布,该模型通过alpha(shape参数)和ncatG参数控制。该参数设置是否给定一个alpha值:0,使用ML方法对alpha值进行计算;1,使用下一个参数设置一个固定的alpha值。

                         * 对于密码子序列,当NSsites参数值不为0或model不为0时,推荐设置fix_alpha = 1且alpha = 0,即不设置alpha值,认为位点间的变异速率一致,否则程序报错。若设置了alpha值,则程序认为不同密码子位点的变异速率不均匀,且同时所有位点的omega值一致,当然各分枝的omega值也会一致,这时要求NSsites和model参数值都设置为0(这一般不是我们需要的分析,它不能进行正选择分析了)。

  alpha = 0.5              * 设置一个固定的alpha值或初始alpha值(推荐设置为0.5)。该值小于1,表示只有少数热点位置的替换率较高;该值越小,表示位点替换率在各位点上越不均匀;若设置fix_alpha = 1且alpha = 0则表示所有位点的替换率是恒定一致的。


  ncatG = 5                  * 序列中不同位点的变异速率服从GAMMA分布的,ncatG是其一个参数,一般设置为5,4,8或10,且序列条数越多,该值设置越大。

  Small_Diff = 0.1e-6

  getSE = 0                  * 设置是否计算并获得各参数的标准误:0,不需要;1,需要。

  method = 0

  clock = 1                 * 设置进化树中各分支上的变异速率是否一致,服从分子钟理论:0,变异速率不一致,不服从分子钟理论(当输入数据中物种相差较远时,各分支变异速率不一致);1,所有分枝具有相同的变异速率,全局上服从分子钟理论;2,进化树局部符和分子钟理论,程序认为除了指定的分支具有不同的进化速率,其它分枝都具有相同的进化速率,这要求输入的进化树信息中使用#来指定分支;3,对多基因数据进行联合分析。

" >codeml.ctl


cp /share/work/biosoft/PAML/latest/dat/wag.dat .





### 将生成的tmp0001.trees文件进行修改,定根(添加一对括号即可),添加化石时间(单个时间点)


### 运行codeml, 查看tmp0001.out结果中的替换速率,

### 替换速率在“Substitution rate is per time unit” 下一行

codeml  codeml.ctl


## step2 使用近似似然法计算分化时间

###调整mcmctree.ctl 中的rgene_gamma 第二个参数b, 使得a/b 约等于之前得到的替换率, 将usedata 设置为3 运行




echo "

          seed = -1        *设置随机数作为seed,-1代表使用系统当前时间作为随机数

       seqfile = input.phy *输入多序列比对文件

      treefile = input.tre *带校准点(化石时间)的有根树文件

       outfile = out.txt   *输出文件

      mcmcfile = mcmc.txt  *输出的mcmc信息文件,可用Tracer软件查看


         ndata = 1   * 输入的多序列比对的数据区域的数量;多个数据phy格式合并

       seqtype = 2   * 设置多序列比对数据类型;0:核酸数据;1:密码子比对数据;2:氨基酸数据;

       usedata = 3   * 0: no data; 1:seq like; 2:use in.BV; 3: out.BV

                     * 是否利用多序列比对数据;

                     * 0: no data,不使用,不会进行likelihood估算,会快速得到mcmc树,但分歧时间不可用;

                     * 1:seq like,使用多序列比对数据进行likelihood估算,正常进行mcmc; usedata=1时model无法选择;

                     * 2:use in.BV, 进行正常的approximation likelihood分析,不读取多序列比对数据,直接读取当前目录的in.BV文件,in.BV是由usedata = 3时生成的out.BV重命名得来;此外,由于程序BUG,当设置usedata = 2时,一定要在改行参数后加 *,否则程序报错 Error: file name empty..;

                     * 3:out.BV,程序利用多序列比对数据调用baseml/codeml命令对数据进行分析,生成out.BV文件。由于mcmctree调用baseml/codeml进行估算的参数设置可能不太好(特别时对蛋白序列进行估算时),推荐自己修改软件自动生成的baseml/codeml配置文件,然后再手动运行baseml/codeml命令,再整合其结果文件为out.BV文件。


         clock = 2       * 设置分子钟算法,1: global clock,表示所有分支进化速率一致; 2: independent rates,各分支的进化速率独立且进化速率的对数log(r)符合正态分布; 3,correlated rates方法,和方法2类似,但是log(r)的方差和时间t相关。

                         *       TipDate = 1 100  *当外部节点由取样时间时使用该参数进行设置,同时该参数也设置了时间单位。具体数据示例请见examples/TipData文件夹。

        RootAge = '<2'  * constraint on root age, used if no fossil for root.设置root节点的分歧时间,一般设置一个最大值。


         model = 4      * models for DNA:

                        *  0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85;*设置碱基替换模型;当设置usedata = 1时,model不能使用超过4的模型,所以usedata = 1时用model = 4;usedata不等于1时,用model = 7,即GTR模型;

                        * models for codons:

                        *  0:one 恒定速率模型, 1:b 中性模型, 2:2 or more dN/dS ratios for branches 选择模型。

                        * models for AAs or codon-translated AAs:

                        *  0:poisson, 1:proportional, 2:Empirical, 3:Empirical+F

                        *  6:FromCodon, 7:AAClasses, 8:REVaa_0, 9:REVaa(nr=189)

         alpha = 0.5    * alpha for gamma rates at sites;*核酸序列中不同位点,其进化速率不一致,其变异速率服从GAMMA分布。一般设置GAMMA分布的alpha值为0.5。若该参数值设置为0,则表示所有位点的进化速率一致。此外,当userdata = 2时,alpha、ncatG、alpha_gamma、kappa_gamma这4个GAMMA参数无效。因为userdata = 2时,不会利用到多序列比对的数据。

         ncatG = 5      * No. categories in discrete gamma;设置离散型GAMMA分布的categories值。


     cleandata = 0    * remove sites with ambiguity data (1:yes, 0:no)?

                      * 设置是否移除不明确字符(N、?、W、R和Y等)或含以后gap的列后再进行数据分析: 0,不需要,但在序列两两比较的时候,还是会去除后进行比较;


       BDparas = 1 1 0.1   * birth, death, sampling;*设置出生率、死亡率和取样比例。若输入有根树文件中的时间单位发生改变,则需要相应修改出生率和死亡率的值。例如,时间单位由100Myr变换为1Myr,则要设置成'.01 .01 0.1'。

   kappa_gamma = 6 2       * gamma prior for kappa;设置kappa(转换/颠换比率)的GAMMA分布参数。

   alpha_gamma = 1 1       * gamma prior for alpha;设置GAMMA形状参数alpha的GAMMA分布参数。只对usedata=1时起作用,其他2,3不起作用


   rgene_gamma = 2 20 1 0   * au bu a prior, gamma prior for rate for genes;物种替换速率可查文献;设置序列中所所有位点平均[碱基/密码子/氨基酸]替换率的Dirichlet-GAMMA分布参数:alpha=2、beta=20、初始平均替换率为每100million年(取决于输入有根树文件中的时间单位)1个替换。若时间单位由100Myr变换为1Myr,则要设置成'2 2000 1'。总体上的平均进化速率为:2 / 20 = 0.1 个替换 / 每100Myr,即每个位点每年的替换数为 1e-9。

  sigma2_gamma = 1 10 1    * gamma prior for sigma^2     (for clock=2 or 3);设置所有位点进化速率取对数后方差(sigma的平方)的Dirichlet-GAMMA分布参数:alpha=1、beta=10、初始方差值为1。当clock参数值为1时,表示使用全局的进化速率,各分枝的进化速率没有差异,即方差为0,该参数无效;当clock参数值为2时,若修改了时间单位,该参数值不需要改变;当clock参数值为3时,若修改了时间单位,该参数值需要改变。


      finetune = 1: .1 .1 .1 .1 .01 .1  * times, rates, mixing, paras, RateParas;冒号前的值设置是否自动进行finetune,一般设置成1,然程序自动进行优化分析;冒号后面设置各个参数的步进值:times, musigma2, rates, mixing, paras, FossilErr。由于有了自动设置,该参数不像以前版本那么重要了,可能以后会取消该参数。


         print = 1      *设置打印mcmc的取样信息:0,不打印mcmc结果;1,打印除了分支进化速率的其它信息(各内部节点分歧时间、平均进化速率、sigma2值);2,打印所有信息。

        burnin = 2000   *将前2000次迭代burnin后,再进行取样(即打印出该次迭代估算的结果信息,各内部节点分歧时间、平均进化速率、sigma2值和各分支进化速率等)。

      sampfreq = 10      *每10次迭代则取样一次

       nsample = 100000  *当取样次数达到该次数时,则取样结束,同时结束程序。


*** Note: Make your window wider (100 columns) when running this program.



">mcmctree2.ctl





mcmctree  mcmctree2.ctl


### 将生成的out.BV 文件重命名成 in.BV

mv out.BV in.BV


### 将mcmctree.ctl配置文件usedata 设置为2 ,运行mcmctree


echo "

          seed = -1        *设置随机数作为seed,-1代表使用系统当前时间作为随机数

       seqfile = input.phy *输入多序列比对文件

      treefile = input.tre *带校准点(化石时间)的有根树文件

       outfile = out.txt   *输出文件

      mcmcfile = mcmc.txt  *输出的mcmc信息文件,可用Tracer软件查看


         ndata = 1   * 输入的多序列比对的数据区域的数量;多个数据phy格式合并

       seqtype = 2   * 设置多序列比对数据类型;0:核酸数据;1:密码子比对数据;2:氨基酸数据;

       usedata = 2   * 0: no data; 1:seq like; 2:use in.BV; 3: out.BV

                     * 是否利用多序列比对数据;

                     * 0: no data,不使用,不会进行likelihood估算,会快速得到mcmc树,但分歧时间不可用;

                     * 1:seq like,使用多序列比对数据进行likelihood估算,正常进行mcmc; usedata=1时model无法选择;

                     * 2:use in.BV, 进行正常的approximation likelihood分析,不读取多序列比对数据,直接读取当前目录的in.BV文件,in.BV是由usedata = 3时生成的out.BV重命名得来;此外,由于程序BUG,当设置usedata = 2时,一定要在改行参数后加 *,否则程序报错 Error: file name empty..;

                     * 3:out.BV,程序利用多序列比对数据调用baseml/codeml命令对数据进行分析,生成out.BV文件。由于mcmctree调用baseml/codeml进行估算的参数设置可能不太好(特别时对蛋白序列进行估算时),推荐自己修改软件自动生成的baseml/codeml配置文件,然后再手动运行baseml/codeml命令,再整合其结果文件为out.BV文件。


         clock = 2       * 设置分子钟算法,1: global clock,表示所有分支进化速率一致; 2: independent rates,各分支的进化速率独立且进化速率的对数log(r)符合正态分布; 3,correlated rates方法,和方法2类似,但是log(r)的方差和时间t相关。

                         *       TipDate = 1 100  *当外部节点由取样时间时使用该参数进行设置,同时该参数也设置了时间单位。具体数据示例请见examples/TipData文件夹。

        RootAge = '<2'  * constraint on root age, used if no fossil for root.设置root节点的分歧时间,一般设置一个最大值。


         model = 4      * models for DNA:

                        *  0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85;*设置碱基替换模型;当设置usedata = 1时,model不能使用超过4的模型,所以usedata = 1时用model = 4;usedata不等于1时,用model = 7,即GTR模型;

                        * models for codons:

                        *  0:one 恒定速率模型, 1:b 中性模型, 2:2 or more dN/dS ratios for branches 选择模型。

                        * models for AAs or codon-translated AAs:

                        *  0:poisson, 1:proportional, 2:Empirical, 3:Empirical+F

                        *  6:FromCodon, 7:AAClasses, 8:REVaa_0, 9:REVaa(nr=189)

         alpha = 0.5    * alpha for gamma rates at sites;*核酸序列中不同位点,其进化速率不一致,其变异速率服从GAMMA分布。一般设置GAMMA分布的alpha值为0.5。若该参数值设置为0,则表示所有位点的进化速率一致。此外,当userdata = 2时,alpha、ncatG、alpha_gamma、kappa_gamma这4个GAMMA参数无效。因为userdata = 2时,不会利用到多序列比对的数据。

         ncatG = 5      * No. categories in discrete gamma;设置离散型GAMMA分布的categories值。


     cleandata = 0    * remove sites with ambiguity data (1:yes, 0:no)?

                      * 设置是否移除不明确字符(N、?、W、R和Y等)或含以后gap的列后再进行数据分析: 0,不需要,但在序列两两比较的时候,还是会去除后进行比较;


       BDparas = 1 1 0.1   * birth, death, sampling;*设置出生率、死亡率和取样比例。若输入有根树文件中的时间单位发生改变,则需要相应修改出生率和死亡率的值。例如,时间单位由100Myr变换为1Myr,则要设置成'.01 .01 0.1'。

   kappa_gamma = 6 2       * gamma prior for kappa;设置kappa(转换/颠换比率)的GAMMA分布参数。

   alpha_gamma = 1 1       * gamma prior for alpha;设置GAMMA形状参数alpha的GAMMA分布参数。只对usedata=1时起作用,其他2,3不起作用


   rgene_gamma = 2 20 1 0   * au bu a prior, gamma prior for rate for genes;物种替换速率可查文献;设置序列中所所有位点平均[碱基/密码子/氨基酸]替换率的Dirichlet-GAMMA分布参数:alpha=2、beta=20、初始平均替换率为每100million年(取决于输入有根树文件中的时间单位)1个替换。若时间单位由100Myr变换为1Myr,则要设置成'2 2000 1'。总体上的平均进化速率为:2 / 20 = 0.1 个替换 / 每100Myr,即每个位点每年的替换数为 1e-9。

  sigma2_gamma = 1 10 1    * gamma prior for sigma^2     (for clock=2 or 3);设置所有位点进化速率取对数后方差(sigma的平方)的Dirichlet-GAMMA分布参数:alpha=1、beta=10、初始方差值为1。当clock参数值为1时,表示使用全局的进化速率,各分枝的进化速率没有差异,即方差为0,该参数无效;当clock参数值为2时,若修改了时间单位,该参数值不需要改变;当clock参数值为3时,若修改了时间单位,该参数值需要改变。


      finetune = 1: .1 .1 .1 .1 .01 .1  * times, rates, mixing, paras, RateParas;冒号前的值设置是否自动进行finetune,一般设置成1,然程序自动进行优化分析;冒号后面设置各个参数的步进值:times, musigma2, rates, mixing, paras, FossilErr。由于有了自动设置,该参数不像以前版本那么重要了,可能以后会取消该参数。


         print = 1      *设置打印mcmc的取样信息:0,不打印mcmc结果;1,打印除了分支进化速率的其它信息(各内部节点分歧时间、平均进化速率、sigma2值);2,打印所有信息。

        burnin = 2000   *将前2000次迭代burnin后,再进行取样(即打印出该次迭代估算的结果信息,各内部节点分歧时间、平均进化速率、sigma2值和各分支进化速率等)。

      sampfreq = 10      *每10次迭代则取样一次

       nsample = 100000  *当取样次数达到该次数时,则取样结束,同时结束程序。


*** Note: Make your window wider (100 columns) when running this program.



">mcmctree3.ctl






mcmctree  mcmctree3.ctl


## 最终结果仍为FigTree.tre

  • 发表于 2023-02-27 20:46
  • 阅读 ( 2491 )
  • 分类:基因组学

0 条评论

请先 登录 后评论
omicsgene
omicsgene

生物信息

652 篇文章

作家榜 »

  1. omicsgene 652 文章
  2. 安生水 321 文章
  3. Daitoue 167 文章
  4. 生物女学霸 120 文章
  5. 红橙子 78 文章
  6. CORNERSTONE 72 文章
  7. rzx 63 文章
  8. xun 61 文章