最近做了一个祖先状态重建的工作，还是蛮有意思的。PS：需要说一下的是，RASP的祖先状态重建主要是分布区重建，用RASP的模型去做性状重建会出现非常奇怪的结果，顾换成R语言来做这个分析。
# R语言祖先状态重建
# 作者：小黑
# 日期：2023.9.19
# R version 4.3.0
library(geiger)
library(phytools)
library(ape)
library(vegan)
library(permute)
library(lattice)
library(picante)
# 读取树文件
setwd("D:/getOganelle/tricholomopsis/Phyllotopsidaceae/all/iqtree/phytools")
tri.tree <- read.nexus(file = "phyllotopsidaceae.nex")
## 画一下看看长啥样
plotTree(tri.tree,fsize=.6)
## 去掉不要的类群
tri.tree <- drop.tip(tri.tree,c("DRR003994", "DRR003993"))
[站外图片上传中...(image-208d98-1695130589217)]
可以看出来我们的系统树大概长这个样子
# 读取数据文件
tri_tab <- read.table("pileus and substrate.txt")
row.names(tri_tab) <- tri_tab[,1]
tri_tab <- tri_tab[,-1]
[站外图片上传中...(image-8a4771-1695130589218)]
我们的特征表差不多长这样，一定要名称要和树上的名称对应上。
# 赋值
pileus <- setNames(tri_tab[,1],rownames(tri_tab))
substrate <- setNames(tri_tab[,2],rownames(tri_tab))
然后就开始做分析了
# 拟合模型
# "ER" is an equal-rates model
# "ARD" is an all-rates-different model
# "SYM" is an symmertrical model
par(mfrow=c(1,3)) # 将三个图放在一起显示
fitER.pileus <- fitDiscrete(tri.tree,pileus,model="ER")
fitER.pileus
plot(fitER.pileus,signif=5)
title(main="Fitted 'ER' model in pileus", line=-10)
fitARD.pileus <- fitDiscrete(tri.tree,pileus,model="ARD")
fitARD.pileus
plot(fitARD.pileus,signif=5)
title(main="Fitted 'ARD' model in pileus", line=-10)
fitSYM.pileus <- fitDiscrete(tri.tree,pileus,model="SYM")
fitSYM.pileus
plot(fitSYM.pileus,signif=5)
title(main="Fitted 'SYM' model in pileus", line=-10)
# 查看三个模型计算出来的AIC值，AIC值越小，表示模型越精确
# 我这里ER=84，ARD=122，SYM=100，所以模拟效果最好
[站外图片上传中...(image-9a35be-1695130589218)]
# 基于模型估计内部节点状态
fitER.V <- ace(pileus,tri.tree,model="ER",type="discrete")
# 可以指定method “ML”，“REML”，“pic”或“GLS”
fitER.V
fitER.V$lik.anc
# Stochastic character mapping，进行特征映射：
# 提出可以从计算的联合贝叶斯后验概率分布中，对单个节点进行重复取样，来重建进化历史。这个过程需要
#（1）计算Q转移矩阵；通过进化树的枝长和状态转移偏好来获得先验信息。
#（2）从联合条件概率分布中对一系列节点状态进行取样；
#（3）基于进化树分支，取样的节点状态模拟进化历史。
## 单次Stochastic mapping
mtree <- make.simmap(tri.tree,pileus,model="ER")
### 如果需要末端对齐的话就将edge.length设置为Null，否则不用
mtree$edge.length <- NULL
### 画出来看看
plot.phylo(mtree,type="phylogram")
add.simmap.legend(colors=cols,prompt=FALSE,x=0.2,y=27,fsize=1)
## 重复100次或1000次
mtrees <- make.simmap(tri.tree,pileus,model="ER",nsim=1000)
mtrees
par(mfrow=c(10,10))
null <- sapply(mtrees,plot,colors=cols,lwd=1,ftype="off")
## 整合100次或1000次树
pd <- summary(mtrees,plot=F)
# 画图
par(mfrow=c(1,1))
cols<-setNames(palette()[1:length(unique(pileus))],sort(unique(pileus)))
plot(pd,ftype="i",pie=pd$ace,piecol=cols,cex=0.35, asp=0.05)
add.simmap.legend(colors=cols,prompt=TRUE,x=1.1*par()$usr[1.5],y=0.01*par()$usr[1],fsize=0.8)
至此，分析和画图的整个过程就全部都画完了
网友评论