R语言 plantlist 0.6.5 + ape 包：14个物种系统发育树构建与性状数据整合实战
在生态学和进化生物学研究中，系统发育树（Phylogenetic Tree）是理解物种间进化关系的重要工具。传统构建系统发育树需要基因序列数据，但对于许多研究场景，特别是植物学研究，我们往往只需要基于分类学关系构建一个合理的系统发育框架。本文将详细介绍如何利用R语言中的plantlist、ape和phylosignal等包，实现从物种名录到系统发育树构建，再到性状数据整合可视化的完整流程。
1. 环境准备与数据获取
1.1 安装必要的R包
在开始之前，确保已安装以下R包。如果尚未安装，可以使用以下命令安装：
install.packages(c("plantlist", "ape", "phylosignal", "phylobase", "readxl", "dplyr"))这些包的主要功能如下：
- plantlist：用于植物学名查询和验证
- ape：系统发育树分析和可视化
- phylosignal：系统发育信号分析
- phylobase：系统发育数据整合
- readxl：读取Excel文件
- dplyr：数据整理和操作
1.2 准备物种名录数据
假设我们需要为以下14个植物物种构建系统发育树：
species_list <- c( "Chimonanthus praecox", "Michelia figo", "Cycas revoluta", "Ginkgo biloba", "Litsea rotundifolia", "Bombax ceiba", "Delonix regia", "Ficus microcarpa", "Mangifera indica", "Dimocarpus longan", "Triadica sebifera", "Liquidambar formosana", "Alstonia scholaris", "Camphora officinarum" )提示：实际研究中，物种名录通常存储在Excel或CSV文件中。可以使用
readxl::read_excel()或read.csv()导入。
2. 物种分类学信息查询
2.1 使用plantlist包验证学名
plantlist包的TPL()函数可以验证植物学名并获取分类学信息：
library(plantlist) taxa_check <- TPL(species_list)查看前几行结果：
head(taxa_check)输出将包含每个物种的科、属、种信息，以及APG III分类系统中的科编号等。
2.2 转换分类学格式
为了后续构建系统发育树，我们需要将查询结果转换为"科/属/种"格式：
taxa_table <- taxa.table(taxa_check) print(taxa_table)输出示例：
[1] "Calycanthaceae/Chimonanthus/Chimonanthus_praecox" [2] "Magnoliaceae/Michelia/Michelia_figo" [3] "Cycadaceae/Cycas/Cycas_revoluta" ...3. 系统发育树构建
3.1 使用Phylomatic生成Newick格式树
虽然可以直接在R中完成所有步骤，但使用Phylomatic工具可以获得更准确的系统发育关系。以下是完整流程：
- 将分类学信息保存为文本文件：
writeLines(taxa_table, "taxa_list.txt")- 使用Phylomatic的在线服务或本地awk版本生成Newick格式树文件（保存为"tree.newick"）
3.2 在R中读取和可视化系统发育树
使用ape包读取和绘制系统发育树：
library(ape) phy_tree <- read.tree("tree.newick") # 或直接使用字符串 plot(phy_tree, cex = 0.8)基础绘图效果可能不够理想，我们可以进行美化：
plot(phy_tree, type = "fan", edge.width = 2, tip.color = "blue", cex = 0.8, main = "14种植物系统发育树")4. 性状数据整合与可视化
4.1 准备性状数据
假设我们有以下性状数据（通常来自实验测量或文献）：
| species | biomass | leaf_size | growth_rate |
|---|---|---|---|
| Chimonanthus praecox | 12.5 | 15.2 | 0.45 |
| Michelia figo | 8.7 | 12.1 | 0.38 |
| ... | ... | ... | ... |
将数据导入R并整理：
library(dplyr) trait_data <- read.csv("trait_data.csv") %>% mutate(species = gsub(" ", "_", species)) %>% group_by(species) %>% summarise(across(everything(), mean, na.rm = TRUE)) %>% column_to_rownames("species")4.2 性状数据与系统发育树整合
使用phylobase包将性状数据映射到系统发育树上：
library(phylobase) library(phylosignal) p4d <- phylo4d(phy_tree, trait_data)可视化整合结果：
barplot.phylo4d(p4d, tree.type = "phylo", tree.ladderize = TRUE, trait.bg.col = "lightblue", scale = FALSE)4.3 系统发育信号分析
phylosignal包可以量化性状的系统发育信号：
phylo_signal <- phyloSignal(p4d) print(phylo_signal)结果将显示每个性状的Blomberg's K和Pagel's λ等指标，帮助我们理解性状在系统发育树上的分布模式。
5. 高级分析与可视化技巧
5.1 系统发育主成分分析（Phylogenetic PCA）
library(phytools) phylo_pca <- phyl.pca(phy_tree, trait_data) biplot(phylo_pca)5.2 性状演化模型拟合
比较不同演化模型对性状数据的拟合优度：
models <- c("BM", "OU", "EB") aic_scores <- sapply(models, function(m) { fitContinuous(phy_tree, trait_data$biomass, model = m)$opt$aic })5.3 使用ggtree增强可视化
library(ggtree) library(ggplot2) ggtree(phy_tree) %<+% trait_data + geom_tiplab(aes(color = biomass), size = 3) + scale_color_gradient(low = "blue", high = "red") + theme_tree2()6. 完整可复现代码示例
以下是一个端到端的可执行脚本，整合了上述所有步骤：
# 加载所需包 library(plantlist) library(ape) library(phylobase) library(phylosignal) library(dplyr) library(ggtree) # 1. 物种名录 species <- c("Chimonanthus praecox", "Michelia figo", "Cycas revoluta", "Ginkgo biloba", "Litsea rotundifolia", "Bombax ceiba", "Delonix regia", "Ficus microcarpa", "Mangifera indica", "Dimocarpus longan", "Triadica sebifera", "Liquidambar formosana", "Alstonia scholaris", "Camphora officinarum") # 2. 验证学名并获取分类信息 taxa_check <- TPL(species) taxa_table <- taxa.table(taxa_check) # 3. 构建系统发育树 writeLines(taxa_table, "taxa_list.txt") # 此处需运行Phylomatic生成tree.newick phy_tree <- read.tree("tree.newick") # 4. 性状数据整合 trait_data <- data.frame( species = gsub(" ", "_", species), biomass = runif(14, 5, 15), leaf_size = runif(14, 10, 20) ) %>% column_to_rownames("species") # 5. 可视化 p4d <- phylo4d(phy_tree, trait_data) barplot.phylo4d(p4d) # 6. 系统发育信号分析 phyloSignal(p4d) # 7. ggtree高级可视化 ggtree(phy_tree) %<+% as.data.frame(trait_data) + geom_tiplab(aes(color = biomass), size = 3, offset = 0.1) + geom_tippoint(aes(size = leaf_size), shape = 16) + scale_color_gradient(low = "blue", high = "red") + theme_tree2()在实际项目中，这套流程已经帮助我快速分析了多个植物功能性状的系统发育模式。特别是在处理50-100个物种的中等规模数据集时，这种方法既保证了科学性，又大大节省了时间。
)
