1. **前言：**目前，主流的估算tMRCA和进化速率的方法和软件在这里不再啰嗦，详见“[使用Treetime分析tMRCA和进化速率](https://www.yuque.com/wusheng/gw7a9p/xmttdb)”，本文抛转引玉，简单介绍使用R语言的**treedater**包估算**tMRCA和平均进化速率**。<br />**    软件地址：**[https://github.com/emvolz/treedater](https://github.com/emvolz/treedater)<br />**    参考文献：**《Scalable relaxed clock phylogenetic dating》

1.文件准备

（1）.nwk格式的树文件；
（2）.csv 格式的序列采样时间文件；
（3）对齐后的序列文件，用来查看长度
其中，序列采样时间.csv文件里内容示例如下：使用treedater估算tMRCA和平均进化速率 - 图1
需要注意的而是：（1）.nwk文件和.csv文件里的序列名称要一致；（2）这里的date（日期）数据要求为十进制格式，不能是“年月日”那种格式。可以将“年月日”格式的日期，在beauti里，通过“tip-date”转变成这种十进制格式的日期，或者通过其他方式转换也行。

2.安装并使用treedater

假设你电脑上已经安装了R和RStudio，现在进入RStudio。
（1）第一步：切换工作目录
我们可以事先创建一个文件夹，假设为CeSi，将.nwk格式的树文件和.csv 格式的采样时间文件全部放在该文件夹下，进入RStudio后，通过点击菜单栏“Tools”→“Global Options..”，将工作目录设置为文件夹CeSi。或者也可以在RStudio中使用代码切换工作目录。

（2）估算tMRCA和进化速率
笔者使用示例数据：树文件ebola.nwk、采样时间文件ebola.metadata.csv、序列文件ebola.fasta（可以查看对齐后序列长度）
示例数据附件：ebola.zip
欲估算ebola病毒的tMRCA和进化速率，R代码和相关注释如下：

  1. library("devtools")  #调用库devtools
  2. install_github("emvolz/treedater") # 在线下载并安装treedater
  3. require(treedater)   #载入treedater
  4. require(ape)   #载入ape，为后续读取树文件做准备
  5. tre<-read.tree("ebola.nwk")
  6. Times <-  read.csv("ebola.metadata.csv",header=TRUE) #读取采样时间文件中的数值，header=TRUE表示将第一行作为表头(因为第一行没有时间信息)
  7. sts <- setNames(Times[,2], Times[,1]) #显示所提取的时间信息
  8. s=19006 #指定序列长度，本次测试序列对齐后长度为19006bp
  9. dtr <- dater( tre , sts, s) #获得tMRCA和平均突变率等信息
  10. dtr #输出结果

—————————————————————————————————————————-
题外话:
注:某次数据量比较大时，笔者的使用ubuntu上的RStudio和R进行上面分析，由于之前没有安装devtools包，install.packages(“devtools”)死活安装不上，提示libssl缺少、httr缺少，在RStudio安装libssl也安装不上，最后通过如下方式解决:
(1) 在ubuntu的命令提示符里(不是在RStudio里)输入: sudo apt-get install libssl-dev
(2)再进入RStudio，输入:

  1. install.packages("httr")
  2. install.packages("devtools")
  1. 成功!<br />----------------------------------------------------------------------------------------

运行后，输出结果如下：

  1. Phylogenetic tree with 362 tips and 361 internal nodes.
  3. Tip labels:
  4.     G3734, LIBR0605, EM098, LIBR10237, LIBR1221, LIBR10196, ...
  5. Node labels:
  6.     NA, , 0.925, 0.365, 0.846, 0.948, ...
  8. Rooted; includes branch lengths.
  10.  Time of common ancestor 
  11. 2013.99556045057 
  13.  Time to common ancestor (before most recent sample) 
  14. 2.25443954942602 
  16.  Mean substitution rate 
  17. 0.000734601146412705 
  19.  Strict or relaxed clock 
  20. relaxed 
  22.  Coefficient of variation of rates 
  23. 0.210634771756986

3.结果解读

  1.  该程序默认基于宽松(relaxed )分子钟假设来估算（可调），从上述结果可以看出，本次测试的ebola病毒数据，**其共同祖先分歧时间（tMRCA）为2013.99556045057，平均进化速率为0.000734601146412705 （/site/year）。**<br />     ~~**Coefficient of variation of rates**值没有接近0，表示这些数据还是适合宽松分子钟。~~<br />     Coefficient of variation of rates 的值为0.21，较低，表明进化速率在不同谱系之间差异不大，反馈这些数据可能适合严格分子钟（没有具体的阈值）。<br />      需要注意的是，共同祖先分歧时间，软件返回的结果是十进制，我们需要将它转化为日期格式，我们可以使用R语言的lubridate包进行转化，代码如下：
  1. library(lubridate)
  2. date_decimal(2013.99556045057)

返回结果如下：

  1. [1] "2013-12-30 09:06:34 UTC"

即2013年12月30日。

4.注意点

  1. 和 Treetime类似，treedater这个软件预先利用每条序列的采样时间进行分子钟校准来计算这些数值。所以，应用treedater估算tMRCA和平均进化速率的前提是这些序列要具有时间结构，**如果这些序列不具备时间结构，那么通过它估算的tMRCA和进化速率则没有意义。**<br />  <br />  **FAQ:**<br />  ** 如果序列不具有时间结构，怎么估算它们的tMRCA和平均进化速率？**<br />**   **如果序列不具备时间结构，Treetime以及treedater等此类通过序列的采样时间进行分子钟校准的方法已经不再适用，此时可以使用beast， 同理，beast里也不能开启Tip-date功能。<br />**   **<br />** 注：treedater完整的使用方法，可参考官网教程：**[treedater官网教程](https://www.yuque.com/wusheng/gw7a9p/gw7ldu)

BY 月明时
2019年3月19日初稿
欢迎转载，请注明出处，如有错误请欢迎与我反馈。