絮语：

病毒（或其他病原）的进化速率（rate）和时间尺度（timescale）可以通过对具有时间结构（Time-structured）的序列进行系统发育分析获得，应用病毒序列的采样时间进行分子钟校准前，通常需要进行日期随机化检验（date-randomization test）判断分析的数据集中是否有足够的时间信号（Temporal signal）可用于校准。

1. 两个检验标准（CR1和CR2）：通过原始（正确的）采样时间进行校准估算获得的进化速率平均值（Mean）与应用日期随机化后的数据获得的进化速率95%置信区间（CI）进行比较（如下图中的小红点与黑色线段的上下限）。如果两者之间没有重合，则通过相对宽松的标准1（即：CR1）；通过原始（正确的）采样时间估算获得进化速率的95%CI与随机化后的数据获得进化速率的任何区间（如下图的红色虚线段与黑色线段的任意一个值）均没有重合，则能通过更为严格的标准2（即：CR2）；

2. 两个不同方法： 标准的Bayesian dating permutation tests 和相同采样时间归组的 Clustered permutation tests。 
本文以标准的方法为示例，应用R包TipDatingBeast进行DRTs简明图解分析。

所需工具：
1. R 包：
官网：https://cran.r-project.org/ ;下载安装对应操作系统的版本；
2. 库文件 TipDatingBeast
安装R后，在界面窗口输入以下命令（红色部分）安装
> install.packages("TipDatingBeast")

图解操作：
Step 1.  配置用于BEAST运行的XML文件
应用 BEAUti 1.8.x 配置时，在"Tip"标签下，启用 Tip-date 功能，将每条序列的采样时间根据不同规则提取出来用于校准（如下图所示），其他根据规范逐一配置完成，最好导出配置好的xml文件（本例为PMMoV_DRTs.xml），注意该xml文件中的时期是正确采样时间（非随机化的）；

Step 2. 随机化 XML 文件的采样日期

打开 R-console （Windows下推荐 R-gui.exe ）前，先将工作目录改为 PMMoV_DRTs.xml 所在的目录，可以通过“Misc”菜单下的“Change Working Directory” 修改工作目录，也可以通过命令 setwd()来修改当前工作目录（如：setwd("/Users/Raindy/Documents/R/DRTs")）。

随后输入下列命令分别载入TipDatingBeast文件和生成20个日期随机化的xml文件：

>library(TipDatingBeast)                      
# 载入TipDatingBeast库
>RandomDates(name="PMMoV_DRTs",reps=20,writeTrees=F)
# 将目录内名称为PMMoV_DRTs的xml文件进行日期随机化，并生成20个不同的重复
# 注意只需要文件名（不需要扩展名xml）

当出现 Replicates done: 20 时，20个 日期随机化后的xml 文件已经在工作目录生成。


Step 3. 应用BEAST分别运行前两步得到的*.xml，获得运行日志文件*.log
在大批量运行这些xml文件前，建议先随机抽个随机化的xml文件，看看程序是否会存在如下类似报错信息。

SEVERE: Parsing error - poorly formed XML (possibly not an XML file):
The string "--" is not permitted within comments.
java.lang.RuntimeException: Terminate
解决办法：
用文本编辑器搜索字符串    <!-- write tree log to file ，将它后一行的 <!--  的字符串删除后删除，保存xml 文件

Step 4. 应用 TipDatingBeast 绘制DTRs结果

将 Step 3运行获得的21个 *.log 复制到当前工作目录，输入以下命令，对DRTs结果进行绘制：
>library(TipDatingBeast)
# 载入TipDatingBeast库
>PlotDRT(name="PMMoV_DRTs",reps=20,burnin=0.1)
# 对PMMoV_DRTs.log和PMMoV_DRTs.rep1.log，PMMoV_DRTs.rep2.log，....PMMoV_DRTs.rep20.log的参数进行绘制，可以输入VAR查看变量
#示例数据最佳分子钟模型为relax clock with lognormal distribution，所以对应的参数VAR为meanRate


绘制完成，工作目录会生成三个文件，meanRate.stats.csv和两个PDF文件，均为DRTs的结果，一个是原数绘制的，另一个是数值取Log10后绘制的结果。

Step 5. 结果美化及解读
可以使用R对meanRate.stats.csv进行绘制，效果图如下图Figure  S2B 所示，DRTs结果表明该数据集通过DRTs检验，具有时间信号，其采样时间可以用于分子钟校准。


Raindy注：

示例图详见参考文献5: https://www.sciencedirect.com/science/article/pii/S0168170218303563

参考文献：

1. Duchêne, S., Duchêne, D., Holmes, E.C., Ho, S.Y.W., 2015. The performance of the date-randomization test in phylogenetic analyses of time-structured virus data. Molecular biology and evolution 32, 1895-1906.

2. Murray, G.G.R., Wang, F., Harrison, E.M., Paterson, G.K., Mather, A.E., Harris, S.R., Holmes, M.A., Rambaut, A., Welch, J.J., 2016. The effect of genetic structure on molecular dating and tests for temporal signal. Methods in Ecology and Evolution 7, 80-89.

3. Rieux, A., Khatchikian, C.E., 2017. TipDatingBeast: an R package to assist the implementation of phylogenetic tip-dating tests using BEAST. Molecular ecology resources 17, 608-613.

4. Tong, K.J., Duchene, D.A., Duchene, S., Geohegan, J.L., Ho, S.Y.W., 2017. A comparison of mnethods for estimating substitution rates from ancient DNA sequence data. bioRxiv. 

5. Guan, X., Yang, C., Fu, J., Du, Z., Ho, S.Y.W., Gao, F., 2018. Rapid evolutionary dynamics of pepper mild mottle virus. Virus research 256, 96-99.