分子钟假说

基因序列中密码子随着时间的推移以几乎恒定的比例相互替换，即具有恒定的演化速率，两个物种之间的遗传距离将与物种的分化时间成正比。我们通常采用单拷贝基因家族中的四重简并位点来估算分子钟（替换速率）以及物种间的分化时间，密码子中的四重简并位点由于第三碱基不改变所编码的氨基酸属于中性进化，其中中性替换速率一般用每个位点每年的变异数来衡量。

1. 分子钟应当被看作是氨基酸或核苷酸突变的随机性所导致的随机钟
2. 分子钟假说允许不同蛋白质间进化速率不同，以不同的速率跳动
3. 速率恒定性未必对所有物种适用

全局分子钟模型：序列间的期望距离随分歧时间线性增加。
局部分子钟模型：对分支赋予不同的进化速率

常用软件：

image.png

分析实操

使用软件：PAML 软件包中的mcmctree 工具
准备文件：

1. 多序列比对文件，Phylip 格式
2. 有化石时间标定的进化树
3. mcmctree 配置文件

使用氨基酸序列进行估计，使用orthofinder 的进化树及单拷贝基因数据做输入。

#将orthofinder的进化树SpeciesTree_rooted.txt和多序列比对结果SpeciesTreeAlignment.fa拷贝过来


#重命名 
mv SpeciesTree_rooted.txt  input.tree
#fasta格式转phylip格式
trimal -in SpeciesTreeAlignment.fa -out  Sequence.phy -phylip_paml

step1 估算位点替换率

mcmctree mcmctree.ctl

• mcmctree.ctl 配置文件

          seed = -1
       seqfile = ./Sequence.phy
      treefile = ./input.tree
       outfile = out

         ndata = 1
       seqtype = 2  * 0: nucleotides; 1:codons; 2:AAs
       usedata = 3    * 0: no data; 1:seq like; 2:use in.BV; 3: out.BV
         clock = 3    * 1: global clock; 2: independent rates; 3: correlated rates
       RootAge = <1.0  * safe constraint on root age, used if no fossil for root.

         model = 0    * 0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85
         alpha = 0    * alpha for gamma rates at sites
         ncatG = 5    * No. categories in discrete gamma

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

       BDparas = 1 1 0    * birth, death, sampling
   kappa_gamma = 6 2      * gamma prior for kappa
   alpha_gamma = 1 1      * gamma prior for alpha

   rgene_gamma = 2 2   * gamma prior for overall rates for genes
  sigma2_gamma = 1 10    * gamma prior for sigma^2     (for clock=2 or 3)

      finetune = 1: 0.1  0.1  0.1  0.01 .5  * auto (0 or 1) : times, musigma2, rates, mixing, paras, FossilErr

         print = 1
        burnin = 2000
      sampfreq = 2
       nsample = 20000

将生成的tmp0001.trees文件进行修改，定根（添加一对括号即可），添加化石时间
修改生成的tmp0001.ctl文件为

seqfile = tmp0001.txt
treefile = tmp0001.trees
outfile = tmp0001.out
noisy = 3
seqtype = 2
model = 0
Small_Diff = 0.1e-6
getSE = 0
method = 1
clock = 1

再次运行

codeml  tmp0001.ctl

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

#调整mcmctree.ctl 中的rgene_gamma 第二个参数b, 使得a/b 约等于s上一步tmp0001.out文件得到的替换率， 将usedata 设置为3 再运行

mcmctree  mcmctree.ctl

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

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

• 结果文件FigTree.tre