如何用vegan包基于计数表与元数据计算并绘制Shannon/Simpson指数

问题背景 #

需要基于提供的原始物种计数表和元数据完成两项任务：

1. 计算Shannon多样性指数与Simpson多样性指数
2. 按元数据中的yield和day分组绘制多样性指数可视化图表

用户此前尝试转置数据后运行代码，得到的是物种层面的多样性结果（而非样本层面），不符合预期需求。

原始数据 #

物种计数表 #

Datasets    R21P2_BS_sortme_non_rRNA_diamond    R22P1_BS_sortme_non_rRNA_diamond    R22P4_BS_sortme_non_rRNA_diamond    R22P6_BS_sortme_non_rRNA_diamond    R23P4_BS_sortme_non_rRNA_diamond    R23P5_BS_sortme_non_rRNA_diamond    R22P2_BS_sortme_non_rRNA_diamond    R22P3_BS_non_rRNA_diamond
Candidatus Koribacter   0   0   0   0   0   0   0   0
Candidatus Koribacter versatilis    6450.7397   3675.8171   2588.1077   2657.2302   2345.8489   0   6476.0293   6215.4624
Candidatus Sulfotelmatomonas    0   0   0   0   0   0   0   0
Candidatus Sulfotelmatomonas gaucii 0   0   0   0   0   2988    0   0
Edaphobacter    2452.491    2315.793    0   0   2334.1106   955 3006.5918   2504.095
Edaphobacter modestus   0   0   0   0   0   2791    0   0
Occallatibacter 0   0   0   0   0   0   0   0
Occallatibacter savannae    3342.1418   0   0   0   0   0   2275.6936   0
unclassified Acidobacteriales   713.793 1381.9713   371.48083   372.63043   413.5484    0   1512.612    1004.95465
Acidobacteriales bacterium  1989.0435   2009.239    2346.533    2188.6824   2098.4421   0   2224.878    2977.275
Bryobacterales  92.96951    63.00996    77.73442    82.42516    65.914925   47  74.10613    63.01696
Bryobacteraceae 632.5847    539.4786    575.83264   574.2777    421.6749    344 999.37415   679.9198
Paludibaculum   0   0   0   0   0   0   0   0
Paludibaculum fermentans    5697.7427   3003.9468   3720.4888   3304.1206   2315.1487   2032    4642.855    4024.7937
Solibacteraceae 43.124413   14.159542   25.562664   24.5313 26.185379   26  37.688263   28.744577

元数据文件 #

id,day,yield
R21P2_BS_sortme_non_rRNA_diamond,0,1
R22P1_BS_sortme_non_rRNA_diamond,30,1
R22P4_BS_sortme_non_rRNA_diamond,30,0
R22P6_BS_sortme_non_rRNA_diamond,30,0
R23P4_BS_sortme_non_rRNA_diamond,60,0
R23P5_BS_sortme_non_rRNA_diamond,60,0
R22P2_BS_sortme_non_rRNA_diamond,30,1
R22P3_BS_non_rRNA_diamond,30,1

用户错误尝试说明 #

错误代码 #

veganData  <- read.csv("./transposed-bs-allseason-230922-taxonomy-phylum.csv", header= TRUE,sep ='\t')
diversity(veganData[, -(1:4)],"simpson",MARGIN = 2)

错误输出 #

candidate.division.NC10        candidate.division.WOR.3
0.9525816                       0.8715260
Candidatus.Aminicenantes        Candidatus.Eisenbacteria
0.9551888                       0.9532059
Candidatus.Rokubacteria                 Calditrichaeota
1.0000000                       0.4999618
Elusimicrobia                   Bacteroidetes
0.9380819                       0.9570427
Ignavibacteriae candidate.division.Zixibacteria
0.9466634                       0.7953193

错误原因 #

MARGIN=2指定按列计算多样性，得到的是每个物种的多样性值，而我们需要按样本行计算每个样本的多样性指数。

完整解决方案 #

步骤1：安装并加载所需包 #

# 首次运行需安装包
install.packages(c("vegan", "ggplot2", "dplyr", "tidyr"))

# 加载包
library(vegan)
library(ggplot2)
library(dplyr)
library(tidyr)

步骤2：读取并预处理数据 #

# 读取物种计数表，将第一列设为行名（物种名）
species_counts <- read.table("物种计数表文件名.txt", header = TRUE, row.names = 1, sep = "\t")

# 转置数据：让样本作为行，物种作为列（vegan要求样本在行、物种在列）
species_counts_t <- t(species_counts)

# 读取元数据
metadata <- read.csv("元数据文件名.csv", header = TRUE)

# 确保样本顺序与元数据一致，避免后续合并出错
species_counts_t <- species_counts_t[metadata$id, ]

步骤3：计算多样性指数 #

# 计算Shannon指数（默认MARGIN=1，按行计算样本多样性）
shannon_index <- diversity(species_counts_t, index = "shannon")

# 计算Simpson指数
simpson_index <- diversity(species_counts_t, index = "simpson")

# 将指数合并到元数据中
metadata$shannon <- shannon_index
metadata$simpson <- simpson_index

步骤4：分组可视化多样性指数 #

方法1：箱线图（展示样本分布）

# 将数据转为长格式，适配ggplot绘图
metadata_long <- metadata %>%
  pivot_longer(cols = c(shannon, simpson), 
               names_to = "指数类型", 
               values_to = "指数值")

# 绘制分组箱线图
ggplot(metadata_long, aes(x = factor(day), y = 指数值, fill = factor(yield))) +
  geom_boxplot(position = position_dodge(0.8)) +
  geom_jitter(width = 0.2, position = position_dodge(0.8), alpha = 0.7) +
  facet_wrap(~指数类型, scales = "free_y") +
  labs(x = "天数(day)", y = "多样性指数值", fill = "产量(yield)") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

方法2：带误差棒的柱状图（展示分组均值与波动）

# 按day和yield分组计算均值与标准差
summary_data <- metadata_long %>%
  group_by(day, yield, 指数类型) %>%
  summarise(均值 = mean(指数值),
           标准差 = sd(指数值))

# 绘制柱状图
ggplot(summary_data, aes(x = factor(day), y = 均值, fill = factor(yield))) +
  geom_col(position = position_dodge()) +
  geom_errorbar(aes(ymin = 均值 - 标准差, ymax = 均值 + 标准差),
                width = 0.2, position = position_dodge(0.9)) +
  facet_wrap(~指数类型, scales = "free_y") +
  labs(x = "天数(day)", y = "平均多样性指数", fill = "产量(yield)") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

内容的提问来源于stack exchange，提问作者sumitra sivaprakasam