使用 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    得到节点的信息表格

使用 CAFE 进行基因家族扩张分析

1. CAFE 简介

CAFE (Computational Analysis of gene Family Evolution) 用于进行基因家族的扩张收缩分析。

2. CAFE 下载和安装

$ wget http://heanet.dl.sourceforge.net/project/cafehahnlab/cafe.linux.x86_64
$ wget http://downloads.sourceforge.net/project/cafehahnlab/CAFEv3.1_Manual.pdf
$ wget http://downloads.sourceforge.net/project/cafehahnlab/CAFEv3.1_Manual.doc
$ mv cafe.linux.x86_64 ~/bin/cafe
$ cafe

3. CAFE 的简单使用

CAFE 需要的输入:

1. 基因家族在各个物种中的数目。该文件内容有多行,以 tab 分割,第一行内容必须如下:
Description    ID    species1    species2    ...
上面 species1 和 species2 为物种名简称。后面的行为相应的值。 Description 的值可以空着; ID 为基因家族的名称;其它的为数字,表示该基因家族在相应物种中的数目。
该文件可以由 OrthoMCL 的结果统计得到。

2. 超度量树,枝长表示分歧时间。可以使用 BEAST 或 r8s 等生成。树中的物种名必须和第一个文件中的物种名对应。

3. 输入参数  λ 的值。该值表示在物种进化过程中,每个单位时间内基因获得与丢失的概率。可以由软件自己进行计算。

直接输入命令 cafe 则进入了软件的命令行界面,也可以将命令写入到 cafe 脚本中,直接运行。一个简单的示例如下:

#!~/bin/cafe
version
date
load -i orthoMCL2cafe.tab -t 16 -l log.txt -p 0.05
tree (((chimp:6,human:6):81,(mouse:17,rat:17):70):6,dog:93)
lambda -s -t (((1,1)1,(2,2)2)2,2)
report out

以下简单介绍如上命令

version
    显示软件的版本

date
    显示当前日期

load
    -i 输入的数据文件
    -t 设置程序运行的线程数,默认为 8
    -l 设置输出的日志文件,默认标准输出
    -p 设置 p_value 的阈值,默认为 0.01

tree 输入超度量树

lambda
    -l 设置 λ 的值
    -s 设置软件自动寻找最优的 λ 值
    -t 默认下所有分支的 λ 值是相同的,若需要不同的分支有不同的 λ 值,则用该参数进行设置。该参数的值和 tree 命令中的树的内容一致,只是去除了分歧时间,并将物种名换成了表示 λ 值的编号。其中,相同的编号表示有相同的 λ 值。例如,该参数的值为 (((1,1)1,(2,2)2)2,2) ,它表示 chimp,human 和 紧邻的分支有相同的 λ 值,其它分支有另外相同的 λ 值。

report out
    设置输出文件的前缀为 out

4. CAFE 的输出结果

CAFE 的输出结果为 out.cafe,该文件内容包含如下几部分:

4.1 Tree

第 1 行为输入的树的信息。

4.2 Lambda

第 2 行为 λ 值。

4.3 节点 ID

第 3 行为节点的 ID。同样是树的文字内容,不过给每个节点进行了编号,有利于后面的数值的对应。

4.4 每个基因家族中扩张或收缩的基因数目

第 4 行给出了一系列的 节点ID 对。CAFE 对这些 ID 对进行了基因家族扩张的统计。
第 6 行的值为平均每个基因家族中扩张的基因数目,负数表示基因家族收缩。

4.5 基因家族数目

第 7 行给出发生了扩张的基因家族数目;
第 8 行给出没有发生改变的基因家族数目;
第 9 行给出发生了收缩的基因家族数目;

4.6 每个基因家族的具体扩张与收缩情况

第 10 行之后是每个基因家族具体的扩张情况。其内容分为 4 列:

第 1 列: 基因家族 ID
第 2 列: 树的信息,其中每个物种名后面多个一个下划线和数字,该数字表示基因家族的数目。特别是每个节点的基因家族数目都计算出来了,从而知道在某一个分化过程中基因家族的扩张情况。
第 3 列: 该基因家族总体上的扩张情况的 p_value,不同物种中的基因数目差异越大,则 p 值越小。
第 4 列: 计算出了每个分枝的基因家族扩张的显著性。

根据 CAFE 的结果,可以自行编写程序提取信息。