R语言生态学：进化树推断物种分化历史：分类单元数与时间关系、支系图可视化

原创

拓端

发布于 2023-01-30 19:16:03

6060

发布于 2023-01-30 19:16:03

文章被收录于专栏：拓端tecdat

全文链接：http://tecdat.cn/?p=31434

原文出处：拓端数据部落公众号

我们围绕进化树技术进行一些咨询，分析生物类群在时间上的多样性是如何变化的。我们将用到分类单元数-时间图（Lineages-through-time plot）,该图可以用来描述物种多样化的总体趋势。

数据

3500trees.nexus是nexus格式的文件，里面有3500棵树。

besttree.nexus也是nexus格式的文件，里面有1颗树，是从3500颗树中筛选出来的一致树。

各支系图示

这棵树总共有4大支系（Lineage），现在我需要做的分析都是需要分别做总的，以及4个支系的，也就是说同样的分析要做5次，针对5组不同的对象。

分析方法

方法主要是物种多样化速率(diversification rate) 相关的内容。

trees=read.nexus("3500trees.nexus")  
besttree=read.nexus("besttree.nexus")

1.mltt plot (multiple lineage through time)

分类单元数-时间图

lingeage的数值取log的，95%置信区间的ltt plot,中间黑色线的是besttree的，要显示出来。分别做总的，以及4个支系的，共5个图。

plot(trees)

,log='y')

besttree

# 95% ltt置信区间

ltt.ci<-function(tree.all){  
  ntip=length(tree.all[[1]]$tip.label)  
  ntree=length(tree.all)

2.gamma statistic

检验分化速率的变化趋势，看γ的值是正的还是负的。结果需要得到每组的γ值及P值。

mmaStat(besttree)

## [1] -3.693285

3. Monte Carlo constant rates test

检验样品不全是否对分化速率的结果有显著的影响，应该也是每组都要做的。

mc.out <- mcmc.pop
line(tree.hiv)  
plot(sk, l

4.对每个组做几个模型的检验，主要包括Pure-birth, birth-death, Yule 2-rate,density-dependent logistic,density-dependent exponential模型。

tdAICr
## --------------Model Summary----------------  
##  
## MODEL pureBirth  
##  
## Parameters:  r1  
##  
## LH 535.1086  
##  
## AIC -1068.217  
##  
## r1 0.1817879  
##  
## a -1068.217  
##  
##  
## --------------------------  
## MODEL bd  
##  
## Parameters:  r1, a  
##  
## LH 535.1086  
##  
## AIC -1066.217  
##  
## r1 0.1817879  
##  
## a 0  
##  
##  
## --------------------------  
## MODEL DDL  
##  
## Parameters:  r1, k  
##  
## LH 542.2213  
##  
## AIC -1080.443  
##  
## r1 0.2537928  
##  
## a -1080.443  
##  
## k 554  
##  
##  
## --------------------------  
## MODEL DDX  
##  
## Parameters:  r1, X  
##  
## LH 536.991  
##  
## AIC -1069.982  
##  
## r1 0.3098342  
##  
## a -1069.982  
##  
## x 0.1131752  
##  
##  
## --------------------------  
##  
## Best Constant Rate Model = pureBirth  AIC  -1068.217  
##  
## Best Rate Variable Model = DDL  AIC  -1080.443  
##  
## delta AICrc =  12.2254

##       model params np mtype       LH        r1 r2         a        xp   k  
## 1 pureBirth     r1  1    RC 535.1086 0.1817879 NA -1068.217        NA  NA  
## 2        bd  r1, a  2    RC 535.1086 0.1817879 NA     0.000        NA  NA  
## 3       DDL  r1, k  2    RV 542.2213 0.2537928 NA -1080.443        NA 554  
## 4       DDX  r1, X  2    RV 536.9910 0.3098342 NA -1069.982 0.113175