时间树

ht • 2024-07-17 09:40 • 杂文

基因序列中的密码子随着时间的推移以几乎恒定的比例进行替换，具有恒定的演化速率，因此两个物种的遗传距离与物种分歧时间成正比。序列比对时，我们通常会利用单拷贝直系同源基因多序列比对文件用于构建进化树。为了获取更多的遗传信息，我们会将蛋白质比对转换成密码子比对，这时我们可以将这些密码子比对文件串联起来形成超级比对文件。而我们绘制时间树，需要稳定的核酸替换率，而三位密码子所受到的选择压力各不相同，因此有必要使用 PartitionFinder2 将串联而成的多序列比对文件进行分区（可以分成三个部分）。将分区的比对文件可以分别计算进化模型，RaxML，MrBayes，Beast，iqtree 都支持分区建树。而推导时间树的程序 mcmctree 和 BEAST 都支持分区信息。

PROGRAM	SUPPORT MULTITHREADING	INPUT FILE（Fossil calibration information）	TIME	EASY INSTALLATION
Mega	YES	A tree without branching information,Multi-sequence comparison file	About 1 hour	YES
R8s	NO	A tree with branching information（root）	Within 10 minutes	NO
MCMCtree	NO	A tree without branching information(root),Multi-sequence comparison file	More than four days	NO
BEAST	YES	Multi-sequence comparison file	More than four days	YES

一.mcmctree

三个输入文件（最好是核酸比对文件）：

1.序列比对文件
2.带有校准点的有根树
3.配置文件mcmctree.ctl

1.序列比对文件

2.带有校准点的有根树

1.获取树文件（树文件通过核酸比对建立更加可靠）
2.去除树文件中分支长度信息
3.访问TimeTree网站或古生物化石数据库获取物种分歧时间信息

TimeTree网站

Cl：置信区间（Confidence Interval）
MYA：百万年以前（Million Years Ago）的缩写
CI: (6.0 - 6.5 MYA)表示置信区间为从600万年前到650万年前的时间段。换句话说，这个置信区间的长度是0.5百万年，也就是50万年。

化石校准点

 7  1

((((human, (chimpanzee, bonobo)) '>.06<.08', gorilla), (orangutan, sumatran)) '>.12<.16', gibbon);

这棵带有校准点的有根树中时间单位是：一亿年 --100个百万年
上图这棵树中：人类(human)和黑猩猩(chimpanzee)之间的分歧时间在600万年-800万年之间
要注意的点是：树在根上没有化石校正点

3.配置文件mcmctree.ctl

使用PAML进行分歧时间计算 | 陈连福的生信博客 (chenlianfu.com)
估算系统树分歧时间 —— paml.mcmctree,r8s | 生信技工 (yanzhongsino.github.io)
【比较基因组】从进化树到分化时间 - 简书 (jianshu.com)

***输入输出参数：**
        seed = -1           *设置一个随机数来运行程序。若设置为-1，表示使用系统当前时间为随机数。
     seqfile = input.txt    *设置输入的多序列比对文件路径
    treefile = input.trees  *设置输入的带校准点信息有根树的文件路径 
    mcmcfile = mcmc.txt     *设置输出的mcmc信息文本文件，可以用Tracer软件查看
     outfile = out.txt      *设置输出文件路径，该文件中记录了超度量树和进化速率等参数信息。

***数据使用说明参数：**
       ndata = 3         *设置输入的多序列比对的数据个数。这里指密码子3个位置的数据。
     seqtype = 0         *设置多序列比对的数据类型：0，表示核酸数据；1，表示密码子比对数据；2，表示氨基酸数据。根据不同的数据类型，程序调用相应的baseml或codeml进行相应的参数计算。
     usedata = 1         *设置是否利用多序列比对的数据：0，表示不使用多序列比对数据，则不会进行likelihood计算，虽然能得到mcmc树且计算速度飞快，但是其分歧时间结果是有问题的；1，表示使用多序列比对数据进行likelihood计算，正常进行MCMC，是一般使用的参数; 2，进行正常的approximation likelihood分析，此时不需要读取多序列比对数据，直接读取当前目录中的in.BV文件。该文件是使用usedata = 3参数生成的out.BV文件重命名而来的。此外，由于程序BUG，当设置usedata = 2时，一定要在改行参数后加 *，否则程序报错 Error: file name empty.. 3，程序利用多序列比对数据调用baseml/codeml命令对数据进行分析，生成out.BV文件。由于mcmctree调用baseml/codeml进行计算的参数设置可能不太好（特别时对蛋白序列进行计算时），推荐自己修该软件自动生成的baseml/codeml配置文件，然后再手动运行baseml/codeml命令，再整合其结果文件为out.BV文件。
   cleandata = 0         *设置是否移除不明确字符（N、？、W、R和Y等）或含以后gap的列后再进行数据分析：0，不需要，但在序列两两比较的时候，还是会去除后进行比较；1，需要。 
       clock = 2         *设置分子种方法：1，global clock方法，表示所有分支进化速率一致；2，independent rates方法，各分支的进化速率独立且进化速率的对数log(r)符合正态分布; 3，correlated rates方法，和方法2类似，但是log(r)的方差和时间t相关。
*    TipDate = 1 100     *当外部节点由取样时间时使用该参数进行设置，同时该参数也设置了时间单位。具体数据示例请见examples/TipData文件夹。
     RootAge = '<1.0'    *设置root节点的分歧时间，一般设置一个最大值。
 
***位点替换模型参数：**
       model = 4         *设置碱基替换模型：0，JC69；1，K80；2，F81；3，F84；4，HKY85。当设置usedata = 1时，不能使用其它碱基替换模型。若需要选择其它碱基替换模型，则考虑使用usedata = 3参数，就可以设置model参数值为其它的碱基替换模型。此时，程序会将数据进行分割，例如，ndata = 3，则程序分别对3个数据进行分析，生成baseml的输入文件，并调用baseml进行分析（若数据量较大，分析较慢，推荐按ctrl + c 强行终止，程序则终止baseml的运行，进行下一个数据的输入文件生成并运行下个数据的baseml命令，再按 ctrl + c 强行终止，这样可以让程序生成3各数据的输入文件和baseml配置文件，然后手动并行化运行，加快时间）；每个分析结果中会得到rst2结果文件；可以在此处手动修改tmp0001.ctl配置文件，例如，修改其model参数值，推荐修改outfile参数值，从而可以手动并行化运行baseml命令；然后将所有数据的rst2结果使用cat命令合并成文件in.BV，用于下一步设置参数usedata = 2，进行MCMC分析。
       alpha = 0.5       *核酸序列中不同位点，其进化速率不一致，其变异速率服从GAMMA分布。一般设置GAMMA分布的alpha值为0.5。若该参数值设置为0，则表示所有位点的进化速率一致。此外，当userdata = 2时，alpha、ncatG、alpha_gamma、kappa_gamma这4个GAMMA参数无效。因为userdata = 2时，不会利用到多序列比对的数据。
       ncatG = 5         *设置离散型GAMMA分布的categories值。
     BDparas = 1 1 0.1   *设置出生率、死亡率和取样比例。若输入有根树文件中的时间单位发生改变，则需要相应修改出生率和死亡率的值。例如，时间单位由100Myr变换为1Myr，则要设置成".01 .01 0.1"。
 kappa_gamma = 6 2       *设置kappa（转换/颠换比率）的GAMMA分布参数。
 alpha_gamma = 1 1       *设置GAMMA形状参数alpha的GAMMA分布参数。 

***进化速率参数：**
 rgene_gamma = 2 20 1    *设置序列中所所有位点平均[碱基/密码子/氨基酸]替换率的Dirichlet-GAMMA分布参数：alpha=2、beta=20、初始平均替换率为每100 million 年（取决于输入有根树文件中的时间单位）1个替换。若时间单位由100Myr变换为1Myr，则要设置成"2 2000 1"。总体上的平均进化速率为：2 / 20 = 0.1 个替换 / 每100Myr，即每个位点每年的替换数为 1e-9。
sigma2_gamma = 1 10 1    *设置所有位点进化速率取对数后方差（sigma的平方）的Dirichlet-GAMMA分布参数：alpha=1、beta=10、初始方差值为1。当clock参数值为1时，表示使用全局的进化速率，各分枝的进化速率没有差异，即方差为0，该参数无效；当clock参数值为2时，若修改了时间单位，该参数值不需要改变；当clock参数值为3时，若修改了时间单位，该参数值需要改变。
    finetune = 1: .1 .1 .1 .1 .1 .1    *冒号前的值设置是否自动进行finetune，一般设置成1，然程序自动进行优化分析；冒号后面设置各个参数的步进值：times, musigma2, rates, mixing, paras, FossilErr。由于有了自动设置，该参数不像以前版本那么重要了，可能以后会取消该参数。

***MCMC参数：**
       print = 1         *设置打印mcmc的取样信息：0，不打印mcmc结果；1，打印除了分支进化速率的其它信息（各内部节点分歧时间、平均进化速率、sigma2值）；2，打印所有信息。 
      burnin = 2000      *将前2000次迭代burnin后，再进行取样（即打印出该次迭代计算的结果信息，各内部节点分歧时间、平均进化速率、sigma2值和各分支进化速率等）。
    sampfreq = 10        *每10次迭代则取样一次。
     nsample = 20000     *当取样次数达到该次数时，则取样结束（程序也将运行结束）。至少要一万个

二.MAGA

CLOCKS -> Compute Timetree -> RelTime-ML
LOAD SEQUENCE DATA : alignment data (fasta format)
LOAD TREE FILE : tree data (nwk format,no branch)
Select Taxa : add outgroup data(click OK)
Select Branch : (click a branch then click OK)
CALIBRATE NOTE -> (click clock icon ) -> SELECT CALIBRATE TYPE(Normal Distribution) -> (click OK)
可以多线程运行，但是由于电脑的内存不足可能无法运行(8条序列，1101462 * 8个碱基，54GB内存)

也可以下载运行linux版本的maga。只不过下载的都是预编译的版本，需要系统内部特定版本的glibc ，为此还得专门去编译 glibc。

三.BEAST

BEAST1与BEAST2比较 - 高芳銮的博文

将多序列比对文件（fasta文件）转换成nexus文件（可以使用mega的另存为）
使用BEAST下面的子程序 beauti （可视化页面）将nex文件转换成xml文件。

# 导入beast 所需要的库文件，否则不能正常运行
git clone https://github.com/beagle-dev/beagle-lib.git
cd beagle-lib/
mkdir build
cd build/
cmake -DCMAKE_INSTALL_PREFIX=/home/bqxiao/program/beagle ..
make
make install
export LD_LIBRARY_PATH=/home/bqxiao/program/beagle/lib:$LD_LIBRARY_PATH

3.1 BEAST 1

BEAST 运行的本身命令比较简单：

nohup beast file.xml -threads N &
# beast（linux版） 本身是一个shell脚本，通过调用java来运行beast.jar程序，在beast脚本内部，可以设置给java虚拟机的内存大小。当选择核算替换模型和线程数时，应当考虑堆内存大小。

beauti 主要参数：
patitions : 在这个页面显示导入的数据（xml，nex）的分区信息。在这里可以设置分区之间是否使用相同的模型，先验树等信息
Taxa : 添加校准节点，可以在序列数据中设置分类单元子集(相同子集拖入 Included Taxa 内)，在后续可以添加校准时间
Tip ：为每一条序列添加注释（如：采集时间，采集地点等）
Sites ：进化模型选择
Clocks ：分子钟模型设置
Tree ：种群或者物种规模随时间的变化模型，主要考虑序列来自物种内还是物种间。
Priors ：为前面的所有模型添加先验信息，在这里可以设置前面的 Taxa 校准时间。
MCMC ：

Length of chain ：MCMC分析运行的总代数。可以通过以增加代数来提高ESS值（Effective Sample Size）使参数收敛。增加采样频率也可以提高ESS值，但是过于采样频繁也不会提高ESS值。
Log parameters every ：采样频率，即多少代生成一个样本。这个值的设置受Length of chain影响，最主要的一点需要保证最后的样本最少得有10,000个。
Echo state to screen ：将状态摘要输出到BEAST控制台窗口的频率。
最后点击生成xml文件

3.2 BEAST 2

BEAST 2中beast的使用方法和BEAST 1 中的beast使用方法一致。而子程序 beauti 有些许偏差。

四.r8s

r8s 分析速度非常快，对于一棵树只需要几秒钟可以得出结果。其中 pyr8s 是 r8s 的重写版本。其中 r8s 中的很多参数的（模型）算法只有部分在 pyr8s 得到的。

4.1 r8s v1.8.1

原本的r8s目前已经不支持了，仅有 r8s download | SourceForge.net 这个网站还有源码下载，但是依然有一些问题。
比如在编译 fortran 一些文件时，有些代码过于古老，特性已被删除，不能通过编译。
因此可以将make报错部分手动编译

# 原始make部分
gfortran -fopenmp -march=nocona -mtune=haswell -ftree-vectorize -fPIC -fstack-protector-strong -fno-plt -O2 -ffunction-sections -pipe -isystem /home/bqxiao/miniconda3/include -c -o tn.o tn.f
# 手动编译，保留已删除的特性： -std=legacy
# 其他标准还有 -std=f95 -std=f2003 -std=f2008
gfortran -std=legacy -fopenmp -march=nocona -mtune=haswell -ftree-vectorize -fPIC -fstack-protector-strong -fno-plt -O2 -ffunction-sections -pipe -isystem /home/bqxiao/miniconda3/include -c -o tn.o tn.f

另一部分编译报错于全局变量重定义 -> 只需要在重定义部分前置 extern 修饰即可。
使用 :

r8s -b -f file.nex > file.out
# -b 使用命令行而不是交互式

配置文件如下：

#NEXUS

begin trees;
tree PAUP_9 = (Marchantia:0.033817,(Lycopodium:0.040281,((Equisetum:0.048533,(Osmunda:0.033640,Asplenium:0.036526):0.000425):0.011806,((((Cycas:0.009460,Zamia:0.018847):0.005021,Ginkgo:0.014702):1.687e-86,((Pinus:0.021500,(Podocarpac:0.015649,Taxus:0.021081):0.006473):0.002448,(Ephedra:0.029965,(Welwitsch:0.011298,Gnetum:0.014165):0.006883):0.016663):0.006309):0.010855,((Nymphaea:0.016835,(((((Saururus:0.019902,Chloranth:0.020151):1.687e-86,((Araceae:0.020003,(Palmae:0.006005,Oryza:0.031555):0.002933):0.007654,Acorus:0.038488):0.007844):1.777e-83,(Calycanth:0.013524,Lauraceae:0.035902):0.004656):1.687e-86,((Magnolia:0.015119,Drimys:0.010172):0.005117,(Ranunculus:0.029027,((Nelumbo:0.006180,Platanus:0.002347):0.003958,(Buxaceae:0.013294,((Pisum:0.035675,(Fagus:0.009848,Carya:0.008236):0.001459):0.001994,(Ericaceae:0.019136,Solanaceae:0.041396):0.002619):1.687e-86):0.004803):1.687e-86):0.006457):0.002918):0.007348,Austrobail:0.019265):1.687e-86):1.687e-86,Amborella:0.019263):0.003527):0.021625):0.012469):0.019372);
End;

begin rates;

blformat nsites=952 lengths=persite;

collapse;

mrca LP marchantia pisum;
mrca ANGIO amborella pisum;

mrca VP lycopodium pisum;
mrca V1 Equisetum pisum;
mrca FERN Osmunda Asplenium;
mrca SP Ginkgo pisum;
mrca GYMNO Ginkgo Gnetum;
mrca CYC Cycas Zamia;
mrca G2 Taxus Podocarpac;
mrca GNET Ephedra Gnetum;

mrca CL Calycanth Lauraceae;
mrca ChS Chloranth Saururus;
mrca MON Acorus Oryza;
mrca PO Palmae Oryza;
mrca MAG Magnolia Drimys;
mrca EUD Ranunculus Carya;
mrca NEL Nelumbo Platanus;
mrca BX Buxaceae Carya;
mrca ERIC Ericaceae Solanaceae;
mrca FC Fagus Carya;

fixage taxon=LP age=450;

constrain taxon=CL min_age=110;
constrain taxon=ChS min_age=125;
constrain taxon=MON min_age=120;
constrain taxon=PO  min_age=85;
constrain taxon=MAG min_age=125;
constrain taxon=NEL min_age=110;
constrain taxon=BX  min_age=105;
constrain taxon=ERIC min_age=91;
constrain taxon=FC  min_age=96;

constrain taxon=VP  min_age=420;
constrain taxon=V1  min_age=390;
constrain taxon=FERN  min_age=255;
constrain taxon=GYMNO  min_age=310;
constrain taxon=CYC  min_age=220;
constrain taxon=GNET  min_age=115;
constrain taxon=G2  min_age=210;

constrain taxon=EUD max_age=125;
constrain taxon=SP max_age=330;

unfixage taxon=lp;

set num_time_guesses=5;
unfixage taxon=platanus;
constrain taxon=platanus maxage=90;

divtime method=LF algorithm=tn;
showage;
describe plot=chronogram;
describe plot=tree_description;

end;

4.2 pyr8s

安装非常简单

pip install git+https://github.com/iTaxoTools/pyr8s.git

4.2.1 python 内部使用

pyr8s 原文档使用方法有问题

import pyr8s.parse as p
a = p.parse('tests/legacy_1')
# 其中 parse.py 中根本没有 parse 函数

通过查看 run.py 可以发现，应当是以下代码：

import pyr8s.parse as p
a=p.from_file("tests/legacy_1")
result=a.run()

4.2.2 shell 环境中使用

pyr8s file.nex > result.txt

配置文件如下：

#NEXUS


begin trees;
tree Apiaceae = [&U] (((Angelica_sinensis:0.0263287636,((Apium_graveolens:0.0460894249,Heracleum_sosnowskyi:0.0322190162):0.0024624336,(Coriandrum_sativum:0.0283896098,Peucedanum_praeruptorum_Dunn:0.0307326278):0.0032486520):0.0020207251):0.0171794647,Oenanthe_sinensis:0.0537308346):0.0110179388,Daucus_carota:0.0638858678,cacao:0.6170394091);
End;


begin rates;

[ 
    nsites : 核苷酸的个数
    persite : 最大似然树
]
blformat nsites=1101462 lengths=persite;

collapse;


mrca DO Daucus_carota Oenanthe_sinensis;
mrca CH Coriandrum_sativum Heracleum_sosnowskyi;


constrain taxon=DO min_age=24 max_age=44;
constrain taxon=CH min_age=5 max_age=40;

set num_time_guesses=5;

divtime method=NP algorithm=PL;

[
    showage : 展示各个分支的节点的年龄
]
showage;

[
    describe plot=chronogram : 在终端画出树的形状
]
describe plot=chronogram;

[
    describe plot=tree_description : 在终端导出这棵树的nwk格式
]
describe plot=tree_description;

end;

版权声明：
作者：ht
链接：https://www.techfm.club/p/138920.html
来源：TechFM
文章版权归作者所有，未经允许请勿转载。

THE END

download GitHub 代码网站随机

二维码

注定要淋这场雨

< <上一篇

夏日之清凉

下一篇>>

搜索内容

时间树