别再只画箱图了！用R的ggpubr玩转α多样性差异分析：Wilcoxon检验与高级可视化技巧

超越基础箱图：用ggpubr打造α多样性分析的科研级可视化

在微生物组研究中，α多样性分析是评估样本内物种丰富度和均匀度的关键步骤。虽然基础的箱线图能够展示不同组间的多样性差异，但当我们需要将结果呈现给学术同行或整合到研究报告中时，这些简单的图表往往显得过于粗糙。本文将带你深入探索如何利用R语言中的ggpubr包，从统计检验方法选择到可视化细节优化，打造出既科学严谨又美观专业的α多样性分析图表。

1. 数据准备与环境配置

在开始高级可视化之前，我们需要确保数据结构和计算方法的准确性。α多样性指数计算通常使用vegan包，而可视化则主要依赖ggplot2和ggpubr的组合。

# 安装必要包（如未安装） if (!require("vegan")) install.packages("vegan") if (!require("ggplot2")) install.packages("ggplot2") if (!require("ggpubr")) install.packages("ggpubr") # 加载包 library(vegan) library(ggplot2) library(ggpubr)

数据导入和预处理是分析的基础。以下是一个典型的数据处理流程：

# 读取OTU表格（假设为制表符分隔的文本文件） otu_table <- read.table("otu_table.txt", header=TRUE, row.names=1, sep="\t") # 转置矩阵使样本为行，物种为列 otu_t <- t(otu_table) # 计算香农多样性指数 shannon_index <- diversity(otu_t, index="shannon") # 读取样本分组信息 metadata <- read.table("metadata.txt", header=TRUE, row.names=1) # 合并多样性指数和分组信息 analysis_data <- data.frame( Sample=names(shannon_index), Shannon=shannon_index, Group=metadata[names(shannon_index), "Group"] )

注意：实际应用中，你可能需要根据数据特点调整参数，如文件分隔符(sep)、行名列名等。

2. 统计检验方法的选择与实施

组间差异的统计检验是α多样性分析的核心。虽然Wilcoxon秩和检验(也称为Mann-Whitney U检验)是微生物组研究中的常用方法，但选择适当的统计方法需要考虑数据特性。

2.1 检验方法比较

对于微生物组数据，由于多样性指数往往不符合正态分布，Wilcoxon检验通常是更稳妥的选择。

2.2 实施组间比较

ggpubr包的stat_compare_means()函数简化了统计检验的实施过程：

# 定义需要比较的组别 comparisons <- list(c("Healthy", "Disease"), c("Healthy", "Treatment"), c("Disease", "Treatment")) # 基础箱线图添加统计检验 ggplot(analysis_data, aes(x=Group, y=Shannon, fill=Group)) + geom_boxplot() + stat_compare_means(comparisons=comparisons, method="wilcox.test", label="p.signif")

参数说明：

• comparisons：指定需要两两比较的组别列表
• method：统计检验方法，常用"wilcox.test"或"t.test"
• label：显示格式，"p.signif"显示显著性符号(*)，"p.format"显示具体p值

提示：当比较组别较多时，可考虑先进行整体检验(Kruskal-Wallis或ANOVA)，再对显著结果进行两两比较，以控制多重检验带来的假阳性。

3. 高级可视化技巧

基础箱线图虽然直观，但缺乏专业图表应有的精细度和信息量。下面介绍几种提升图表质量的实用技巧。

3.1 图层叠加与美化

将箱线图与散点图结合，既能展示数据分布，又能呈现个体样本信息：

# 定义颜色方案 group_colors <- c("Healthy"="#4E79A7", "Disease"="#E15759", "Treatment"="#59A14F") # 构建高级图表 advanced_plot <- ggplot(analysis_data, aes(x=Group, y=Shannon)) + # 半透明箱线图 geom_boxplot(aes(fill=Group), width=0.6, alpha=0.7, outlier.shape=NA) + # 添加散点并添加随机抖动避免重叠 geom_jitter(aes(color=Group), width=0.2, size=3, alpha=0.6) + # 自定义颜色 scale_fill_manual(values=group_colors) + scale_color_manual(values=group_colors) + # 添加统计检验 stat_compare_means(comparisons=comparisons, method="wilcox.test", label="p.signif", tip.length=0.01) + # 主题美化 theme_minimal() + theme( panel.border=element_rect(color="black", fill=NA, size=0.8), axis.text=element_text(size=12), axis.title=element_text(size=14), legend.position="none" ) + labs(x=NULL, y="Shannon Diversity Index") print(advanced_plot)

3.2 主题定制与细节调整

ggplot2的主题系统允许我们对图表每个元素进行精细控制：

# 自定义主题函数 custom_theme <- function() { theme( panel.background=element_blank(), panel.grid.major=element_line(color="gray90", size=0.3), panel.grid.minor=element_blank(), axis.line=element_line(color="black", size=0.5), axis.text=element_text(color="black", size=11), axis.title=element_text(size=13, face="bold"), legend.title=element_text(size=12), legend.text=element_text(size=11), plot.title=element_text(size=16, hjust=0.5, face="bold"), plot.margin=unit(c(1,1,1,1), "cm") ) } # 应用自定义主题 advanced_plot + custom_theme()

3.3 多指标可视化

当需要同时展示多个α多样性指数(如Shannon、Simpson、Chao1)时，可以采用分面或组合图表：

# 计算多个多样性指数 analysis_data$Simpson <- diversity(otu_t, index="simpson") analysis_data$Chao1 <- estimateR(otu_t)["S.chao1",] # 转换数据格式便于绘图 library(tidyr) diversity_long <- gather(analysis_data, key="Index", value="Value", Shannon, Simpson, Chao1) # 分面绘制多个指数 ggplot(diversity_long, aes(x=Group, y=Value, fill=Group)) + geom_boxplot(alpha=0.7) + geom_jitter(width=0.2, size=2, alpha=0.5) + facet_wrap(~Index, scales="free_y", nrow=1) + scale_fill_manual(values=group_colors) + custom_theme() + stat_compare_means(comparisons=comparisons, method="wilcox.test", label="p.signif")

4. 出版级图表导出

科研图表通常需要满足出版物的高分辨率要求，并支持矢量格式以便后期编辑。

4.1 导出参数优化

# 导出为PDF（矢量图） ggsave("alpha_diversity.pdf", plot=advanced_plot, width=8, height=6, units="in", device=cairo_pdf) # 确保字体嵌入 # 导出为TIFF（位图，适合期刊投稿） ggsave("alpha_diversity.tiff", plot=advanced_plot, width=180, height=135, units="mm", dpi=600, compression="lzw")

4.2 图表元素调整清单

在最终导出前，检查以下元素是否优化：

• 字体大小：轴标签通常10-12pt，标题14-16pt
• 图例位置：避免遮挡数据，通常右上或底部
• 颜色对比：确保在黑白打印时仍可区分
• 比例协调：宽度与高度比例通常1:1到2:1之间
• 边距设置：避免图表元素过于靠近边缘

4.3 自动化报告生成

对于需要频繁更新的分析，可以结合R Markdown创建自动化报告：

```{r setup, include=FALSE} knitr::opts_chunk$set(echo=FALSE, warning=FALSE, message=FALSE) ``` # α多样性分析报告 ## 样本信息 - 总样本数: `r nrow(analysis_data)` - 组别: `r paste(unique(analysis_data$Group), collapse=", ")` ## 多样性分析结果 ```{r diversity-plot, fig.width=8, fig.height=6} # 插入前面创建的图表 print(advanced_plot) ``` 统计检验结果显著性： - Healthy vs Disease: `r format(wilcox.test(Shannon ~ Group, data=subset(analysis_data, Group %in% c("Healthy","Disease")))$p.value, scientific=TRUE)` - Healthy vs Treatment: `r format(wilcox.test(Shannon ~ Group, data=subset(analysis_data, Group %in% c("Healthy","Treatment")))$p.value, scientific=TRUE)`

在实际项目中，我发现将颜色方案与分组变量关联存储在数据框中特别有用，这样在修改分组名称或添加新组时，颜色映射会自动保持一致。另一个实用技巧是在geom_jitter()中使用position=position_jitterdodge()参数，当结合分组箱线图时，它可以更精确地控制散点位置。