3、从宏基因组数据到α多样性指数：R实战解析与可视化-拓冰建站

1. 宏基因组数据与α多样性基础概念

宏基因组测序技术让我们能够直接获取环境样本中所有微生物的基因信息。这种高通量测序方法产生的数据通常以RPKM（Reads Per Kilobase per Million mapped reads）或TPM（Transcripts Per Million）等标准化形式呈现，表示不同基因或物种的相对丰度。

α多样性是描述单个样本内微生物群落多样性的重要指标。就像我们评估一个森林的生态价值时，不仅会看树木数量，还会考察树种丰富度和分布均匀度一样，α多样性通过三个维度反映微生物群落的健康状况：

丰富度指数：简单统计样本中物种的数量，如同计算森林里不同树种的总数
均匀度指数：衡量各物种分布的均衡程度，好比评估森林中各种树木的数量是否接近
综合指数：同时考虑丰富度和均匀度的复合指标

在实际研究中，我们最常用的α多样性指数包括：

Shannon指数：对稀有物种敏感，值越大表示多样性越高
Simpson指数：更关注优势物种，值越小多样性越高（注意不同文献中的计算方式差异）
Sobs：实际观测到的物种数目（又称Richness）
Chao1：通过数学模型估计样本中可能存在的物种总数

理解这些指数的生物学意义非常重要。比如在肠道微生物研究中，Shannon指数降低往往与疾病状态相关；而在环境微生物调查中，Simpson指数能更好反映优势菌群的变化。

2. R环境准备与数据导入

2.1 必要R包安装与加载

进行α多样性分析前，需要确保安装了以下关键R包：

# 安装核心包 install.packages(c("vegan", "ggplot2", "dplyr", "tidyr")) # 加载包 library(vegan) # 多样性分析核心工具 library(ggplot2) # 高级可视化 library(dplyr) # 数据整理 library(tidyr) # 数据重塑

vegan包是生态学分析的瑞士军刀，提供了从多样性计算到排序分析的全套函数。ggplot2则能创建出版级质量的统计图形。如果遇到包安装问题，可以尝试先更新R基础环境：

# 更新所有已安装包 update.packages(ask = FALSE)

2.2 数据导入与预处理

典型的宏基因组丰度表是二维矩阵，行代表物种/基因，列代表样本。假设我们有一个TSV格式的RPKM丰度表：

# 读取丰度表（注意修改为实际路径） abundance_table <- read.delim("path/to/your/RPKM_table.txt", row.names = 1, # 第一列为行名 header = TRUE, # 保留列名 check.names = FALSE) # 保持特殊字符 # 查看数据结构 dim(abundance_table) head(abundance_table[, 1:5]) # 显示前5列的部分数据

处理常见数据问题：

零值过多：建议过滤低丰度物种，保留至少在20%样本中出现的物种

abundance_filtered <- abundance_table[rowSums(abundance_table > 0) >= 0.2*ncol(abundance_table), ]

标准化：虽然RPKM已是标准化指标，但仍建议检查样本间测序深度差异

# 计算样本总reads col_sums <- colSums(abundance_filtered) # 可视化测序深度差异 barplot(col_sums, las = 2, cex.names = 0.7)

2.3 分组信息整合

统计分析通常需要样本分组信息，假设有分组文件group.txt：

group_info <- read.delim("path/to/group.txt", header = TRUE) # 确保样本顺序一致 rownames(group_info) <- group_info$SampleID group_info <- group_info[colnames(abundance_filtered), ]

3. α多样性指数计算实战

3.1 核心指数计算

使用vegan包的diversity()函数可以轻松计算多种α多样性指数：

# 计算Shannon指数（以自然对数为底） shannon <- diversity(t(abundance_filtered), index = "shannon") # 计算Simpson指数（注意这是1-D，值越大多样性越高） simpson <- diversity(t(abundance_filtered), index = "simpson") # 计算观测物种数(Sobs) richness <- specnumber(t(abundance_filtered)) # 计算Pielou均匀度指数 pielou <- shannon / log(richness)

为什么需要转置矩阵？因为vegan包默认期望行是样本、列是物种。如果数据是样本为行，则不需要t()转换。

3.2 结果整合与导出

将计算结果整理为数据框便于后续分析：

alpha_df <- data.frame( SampleID = names(shannon), Shannon = shannon, Simpson = simpson, Richness = richness, Pielou = pielou ) # 合并分组信息 alpha_df <- merge(alpha_df, group_info, by = "SampleID") # 保存结果 write.csv(alpha_df, "alpha_diversity_results.csv", row.names = FALSE)

3.3 指数选择建议

不同研究问题适合不同的α多样性指数：

关注稀有物种：优先使用Shannon指数
关注优势物种：选择Simpson指数
简单物种计数：使用Sobs或Chao1
临床研究：推荐同时报告Shannon和Richness

在环境梯度研究中（如pH梯度），Shannon指数通常能更好反映连续变化；而在比较明显不同的生境时，Simpson指数可能更敏感。

4. 统计分析与可视化

4.1 组间差异检验

常用的非参数检验方法：

# Shapiro检验正态性 shapiro.test(alpha_df$Shannon) # 若符合正态分布且方差齐性（p>0.05） t.test(Shannon ~ Group, data = alpha_df) # 更常用的非参数检验 kruskal.test(Shannon ~ Group, data = alpha_df) # 多组比较 wilcox.test(Shannon ~ Group, data = alpha_df) # 两组比较 # 事后检验（多组比较时） pairwise.wilcox.test(alpha_df$Shannon, alpha_df$Group, p.adjust.method = "BH")

4.2 ggplot2可视化实战

创建出版级质量的箱线图：

library(ggplot2) ggplot(alpha_df, aes(x = Group, y = Shannon, fill = Group)) + geom_boxplot(width = 0.6, outlier.shape = NA) + geom_jitter(width = 0.2, size = 2, alpha = 0.6) + scale_fill_brewer(palette = "Set2") + labs(x = "Experimental Group", y = "Shannon Diversity Index", title = "Comparison of Microbial Alpha Diversity") + theme_classic(base_size = 14) + theme(legend.position = "none", plot.title = element_text(hjust = 0.5))

进阶技巧：

添加统计显著性标记：

library(ggpubr) ggplot(...) + stat_compare_means(method = "kruskal")

分面显示多个指数：

alpha_long <- pivot_longer(alpha_df, cols = c(Shannon, Simpson, Richness)) ggplot(alpha_long, aes(x = Group, y = value, fill = Group)) + geom_boxplot() + facet_wrap(~name, scales = "free_y", ncol = 1)

4.3 结果解读要点

在论文中报告α多样性结果时应注意：

明确说明使用的指数及其生物学意义
报告具体的统计方法和p值
箱线图中标明样本量和中位数
对显著差异给出生物学解释
注意坐标轴刻度的一致性以便比较

常见错误：

仅使用Richness而忽略均匀度指数
未说明是否进行过测序深度标准化
忽略组内变异而过度解读均值差异
未校正多重假设检验

5. 完整流程自动化与进阶技巧

5.1 封装为可复用函数

将整个流程封装成函数提高效率：

calculate_alpha <- function(abundance, group, indices = c("shannon", "simpson", "richness")) { require(vegan) require(dplyr) result <- list() if("shannon" %in% indices) { result$shannon <- diversity(t(abundance), index = "shannon") } if("simpson" %in% indices) { result$simpson <- diversity(t(abundance), index = "simpson") } if("richness" %in% indices) { result$richness <- specnumber(t(abundance)) } result_df <- bind_cols(result) %>% mutate(SampleID = names(result[[1]])) %>% left_join(group, by = "SampleID") return(result_df) }

5.2 并行计算加速

对于大数据集可使用并行计算：

library(parallel) cl <- makeCluster(4) # 使用4个核心 clusterExport(cl, c("abundance_filtered", "diversity")) parLapply(cl, 1:100, function(i) { # 重采样计算 subsample <- abundance_filtered[, sample(ncol(abundance_filtered), replace = TRUE)] diversity(t(subsample), "shannon") }) stopCluster(cl)

5.3 与其他工具的衔接

将结果导出为其他分析工具格式：

# 导出为QIIME2格式 write.table(alpha_df[, c("SampleID", "Shannon")], "alpha.txt", sep = "\t", quote = FALSE, row.names = FALSE) # 导出为LEfSe输入格式 lefse_input <- cbind(SampleID = rownames(t(abundance_filtered)), t(abundance_filtered)) write.table(lefse_input, "lefse_input.txt", sep = "\t", quote = FALSE)

5.4 常见问题排查

遇到错误时的检查清单：

数据方向是否正确（样本为列还是行）
是否有NA或Inf值
分组信息是否与样本完全匹配
包版本是否兼容（特别是vegan包的更新可能改变参数）
矩阵是否包含非数值列

调试技巧：

# 检查数据异常值 summary(colSums(abundance_table)) # 检查分组匹配 setdiff(colnames(abundance_table), group_info$SampleID) # 小规模测试 test <- abundance_table[1:50, 1:5] diversity(t(test), "shannon")