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工具可以获得更准确的系统发育关系。以下是完整流程：

1. 将分类学信息保存为文本文件：

writeLines(taxa_table, "taxa_list.txt")

2. 使用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 准备性状数据

假设我们有以下性状数据（通常来自实验测量或文献）：

将数据导入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个物种的中等规模数据集时，这种方法既保证了科学性，又大大节省了时间。