功能性状的谱系保守性分析是架起进化历程与生态过程的重要桥梁，是很多高阶分析的基础。谱系保守性是很多基于谱系的预测模型的基本前提，比如生态位模型、达尔文归化假说、预适应假说。去除谱系保守性的影响，是跨物种进行性状相关性分析的前提。

为获得更好的代码阅读体验，请移步我的博客园博客。

示例数据

Species	H_N0	H_N2	H_N5	H_N10	H_N20	H_N50
Apocynum venetum	16.42	30.52	37.38	39.98	32.75	34.63
Cynanchum chinense	21.5	45.35	72.96	85.39	82.53	82
Limonium sinense	6.41	10.63	14.34	17.08	17.13	15.86
Chloris virgata	37.06	53.33	72.92	81.41	83.06	78.25
Setaria viridis	31.18	56.84	68.36	73.67	69.77	58.54
Suaeda glauca	24.54	38.52	41.47	45.49	46.29	43.86
Suaeda salsa	18.93	31.84	39.89	41.46	40.96	32.93
Phragmites australis	34.22	44.26	48.48	54.84	54.23	54.49
Xanthium strumarium	17.44	33.25	43	46.7	43.34	30.93
Xanthium sibiricum	16.34	23.03	33.4	43.43	44.03	37.32
Oenothera biennis	6.32	13.14	16.93	16.84	18.26	14.6
Chenopodium ficifolium	11.03	29.57	48.94	39.71	38.36	32.28
Salsola tragus	21.92	37.82	37.02	37.53	34.97	33.21
Polygonum sibiricum	13.92	22.59	27.43	23.57	21.05	22.59
Leymus mollis	30.89	39.74	44.83	38.86	38.86	42.49
Leonurus japonicus	1.98	6.45	10.12	11.72	13.31	8.73
Lactuca indica	10.22	21.82	29.01	28.99	26.33	23.14

 示例代码

#2021年4月1日 刘乐乐
#组会分享:功能性状的谱系信号与谱系独立比较的分析方法
#数据集来自祁鲁玉：黄河三角洲湿地植物不同物种对氮沉降响应实验

#V.PhyloMaker包的安装需要借助devtools::install_github,如下
library("devtools") #提供安装V.PhyloMaker包和plantlist物种名录所需要的工具
install_github("jinyizju/V.PhyloMaker")#安装V.PhyloMaker
install_github("helixcn/plantlist")#安装plantlist

####载入需要的软件包####
library("picante") 
library("V.PhyloMaker") 
library("plantlist") #物种名录，匹配分类地位（界门纲目科属种）

#载入数据集
trait <- read.table("data.txt", header = TRUE, sep = "\t", row.names = 1) #第一行为性状名称，第一列为物种名；务必保证没有多余空格

####建立谱系树####
#搜索物种的属和科
species_list <- rownames(trait)
sp_lis <- subset(TPL(species_list),select = c("YOUR_SEARCH", "POSSIBLE_GENUS", "FAMILY"))
colnames(sp_lis) <- c("species", "genus", "family")
#建树
sp_tree <- phylo.maker(sp.lis=sp_lis) #需要运行一小段时间（几分钟），需要耐心等待
sp_tr <- sp_tree$scenario.3
plot(sp_tr)
write.tree(sp_tr, "sp_tree.newick") #写入文件存储

####谱系保守性分析####
row.names(trait) <- gsub(" ", "_", rownames(trait)) #物种名与谱系树的label保持一致
#单个性状 phylosignal\Kcalc
trait1 <- trait[["H_N0"]]
names(trait1) <- rownames(trait)
phylosignal(trait1,sp_tr)
Kcalc(trait1, sp_tr)
#多个性状multiPhylosignal\Kcalc
multiPhylosignal(trait, sp_tr)
apply(trait,2,Kcalc, sp_tr)
#校对排序后，顺序一致，无需进行此步骤
data_matched <- match.phylo.data(phy = sp_tr, data = trait)
multiPhylosignal(data_matched$data, data_matched$phy) #结果同上
#可视化
trait_full <- trait
trait_full$sp <- row.names(trait_full)
color.plot.phylo(sp_tr,df = trait_full, trait = "H_N0", taxa.names = "sp", num.breaks = 3)
#trait_full$repsonse <- (trait_full$H_N0 + trait_full$H_N50)/(trait_full$H_N0)
#multiPhylosignal(trait_full[,-7], sp_tr)
#color.plot.phylo(sp_tr,df = trait_full, trait = "repsonse", taxa.names = "sp", num.breaks = 3)

####谱系独立比较###
#谱系独立的相关性分析
trait1_pic <- pic(trait1,sp_tr) #对单个性状计算pic值
trait_pic <- apply(trait, 2, pic, sp_tr)#对所有性状计算pic值
cor.test(trait[[1]],trait[[5]])
cor.test(trait_pic[,1],trait_pic[,5])
par(mfrow=c(1,2))
corrplot::corrplot(cor(trait), type = "upper")
corrplot::corrplot(cor(trait_pic), type = "upper")

 R语言及R studio使用小技巧

• 查看数据结构 str、class、head

• 查看帮助文件?、help、F1

• 运行例子 example

• 查看原始文献 citation  

• 查看源代码 直接输入函数名

• 需要学习的软件包 ade4、ape、picante、vegan

• 快捷键 ctrl+Enter

PPT及代码、数据文件可在百度网盘下载
链接：https://pan.baidu.com/s/1RUfmZAMkF1b9vuM3EwA7DQ
提取码：phyl

data.txt

script.R

谱系信号分析_LeleLiu.pptx