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

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

source ./env.sh

else

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

exit

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无法选择；

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

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

model = 4 * models for DNA:

* 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)

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，不需要，但在序列两两比较的时候，还是会去除后进行比较；

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不起作用

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无法选择；

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

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

model = 4 * models for DNA:

* 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)

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，不需要，但在序列两两比较的时候，还是会去除后进行比较；

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不起作用

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

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

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

">mcmctree3.ctl

mcmctree mcmctree3.ctl

## 最终结果仍为FigTree.tre