最近做了一個(gè)祖先狀態(tài)重建的工作，還是蠻有意思的。PS：需要說一下的是，RASP的祖先狀態(tài)重建主要是分布區(qū)重建，用RASP的模型去做性狀重建會(huì)出現(xiàn)非常奇怪的結(jié)果，顧換成R語言來做這個(gè)分析。
# R語言祖先狀態(tài)重建
# 作者：小黑
# 日期：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)]
可以看出來我們的系統(tǒng)樹大概長這個(gè)樣子
# 讀取數(shù)據(jù)文件
tri_tab <- read.table("pileus and substrate.txt")
row.names(tri_tab) <- tri_tab[,1]
tri_tab <- tri_tab[,-1]
[站外圖片上傳中...(image-8a4771-1695130589218)]
我們的特征表差不多長這樣，一定要名稱要和樹上的名稱對(duì)應(yīng)上。
# 賦值
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)) # 將三個(gè)圖放在一起顯示
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)
# 查看三個(gè)模型計(jì)算出來的AIC值，AIC值越小，表示模型越精確
# 我這里ER=84，ARD=122，SYM=100，所以模擬效果最好
[站外圖片上傳中...(image-9a35be-1695130589218)]
# 基于模型估計(jì)內(nèi)部節(jié)點(diǎn)狀態(tài)
fitER.V <- ace(pileus,tri.tree,model="ER",type="discrete")
# 可以指定method “ML”，“REML”，“pic”或“GLS”
fitER.V
fitER.V$lik.anc
# Stochastic character mapping，進(jìn)行特征映射：
# 提出可以從計(jì)算的聯(lián)合貝葉斯后驗(yàn)概率分布中，對(duì)單個(gè)節(jié)點(diǎn)進(jìn)行重復(fù)取樣，來重建進(jìn)化歷史。這個(gè)過程需要
#（1）計(jì)算Q轉(zhuǎn)移矩陣；通過進(jìn)化樹的枝長和狀態(tài)轉(zhuǎn)移偏好來獲得先驗(yàn)信息。
#（2）從聯(lián)合條件概率分布中對(duì)一系列節(jié)點(diǎn)狀態(tài)進(jìn)行取樣；
#（3）基于進(jìn)化樹分支，取樣的節(jié)點(diǎn)狀態(tài)模擬進(jìn)化歷史。
## 單次Stochastic mapping
mtree <- make.simmap(tri.tree,pileus,model="ER")
### 如果需要末端對(duì)齊的話就將edge.length設(shè)置為Null，否則不用
mtree$edge.length <- NULL
### 畫出來看看
plot.phylo(mtree,type="phylogram")
add.simmap.legend(colors=cols,prompt=FALSE,x=0.2,y=27,fsize=1)
## 重復(fù)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)
至此，分析和畫圖的整個(gè)過程就全部都畫完了
