分化时间计算

分化时间或称分歧时间估计，依据的是分子钟假说，即氨基酸或核苷酸替代速率在进化过程中随时间或进化谱系近似地保持恒定，以某几个特定类群的化石时间作为校正，然后通过基因序列间的分歧程度以及分子钟来估算物种间的分歧时间，同时估算系统发育树上其它节点的发生时间，从而推断相关类群的起源和不同类群的分歧时间。目前，采用依据BI树和ML树估计物种分歧时间的程序很多，例如R8S、MCMCTREE、MULTIDIVTIME、BEAST、MEGA等，不同软件通过不同的策略将化石时间信息整合到一个系统发育树中，从而计算得到Divergence time Tree。

分子钟假说成立的条件DNA或者蛋白质序列的替代速率是恒定的。20世纪80年代以来，随着DNA序列数据快速积累，大量的证据表明：在长期进化过程中，很多类群的绝大多数基因或蛋白质的序列替换速率根本不符合分子钟假说。对于蛋白质序列，在物种适应辐射过程中，其进化速度可能会大大加快。因此，以蛋白质为基础的恒定进化速率并非理想的分子钟；对于核酸分子，不同基因的分子钟速率不同，但是我们可以认为相同物种某一片段的基因或者蛋白质的替代速率是一定的，在此基础上可以将分子钟模型分为全局分子钟（序列间的期望距离随分歧时间线性增加），局部分子钟模型（不同分支的进化速率不同）
在速率恒定的假设下，遗传距离是时间的线性函数，为了将遗传距离转化为分歧时间，至少需要一个能够提供时间信息的标定点（calibration point）。常用的校准信息可以分为：(1)已知的碱基替代速率；(2)化石校准点；(3)生物地理事件校准点；(4)二次校准点
本篇教程主要参考：https://yanzhongsino.github.io/2021/03/25/bioinfo_phylogeny_caculate.divergence.time/

1.安装软件

PAML (Phylogenetic Analysis by Maximum Likelihood)，网址：https://github.com/abacus-gene/paml
官方安装教程网址：https://github.com/abacus-gene/paml/wiki/Installation#exporting-paths-linux ，因此不过多介绍安装流程

2.准备文件

多序列比对文件，phylip格式
有化石时间标定的进化树（有根树）
mcmctree配置文件

2.1检查进化树的结构是否有问题

ncbi->taxonomy->taxonomy common tree（按照物种分类展示的树结构），输入物种的拉丁名，可以看物种的大致进化关系

2.2寻找相应物种的化石时间

在线网站：https://timetree.org/home ，网站主要记录物种的进化时间和进化关系，会定期更新，有三种搜索模式，1）输入两个物种名称，即可计算得到二者的进化时间，2）提供物种名称，给出该物种的进化事件，3）给出一组物种，他来标定化石时间，一般使用第三种
注意查看引用文件，了解使用的是分子时间还是化石时间

展示结果：
1）直接输入
左侧展示所有研究的分化时间可视化，横坐标文字部分是地质年代，数字为地质年代，曲线从左往右依次是光照，二氧化碳含量，氧含量，陨石撞击事件的陨石坑大小
右侧会给出Median Time：XXXMYA（百万年）、大致区间以及所有相关研究给出的大致时间的正态拟合
下面表格给出的是相关文献以及对应事件
2）提供物种
左侧上同
右侧展示出物种的进化历史，界门纲目科属种出现的时间
3）上传物种名称
先显示进化关系，空心的点表示研究相对较少的时间节点，实心点表示研究较多的时间节点，*号表示物种名称被替换成收录名称
其他大致一样不过多介绍

展示出化石时间后，有两种使用方式：
1.直接下载相应的树文件，作为有化石时间标定的树文件（左侧最下角，To Newick File）
2.copy想要的化石时间，加到自己的树文件中，使用的是RANGE TIME

一般使用第二种，将之前建树所得的newick树文件中的枝长自展值删去，添加上化石时间范围以及物种数和树的数量（中间用空格隔开），添加之后如下：

1
2
3

25 1
((LZ,WX),(JB,(SJ,(YX,(WYZ,(SL,(SD,((DX,(PT,((MGY,NNJ),(YJ,(nigra,CH)))),(MDL,((EL,XZ)'>1.14<1.30',(EZQ,(HGM,((TN,HYM),(XH,(WC,JN))))))))))))))));
((SL,(WYZ,(YX,(SJ,((WX,LZ),JB))))),((((XZ,EL)'>1.14<1.30',(EZQ,(HG,((HYM,TN),(XH,(JN,WC))))))'>1.3<1.42',MDL),(DX,(PT,((NNJ,MGY),(YJ,(CH,nigra)))))),SD);

注意：1代表100MYA，以及之前建树是在ITOL中定的根，一定得是有根树
删去自展值可以使用

1	sed 's/:[^,)(]\+//' Speciestree.txt

2.3获得phylip格式多序列比对文件

1
2

cat SpeciesTreeAlignment.fa |tr '\n' '\t'|sed 's/>/\n/g' |sed 's/\t/      /'|sed 's/\t//g'| awk 'NF > 0' > supergene.phy.tmp
awk '{print "  "NR"  "length($2)}'  supergene.phy.tmp|tail -n 1 | cat -  supergene.phy.tmp >  supergene.phy

cat SpeciesTreeAlignment.fa |tr ‘\n’ ‘\t’ (将换行符替换为制表符) |
sed ‘s/>/\n/g’ (将每个序列名前面的>符号替换为换行符) |
sed ‘s/\t/ /‘ (将每行第一个的制表符替换为多个空格) |
sed ‘s/\t//g’ (删除剩余的制表符,使序列连成一条线) |
awk ‘NF > 0’ (删除空行)> supergene.phy.tmp (临时保存)
awk ‘{print “ “NR” “length($2)}’ supergene.phy.tmp (计算序列长度)|
tail -n 1 (返回最后一行,包括了行数与序列长度) |
cat - supergene.phy.tmp (-为上一管道数据) > supergene.phy(最终文件)（代码来源：https://www.jianshu.com/p/9f46341d2906 ）

或者使用软件fast2phy(https://github.com/davidmnoriega/fast2phy )
亦或者使用

1	trimal -in SpeciesTreeAlignment.fa -out supergene.phy -phylip_paml

如果有多个区域的序列，比如exon和intron，LSC、SSC和IR，不同的基因，密码子的第一二三位，需要不同的模型分开估算，那可以把各自区域分别align之后制作多个phy文件，再合并到一起，用空行隔开，组成input.phy文件。（此时mcmctree.ctl的ndata值为区域的个数）

#OG.length包含两列，一列OGID，一列比对后序列长度
cat ./OG.length | while read line
do
        sample_id=$(echo $line |awk '{print $1}') #获取OG.ID
        sample_a=$(echo $line |awk '{print $2}') #获取序列长度
        sed "s/_.*//g" ../singlegenetree/${sample_id}/${sample_id}.mafft.pep|seqkit seq -w 0|sed -E ":a;N;s/\n/ /g;ta" |sed "s/ >/\n/g" |sed "s/>//g"|sed "s/ /  /g"|sed "1i\5  ${sample_a}" |sed '1i\ ' > ./phy/${sample_id}.phy #把上一步获取的${sample_id}.mafft.pep改为phylip格式，并在首行添加空格行（为了合并后每个OG用空行隔开）

done
cat ./phy/*phy >input.phy #合并所有phy文件为一个.phy

（代码来源：https://yanzhongsino.github.io/2021/03/25/bioinfo_phylogeny_caculate.divergence.time/ ）

2.4mcmctree配置文件

mcmctree配置文件为mcmctree.ctl，具体参数及介绍如下

seed = -1 *设置随机数作为seed，-1代表使用系统当前时间作为随机数
seqfile = supergene.phy *输入多序列比对文件
treefile =MCMTREE.tree *带校准点（化石时间）的有根树文件
outfile = mcmc.out *输出文件
mcmcfile = mcmc.txt *输出的mcmc信息文件，可用Tracer软件查看

seqtype = 0  * 设置多序列比对数据类型；0：核酸数据；1：密码子比对数据；2：氨基酸数据；
usedata = 3
* 是否利用多序列比对数据；
* 0: no data不使用，不会进行likelihood估算，会快速得到mcmc树，但分歧时间不可用; 
* 1:seq like，使用多序列比对数据进行likelihood估算，正常进行mcmc; usedata=1时model无法选择；
* 2:normal 进行正常的approximation likelihood分析，不读取多序列比对数据，直接读取当前目录的in.BV文件，in.BV是由usedata = 3时生成的out.BV重命名得来；此外，由于程序BUG，当设置usedata = 2时，一定要在改行参数后加 *，否则程序报错 Error: file name empty..；
* 3：程序利用多序列比对数据调用baseml/codeml命令对数据进行分析，生成out.BV文件。由于mcmctree调用baseml/codeml进行估算的参数设置可能不太好（特别时对蛋白序列进行估算时），推荐自己修改软件自动生成的baseml/codeml配置文件，然后再手动运行baseml/codeml命令，再整合其结果文件为out.BV文件。

ndata = 1    * 输入的多序列比对的数据区域的数量；
clock = 2    * 设置分子钟算法，1: global clock，表示所有分支进化速率一致; 2: independent rates，各分支的进化速率独立且进化速率的对数log(r)符合正态分布; 3，correlated rates方法，和方法2类似，但是log(r)的方差和时间t相关。
*       TipDate = 1 100  *当外部节点由取样时间时使用该参数进行设置，同时该参数也设置了时间单位。具体数据示例请见examples/TipData文件夹。
RootAge = <10  * constraint on root age, used if no fossil for root.设置root节点的分歧时间，一般设置一个最大值。

model = 0    * models for DNA:
                        * 0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85；*设置碱基替换模型；当设置usedata = 1时，model不能使用超过4的模型，所以usedata = 1时用model = 4
                      * 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：经验模型
			*3：经验模型 + F
			*6：基于密码子模型
			*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时，不会利用到多序列比对的数据。model=4时一般设为0.5，否则一般设为0
ncatG = 5    * No. categories in discrete gamma；设置离散型GAMMA分布的categories值。

cleandata = 0    * remove sites with ambiguity data (1:yes, 0:no)?是否去除gap

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分布参数。

rgene_gamma = 2 20 1   * 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 .1 .1  * times, rates, mixing, paras, RateParas；冒号前的值设置是否自动进行finetune，一般设置成1，然程序自动进行优化分析；冒号后面设置各个参数的步进值：times, musigma2, rates, mixing, paras, FossilErr。由于有了自动设置，该参数不像以前版本那么重要了，可能以后会取消该参数。

print = 1  *设置打印mcmc的取样信息：0，不打印mcmc结果；1，打印除了分支进化速率的其它信息（各内部节点分歧时间、平均进化速率、sigma2值）；2，打印所有信息。 
burnin = 1000000  *将前1000000次迭代burnin后，再进行取样（即打印出该次迭代估算的结果信息，各内部节点分歧时间、平均进化速率、sigma2值和各分支进化速率等）。
sampfreq = 10  *每10次迭代则取样一次
nsample = 500000  *当取样次数达到该次数时，则取样结束，同时结束程序。

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

3.运行

完成后运行
mcmctree mcmctree.ctl

4.其他方法

4.1氨基酸序列计算分化时间

一般不建议使用氨基酸序列去计算分化时间，因为一般需要计算氨基酸替换速率，并选择相应的模型，而使用核酸序列相对简单，但是如果想要直接使用之前Orthofinder的结果，也就是氨基酸序列进行分化时间的计算，则需要进行一些处理

#按照上面的处理之后先计算氨基酸替换速率，修改tmp0001.ctl的内容
cp tmp0001.ctl codeml.ctl
#添加参数，增加gamma分布的先验数据
echo "clock =1" >>codeml.ctl

#然后修改一些参数
model = 2 * model=0使用的是简单的一致模型计算氨基酸替换速率，一般建议使用2，使用WAG模型进行计算
aaRatefile = /home/walnut168/mnt/wf/program/paml-4.10.7/dat/wag.dat * paml目录下的文件一般在paml/dat/wag.dat

修改树文件
将tmp0001.trees文件改为有根树，并添加上化石标定时间，需要添加具体时间，一般取范围的中间数就行
我的原文件为

1
2
3

  1

(LZ, WX, (JB, (SJ, (YX, (WYZ, (SL, (SD, ((DX, (PT, ((MGY, NNJ), (YJ, (nigra, CH))))), (MDL, ((EL, XZ), (EZQ, (HG, ((TN, HYM), (XH, (WC, JN)))))))))))))));

修改后为

1
2
3

  1

((LZ, WX), (JB, (SJ, (YX, (WYZ, (SL, (SD, ((DX, (PT, ((MGY, NNJ), (YJ, (nigra, CH))))), (MDL, ((EL, XZ)'@1.22', (EZQ, (HG, ((TN, HYM), (XH, (WC, JN)))))))))))))));

修改完成后运行

1	codeml codeml.ctl

就算完成查看tmp0001.out，替换速率在“Substitution rate is per time unit” 下一行
完成后修改mcmtree.ctl文件中的rgene_gamma = 2 20，根据氨基酸替换速率修改，例如我的氨基酸替换速率为0.31，则大致修改为rgene_gamma = 2 6
计算完成之后重新运行mcmctree mcmctree.ctl，重新生成out.BV文件
后续操作等于重新进行一次分化时间的计算，最后获得目标文件

注意，运行时，程序运行日志中有进度，5%-100%后面第一列数字为模型接受程度，一般30%左右较好，最大波动为15%-70%，如果超过，最好修改burnin参数，并适当增加迭代次数等等，第二列为预估的分化时间，正常情况下应该是不断稳定下来，如果靠近100%处还波动较大，应该增加迭代次数，增加nsample，增加sampfreq

 5% 0.42 0.29 0.45 0.40 0.25  9.662 7.289 9.250 8.390 7.561 6.124 - 0.067 0.902 -23.42  0:02
 10% 0.38 0.30 0.36 0.41 0.24  9.856 7.443 9.520 8.836 7.780 6.542 - 0.064 0.831 -23.58  0:03
 15% 0.40 0.30 0.39 0.41 0.19  9.811 7.373 9.498 8.704 7.693 6.664 - 0.063 0.812 -23.57  0:03
 20% 0.40 0.30 0.40 0.41 0.17  9.807 7.439 9.465 8.616 7.630 6.722 - 0.063 0.786 -23.50  0:04
 25% 0.38 0.29 0.40 0.40 0.15  9.754 7.513 9.436 8.632 7.685 6.842 - 0.063 0.753 -23.49  0:05
 30% 0.38 0.30 0.40 0.39 0.15  9.778 7.412 9.467 8.672 7.747 6.909 - 0.063 0.722 -23.49  0:05
 35% 0.37 0.30 0.39 0.39 0.16  9.800 7.532 9.493 8.722 7.812 6.939 - 0.064 0.696 -23.53  0:06
 40% 0.37 0.30 0.39 0.39 0.16  9.835 7.668 9.543 8.770 7.866 6.919 - 0.064 0.677 -23.53  0:06
 45% 0.35 0.29 0.41 0.41 0.16  9.853 7.787 9.546 8.745 7.809 6.883 - 0.064 0.659 -23.49  0:07
 50% 0.34 0.29 0.40 0.42 0.16  9.811 7.843 9.483 8.703 7.756 6.855 - 0.064 0.644 -23.51  0:08
 55% 0.34 0.29 0.40 0.41 0.16  9.800 7.875 9.478 8.706 7.783 6.884 - 0.065 0.633 -23.51  0:08
 60% 0.33 0.29 0.40 0.41 0.16  9.829 7.960 9.509 8.740 7.808 6.918 - 0.064 0.621 -23.52  0:09
 65% 0.33 0.29 0.40 0.42 0.15  9.843 8.032 9.526 8.760 7.800 6.944 - 0.064 0.614 -23.52  0:09
 70% 0.33 0.28 0.40 0.42 0.15  9.852 8.055 9.540 8.775 7.824 6.978 - 0.064 0.609 -23.52  0:10
 75% 0.32 0.28 0.40 0.42 0.15  9.866 8.113 9.560 8.789 7.842 6.991 - 0.064 0.604 -23.53  0:10
 80% 0.32 0.28 0.41 0.42 0.15  9.886 8.172 9.576 8.803 7.832 6.993 - 0.064 0.596 -23.54  0:11
 85% 0.31 0.28 0.40 0.43 0.15  9.891 8.207 9.587 8.825 7.832 6.981 - 0.064 0.586 -23.54  0:12
 90% 0.31 0.28 0.40 0.43 0.16  9.907 8.214 9.609 8.849 7.829 6.972 - 0.064 0.580 -23.55  0:12
 95% 0.32 0.28 0.40 0.43 0.16  9.919 8.234 9.621 8.856 7.843 6.968 - 0.064 0.572 -23.56  0:13
100% 0.32 0.28 0.41 0.43 0.16  9.924 8.218 9.627 8.854 7.837 6.950 - 0.064 0.565 -23.55  0:13

除此之外还需要观察mcmc.out文件最后部分，观察t_n部分，该部分为内部节点的分化时间，观察一下是否波动不大

Posterior mean (95% Equal-tail CI) (95% HPD CI) HPD-CI-width

t_n26          9.7170 ( 8.0432, 10.8710) ( 8.0329, 10.8514)  2.8185 (Jnode 48)
t_n27          8.8041 ( 6.8167, 10.2398) ( 6.8624, 10.2803)  3.4178 (Jnode 47)
t_n28          9.4250 ( 7.7433, 10.6173) ( 7.8992, 10.7086)  2.8094 (Jnode 46)
t_n29          8.6206 ( 6.9800,  9.8600) ( 7.1669, 10.0149)  2.8481 (Jnode 45)
t_n30          7.3818 ( 5.9074,  8.4911) ( 5.9742,  8.5341)  2.5599 (Jnode 44)
t_n31          6.1850 ( 4.7974,  7.4212) ( 4.8173,  7.4297)  2.6125 (Jnode 43)
t_n32          5.6159 ( 4.2569,  6.9734) ( 4.2563,  6.9714)  2.7152 (Jnode 42)
t_n33          4.8949 ( 3.7075,  6.2810) ( 3.5774,  6.0312)  2.4538 (Jnode 41)
t_n34          4.3230 ( 3.1748,  5.2996) ( 3.2144,  5.3294)  2.1149 (Jnode 40)
t_n35          3.1460 ( 2.3840,  3.8388) ( 2.3812,  3.8352)  1.4540 (Jnode 39)
t_n36          2.5618 ( 1.9859,  3.1193) ( 1.9480,  3.0724)  1.1243 (Jnode 38)
t_n37          2.1956 ( 1.6387,  2.6295) ( 1.6494,  2.6355)  0.9861 (Jnode 37)
t_n38          1.6816 ( 0.9820,  2.2494) ( 1.0505,  2.2811)  1.2306 (Jnode 36)
t_n39          1.7495 ( 1.2646,  2.2873) ( 1.2470,  2.2652)  1.0182 (Jnode 35)
t_n40          0.1633 ( 0.0571,  0.3442) ( 0.0426,  0.2922)  0.2496 (Jnode 34)
t_n41          3.7507 ( 2.8133,  4.5681) ( 2.8132,  4.5679)  1.7547 (Jnode 33)
t_n42          2.9491 ( 2.2476,  3.6595) ( 2.2476,  3.6593)  1.4117 (Jnode 32)
t_n43          1.1929 ( 1.1358,  1.2913) ( 1.1322,  1.2855)  0.1533 (Jnode 31)
t_n44          0.7227 ( 0.4832,  1.0478) ( 0.4510,  0.9966)  0.5456 (Jnode 30)
t_n45          0.3438 ( 0.2135,  0.4484) ( 0.2025,  0.4341)  0.2315 (Jnode 29)
t_n46          0.3157 ( 0.1971,  0.4061) ( 0.1916,  0.3997)  0.2082 (Jnode 28)
t_n47          0.2863 ( 0.1536,  0.3715) ( 0.1659,  0.3797)  0.2139 (Jnode 27)
t_n48          0.1475 ( 0.0932,  0.2344) ( 0.0897,  0.2286)  0.1389 (Jnode 26)
t_n49          0.1210 ( 0.0707,  0.2043) ( 0.0671,  0.1954)  0.1283 (Jnode 25)
mu             0.0954 ( 0.0759,  0.1246) ( 0.0744,  0.1216)  0.0472
sigma2         0.2907 ( 0.1821,  0.4506) ( 0.1740,  0.4368)  0.2627
lnL          -23.5714 (-34.0770, -15.1300) (-33.1800, -14.4370) 18.7430

一般近似似然法计算就是会有较大的误差，所以还是推荐使用usedata = 1

4.2 4D位点计算分化时间

什么是4D位点？一个遗传密码子通常由三个核苷酸构成，从左到右依次为第一个位点、第二个位点、第三个位点。如果密码子的某个位点上无论是哪种核苷酸，均编码同样的氨基酸，则称这个位点为 4 倍简并位点，这个位点一般是密码子的第三位。4DTV( four-fold synonymous (degenerative) third-codon transversion) 表示4D位点（fourfold degenerate site）上发生颠换（嘌呤突变为嘧啶或者嘧啶突变为嘌呤）的位点替换率。

使用单拷贝基因家族的4D位点计算分化时间，则还需要准备以下文件

1）CDS文件
2）单拷贝基因家族列表
Orthogroups_SingleCopyOrthologues.txt（基因家族ID）
3）单拷贝基因家族的序列
4）基因家族统计列表
Orthogroups.tsv（基因ID）

首先根据单拷贝基因家族列表提取出基因id

1	less Orthogroups_SingleCopyOrthologues.txt \| while read aa ;do grep $aa Orthogroups.tsv ;done > id.txt

或者

1
2
3

awk '{if(NR==FNR){A[$1]=1}else{ if($1 in A || FNR==1){print $0}} }' \
Orthogroups_SingleCopyOrthologues.txt Orthogroups.tsv \
> single_copy.txt

或者

1	awk '{if(NR==1){f=NF}; if( $0 !~ /,/ && NF==f ){print $0}}' Orthogroups.tsv > single_copy.txt

在每个ID前加上物种前缀

1
2
3

awk '{if(NR==1){ n=split($0,A, "\t")} else { for(i=2; i<= n; i++){ \
printf A[i]"_"$i"\t" }; printf "\n" } }' single_copy.txt | \
sed 's/\s\+$//' > single_copy.txt.change

给每一个cds序列添加上物种标签

for i in cds/*.fa;do
base=$(basename "$i" .fa);
sed "s/^>/>${base}_/" "$i" >"$base.fa.change";
done

protein序列进行类似操作

合并cds和protein文件

1 2	cat changecds/.fa.change > all.cds.fa cat changepep/.fa.change >all.pep.fa

使用ParaAT（https://github.com/wonaya/ParaAT ）提取相应的cds文件

echo "6" > proc.txt
ParaAT.pl -h single_copy.txt.change -a all.pep.fa -n all.cds.fa -o Para_out -p proc.txt

-h single_copy.txt \ #基因家族信息，每一行一个家族
-a all.pep.fa \ #全部蛋白序列
-n all.cds.fa \ #全部cds序列
-o Para_out \ #输出结果目录
-p proc.txt #指定线程数文件

注意protein文件和cds文件名称必须与single_copy.txt.change一致，为保证这一点可以进行以下操作

#简化序列名
sed -i '/^>/ s/ .*//' changepep/*.fa.change
#直接将cds翻译成protein保证名称一致
for i in cds/*.fa ;do
base=$(basename "$i" .fa);
seqkit translate "$i" > pep/"$base.fa";
done

除此之外，序列名称不允许包含|，：，-等特殊字符，并且序列数也有限制，因此只推荐数量较少的进行

序列数较多可以使用Hyphy（https://github.com/veg/hyphy-analyses/tree/master/codon-msa ）具体流程较为复杂，需要重新编写脚本，不推荐使用

到此分化时间估计就完成了，可视化也是使用ITOL网站，但是由于文件格式原因，不能直接像之前的树文件一样拖入即可，因此需要进行一些修改，此部分在下一篇教程《基因家族扩张收缩》中介绍