使用 r8s 估算物种分歧时间

1. r8s 简介

r8s用于通过系统发育树估计分歧时间和分子演化速率。
软件运行需要:

一个带有枝长的有根树;
比对的序列长度。

2. r8s 下载与安装

$ wget http://loco.biosci.arizona.edu/r8s/r8s.dist.tgz
$ tar zxf r8s.dist.tgz -C /opt/biosoft/
$ mv /opt/biosoft/dist /opt/biosoft/r8s
$ echo 'PATH=$PATH:/opt/biosoft/r8s/' >> ~/.bashrx
$ source ~/.bashrc

2. r8s 的使用

软件解压缩后,其中带有一个 mannual 的 PDF 文件。

2.1 r8s 的使用方法

直接输入命令 r8s 则会进入该软件的命令行界面,推荐编写了 r8s 的脚本文件后,直接运行。运行 r8s 的方式如下:

$ r8s -b -f r8s_in.txt > r8s_out.txt

-b    运行 r8s 后推出软件的命令行界面
-f    输入的 r8s 脚本文件,该文件包含 r8s 的命令行

r8s_in.txt 的一个示例如下:

#nexus
begin trees;
tree tree_1 = [&R] ((gluc:0.46,cneg:0.54):4.3,(scer:1.07,tree:0.68):0.52);
end;

begin r8s;
blformat lengths=persite nsites=300000 ulrametric=no;
MRCA asc tree scer;
MRCA bas gluc cneg;
fixage taxon=asc age=520;
constrain taxon=bas min_age=350 max_age=410;
#divtime method=PL crossv=yes cvstart=0 cvinc=1 cvnum=18;
set smoothing=100;
divtime method=PL algorithm=TN;
showage;
describe plot=cladogram;
describe plot=chrono_description;
end;

按照上述示例中,第一部分是输入进化树数据;第二部分是运行 r8s 的命令。

2.2 r8s 命令

按照上述示例,需要依次输入上面的 r8s 命令:

blformat

length      设置树的枝长单位。若枝长单位是位点替换率,则其值为 persize,则 枝长 * 序列长度 = 替换数;若枝长单位是替换数,则该参数值为 total。默认参数是 total。
nsites      设置多序列比对的序列长度。
ultrametric 是否是超度量树,一般进化树选 no。默认参数是 no。

MRCA

该命令用来设置内部节点名。示例中设置了 tree 和 scer 的 most recent common ancestor 的节点名为 asc。

fixage

该命令用来设置 MRCA 指定的节点的分歧时间,使用该指定的分歧时间作为校正,来预测其它节点的分歧时间。
r8s 需要至少有一个内部节点进行 fixage。

constrain

该命令用来约束 MRCA 指定的节点的分歧时间,设置该节点允许的最大或最小分歧时间。

divtime

该命令用来进行分歧时间和替换速率计算。总共有 4 种计算方法(LF | NPRS | PL)和 3 种数学算法(Powell | TN | QNEWT)。 一般常用且较优,是使用 PL 和 TN。
但是使用 PL 方法,则需要设置参数 smoothing 的值。 通过设置多个 smoothing 的值来进行一些计算,选择最优的值即可。一般情况下,该值位于 1~1000 能得到一个最佳(最低)的分值。

divtime method=PL crossv=yes cvstart=0 cvinc=1 cvnum=18;
上述命令,是设置 smoothing 的值从 1, 10, 100, 1000 ... 1e17, 来计算,最后得到最佳的 smoothing 值。

若使用 fixage 对 2 个节点的分歧时间进行了固定,则可以运行命令:
divtime method=PL crossv=yes fossilfixed=yes cvstart=0 cvinc=1 cvnum=18;

若使用 fixage 对 1 个节点进行分歧时间固定,同时使用 constrain 对 2 个节点进行了约束,则可以运行命令:
divtime method=PL crossv=yes fossilconstrained=yes cvstart=0 cvinc=1 cvnum=18;

得到最优的 smoothing 值后,使用 set 进行设置,然后运行 divtime 进行分歧时间和替换速率的计算。

showage

使用该命令能得到各个节点的分歧时间和替换速率。

describe

使用该命令得到进化树的图和文字结果。 其 plot 参数如下:

cladogram    得到分支树的图,图上有各个节点的编号,和 showage 的结果结合观察。
phylogram    得到进化树的图,枝长表示替换数。
chronogram   得到超度量树的图,枝长表示时间。
ratogram     得到树图,枝长表示替换速率。

phylo_description     得到树的ASCII文字结果,枝长表示替换数。
chrono_description    得到树的ASCII文字结果,枝长表示时间。
rato_description      得到树的ASCII文字结果,枝长表示替换速率。

node_info    得到节点的信息表格

发表评论

您的电子邮箱地址不会被公开。 必填项已用*标注

此站点使用Akismet来减少垃圾评论。了解我们如何处理您的评论数据