使用Treetime分析tMRCA和进化速率

前言：分析一组序列的tMRCA和进化速率，主要这几类方法：①Regression of root-to-tip distances (软件TempEst 或者软件 Treetime)，②least-squares dating（LSD，速度快，适合大数据）, ③Scalable relaxed clock phylogenetic dating (Treedater)，④Bayesian inference（贝叶斯分析,通过软件beast完成）。
其中通过软件Beast分析tMRCA和进化速率在本人之前的《Beast v1.8使用手册》中已经提到，这里重点介绍Regression of root-to-tip distances（RTT）里的**treetime软件包**分析tMRCA和进化速率。

参考文献：《Estimating evolutionary rates using time-structured data: a general comparison of phylogenetic methods》，《TreeTime: Maximum-likelihood phylodynamic analysis》，《A comparison of methods for estimating substitution rates from ancient DNA sequence data 》

1.这几类方法的区别

（1）假设前提

TempEst 和 Treetime以及LSD假设序列符合严格分子钟，而BEAST默认是假设序列符合宽松分子钟。

（2）适用范围

①RTT方法的适用于进化速率低的情况；估计分子钟速率时，当物种多样性低时TreeTime和LSD比BEAST表现好；当物种多样性高时，BEAST表现好，TreeTime表现最差，物种多样性对beast的影响不大。

②序列数量对Treetime、LSD、BEAST的影响

③LSD可以处理在不同谱系之间具有明显水平速率变化的数据，LSD分析大数据集时速度快，而BEAST则非常慢。

2.Treetime的使用

（1）在线版https://treetime.biozentrum.unibas.ch/

①网站提供了示例数据，我们可以用示例数据看下效果：

①我们选择第一个“Time-tree inference”

这里需要我们准备三个文件，newick格式的树文件（通过其他软件建好的ML树）、fasta格式的对齐序列文件、时间文件，示例数据文件格式如下附件：
data.zip
需要注意的是，三个文件里的序列名称要一致。
默认是GTR模型，也可以进行其他模型设置：

我们先用网站提供的示例数据跑下看下效果：（通过点击“Example datasets”预览示例数据）

我们选择流感H3N2的HA基因，测试下，Load后点击下面的Run.
耐心等待几分钟：

②结果：（点击放大更清晰）

③结果解读

左下角的图，μ 代表进化速率，如示例数据是3.921e-3(/site/year) ，R² 代表相关系数，示例数据为0.586。在上图中，我们可以看到，随着时间的增大，不同序列的Distance to root值也大致呈线性变化，表明这些序列之间的差异和采样时间大致具有的线性关系，反馈出这些序列具有时间结构；同时，不同谱系的序列之间良好的线性关系表示在不同谱系之间进化速率差别不大，说明这些序列可能适合严格分子钟。
遗憾的是虽然图很美，似乎还提供了结果打包的zip文件，但是，zip文件里没有提供图片。没有图片倒问题不大，如果结果里提供数据也可以自己做，比如zip里的树文件倒是可以用FigTree美化下，效果也还过得去。可是，zip文件里没有上图左下角的图片的数据（这张可以看进化速率和分子钟类型的图片非常重要），这是在线版的一个非常大的bug，所以，建议用本地版的Treetime。
**

（2）在Linux上运行Treetime

FAQ：我的电脑是windows系统，没有linux怎么办?
如果你的电脑是win10系统，那就好办了，可以免费安装win10的分发版linux子系统，安装教程见下：
https://www.windows10.pro/bash-on-ubuntu-on-windows/
这种方式，你可以简单理解为在你的win10系统上多装了一个软件。

其实笔者也不想在liunx上运行这个程序，遗憾的是，windows版本的Treetime笔者没有找到和结果输出相关的参数，换句话说，windows版本的Treetime可以跑数据，但是不知道结果输出在哪里？实时界面虽然有一部分结果，但是不是想要的，于是先以linux版本的Treetime为例，以后找到了再补充windows版本的Treetime使用说明。
**Treetime这个程序是python包。**

①Linux上安装Treetime

安装方法非常简单，但是前提是你Linux系统上已经安装了pip和biopython（pip和biopython的安装方式见本人的“Linux学习笔记”知识库）
先切换到管理员账号（非win10的分发版linux子系统不需要切换到root），再在命令端输入下面命令：

pip install phylo-treetime

也可以指定版本进行安装，如安装treetime 0.5.2版本：

pip install phylo-treetime==0.5.2

系统提示安装完成后，再在python里输入import treetime检查下是否安装成功，如下：

如果在python里输入import treetime没有报错则表示treetime包已经安装成功。

②运行Treetime

在linux系统命令终端输入treetime -h 查看帮助文档（注意不是在python里）
部分截图如下：

里面有每个参数的介绍，其中指定输出结果的位置的参数就是下面这个：

比如还有指定进化模型的参数是—gtr等等。

③功能

（i）Timetrees
The to infer a timetree, i.e. a phylogenetic tree in which branch length reflect time rather than divergence, TreeTime offers implements the command:
$ treetime —aln —tree —dates
```
This command will infer a time tree, ancestral sequences, a GTR model, and optionally confidence intervals and coalescent models.
A detailed explanation is of this command with its various options and examples are available at treetime_examples/timetree.md

(ii)Rerooting and substitution rate estimation
To explore the temporal signal in the data and estimate the substitution rate (instead if full-blown timetree estimation), TreeTime implements a subcommand clock that is called as follows
$ treetime clock —tree —aln —dates —reroot least-squares
``<br />The full list if options is available by typingtreetime clock -h.<br />Instead of an input alignment,—sequence-length ` can be provided.
Documentation of additional options and examples are available at treetime_examples/clock.md

(iii)Ancestral sequence reconstruction
The subcommand
$ treetime ancestral —aln input.fasta —tree input.nwk
``<br />will reconstruct ancestral sequences at internal nodes of the input tree.<br />The full list if options is available by typingtreetime ancestral -h.<br />A detailed explanation oftreetime ancestral` with examples is available at treetime_examples/ancestral.md

(iv) Homoplasy analysis
Detecting and quantifying homoplasies or recurrent mutations is useful to check for recombination, putative adaptive sites, or contamination.
TreeTime provides a simple command to summarize homoplasies in data
$ treetime homoplasy —aln —tree
``<br />The full list if options is available by typingtreetime homoplasy -h`.
Please see treetime_examples/homoplasy.md for examples and more documentation.

(v)Mugration analysis
Migration between discrete geographic regions, host switching, or other transition between discrete states are often parameterized by time-reversible models analogous to models describing evolution of genome sequences.
Such models are hence often called “mugration” models.
TreeTime GTR model machinery can be used to infer mugration models:
$ treetime mugration —tree —states —attribute
注意：以上黑体部分只是各功能最简单的模板，具体参数请参考help文档。

④示例数据

以官网提供的示例数据包里的ebola数据为例

打开linux终端，创建文件夹out为用来输出结果的文件夹。再输入treetime的运行命令：

treetime --aln ebola.fasta --tree ebola.nwk --dates ebola.metadata.csv --outdir out

正常运行效果如下：

耐心等待。运行完成后界面有结果提示：

使用上面简单的参数实际输出结果中含有下面这些：

root_to_tip_regression.pdftimetree.pdf
（a） root_to_tip_regression.pdf是下面这个：

可以看出这些数据具有时间结构，且不同谱系序列之间良好的线性关系表明其可能适合严格分子钟。其中进化速率为7.98e-04+/- 1e-04(/site/year),tMRCA为：2013.9+/- 0.13（year）

(b)timetree.pdf是下面这个：


PDF文件是矢量图，可用软件AI进行修改微调。

(c) 觉得上面这棵时间树不好看，可以用结果里的timetree.nexus文件，使用FigTree打开，重新绘制下。

（d）ancestral_sequences.fasta里有预测的每个节点的祖先序列。

3.注意点

(1) Treetime这个软件预先利用每条序列的采样时间进行分子钟校准来计算这些数值，所以可以用来查看这些序列是否具有时间结构，如果具有时间结构，那么它估算的tMRCA和进化速率是有意义的；如果这些序列不具备时间结构，那么它估算的tMRCA和进化速率则没有意义。详情请参考文首第一篇参考文献。

(2) Linux版本的treetime 0.5.2跑数据时出现报错，信息如下:

File "/home/gong/.pyenv/versions/stackless-3.7.5/bin/treetime", line 273, in <module>
    return_code = params.func(params)
  File "/home/gong/.pyenv/versions/stackless-3.7.5/bin/treetime", line 198, in toplevel
    timetree(params)
  File "/home/gong/.pyenv/versions/stackless-3.7.5/lib/python3.7/site-packages/treetime/wrappers.py", line 539, in timetree
    fixed_pi=fixed_pi)
  File "/home/gong/.pyenv/versions/stackless-3.7.5/lib/python3.7/site-packages/treetime/treetime.py", line 175, in run
    self.make_time_tree(**tt_kwargs)
  File "/home/gong/.pyenv/versions/stackless-3.7.5/lib/python3.7/site-packages/treetime/clock_tree.py", line 333, in make_time_tree
    self.init_date_constraints(clock_rate=clock_rate, **kwargs)
  File "/home/gong/.pyenv/versions/stackless-3.7.5/lib/python3.7/site-packages/treetime/clock_tree.py", line 286, in init_date_constraints
    one_mutation=self.one_mutation, branch_length_mode=self.branch_length_mode)
  File "/home/gong/.pyenv/versions/stackless-3.7.5/lib/python3.7/site-packages/treetime/branch_len_interpolator.py", line 39, in __init__
    grid_left = mutation_length * (1 - np.linspace(1, 0.0, n_grid_points/3)**2.0)
  File "<__array_function__ internals>", line 6, in linspace
  File "/home/gong/.pyenv/versions/stackless-3.7.5/lib/python3.7/site-packages/numpy/core/function_base.py", line 113, in linspace
    num = operator.index(num)
TypeError: 'float' object cannot be interpreted as an integer

解决方案:根据上面的路径找到对应的 function_base.py这个文件，将里面的num = operator.index(num) 这行代码修改为:num = operator.index(int(num)),保存即可。

                                                                                                            2019年1月15日初稿<br />转载请注明出处。如有错误，欢迎及时与我反馈。