如何用R语言vegan包计算组内相异度？替代PRIMER方案

问题描述 #

我正在处理海洋调查的丰度数据，用vegan包做分析，已经能用simper()计算组间相异度：

resultados_N_simper <- simper(especies_pw_hellinger, cluster)

但我之前觉得组内相异度只能用PRIMER软件计算，想问这个说法对不对？能不能用R实现？找了很久没找到方法，不确定是做不到还是没找对路子。

附示例数据结构：

structure(list(grupo = c("Clust1", "Clust1", "Clust1", "Clust2", 
"Clust2", "Clust2"), C_cae = c(0, 0, 0, 0, 0, 0), G_arg = c(0, 
0, 0, 1261, 1581, 264), G_mac = c(0, 0, 0, 0, 0, 0), L_lep = c(0, 
0, 0, 0, 0, 0), M_lae = c(0, 0, 0, 0, 0, 0), M_adu = c(7, 13, 
44, 5, 1, 6), M_juv = c(129, 60, 104, 59, 38, 136), M_pou = c(4, 
0, 13, 166, 453, 194), M_mac = c(0, 0, 0, 0, 0, 0), N_aeq = c(0, 
0, 0, 0, 0, 0), P_ble = c(0, 0, 0, 0, 0, 0), T_sca = c(0, 0, 
0, 0, 0, 0), T_lus = c(11, 20, 6, 4, 15, 1), T_min = c(4, 0, 
0, 0, 0, 0)), class = c("tbl_df", "tbl", "data.frame"), row.names = c(NA, -6L))

解答 #

组内相异度完全可以用R语言实现，无需依赖PRIMER。以下是两种常用方法，均基于vegan包完成：

方法1：计算组内平均相异度 #

核心思路是按分组提取样本，计算组内两两样本的相异度后取平均值：

library(vegan)

# 加载示例数据
dat <- structure(list(grupo = c("Clust1", "Clust1", "Clust1", "Clust2", 
"Clust2", "Clust2"), C_cae = c(0, 0, 0, 0, 0, 0), G_arg = c(0, 
0, 0, 1261, 1581, 264), G_mac = c(0, 0, 0, 0, 0, 0), L_lep = c(0, 
0, 0, 0, 0, 0), M_lae = c(0, 0, 0, 0, 0, 0), M_adu = c(7, 13, 
44, 5, 1, 6), M_juv = c(129, 60, 104, 59, 38, 136), M_pou = c(4, 
0, 13, 166, 453, 194), M_mac = c(0, 0, 0, 0, 0, 0), N_aeq = c(0, 
0, 0, 0, 0, 0), P_ble = c(0, 0, 0, 0, 0, 0), T_sca = c(0, 0, 
0, 0, 0, 0), T_lus = c(11, 20, 6, 4, 15, 1), T_min = c(4, 0, 
0, 0, 0, 0)), class = c("tbl_df", "tbl", "data.frame"), row.names = c(NA, -6L))

# 分离物种数据与分组信息
species_dat <- dat[, -1]
# 计算相异度矩阵（这里用Bray-Curtis，可替换为Hellinger转换后的欧氏距离等）
dist_mat <- vegdist(species_dat, method = "bray")

# 定义函数：计算每个分组的平均组内相异度
calc_intra_dist <- function(group_vec, dist_matrix) {
  unique_groups <- unique(group_vec)
  sapply(unique_groups, function(g) {
    group_idx <- which(group_vec == g)
    if (length(group_idx) < 2) return(NA) # 样本数不足2无法计算
    intra_distances <- as.dist(dist_matrix[group_idx, group_idx])
    mean(intra_distances)
  })
}

# 运行计算
intra_mean_dist <- calc_intra_dist(dat$grupo, dist_mat)
print(intra_mean_dist)

方法2：类似SIMPER的组内物种贡献分析 #

如果需要像simper()那样，拆解每个物种对组内相异度的贡献，可以对每个分组单独调用simper()，通过虚拟分组实现组内两两样本的差异分析：

# 分组计算组内SIMPER分析
intra_simper_results <- lapply(unique(dat$grupo), function(g) {
  sub_species <- species_dat[dat$grupo == g, ]
  # 构造虚拟分组：每个样本为单独一组，让simper计算所有两两差异
  dummy_groups <- factor(seq_len(nrow(sub_species)))
  simper_out <- simper(sub_species, dummy_groups)
  
  # 整理物种贡献的平均值
  contrib_df <- do.call(rbind, simper_out)
  avg_contrib <- aggregate(contrib ~ species, data = contrib_df, mean)
  avg_contrib <- avg_contrib[order(-avg_contrib$contrib), ]
  
  list(group = g, species_contribution = avg_contrib)
})

# 查看结果
for (res in intra_simper_results) {
  cat("分组", res$group, "的物种贡献排名:\n")
  print(res$species_contribution)
  cat("\n")
}

注意事项 #

• 相异度方法可灵活替换：若使用Hellinger转换后的数据，建议用vegdist(hellinger_data, method = "euclidean")计算欧氏距离，这是Hellinger转换后的标准搭配。
• 若分组内样本量小于2，无法计算组内相异度，函数会返回NA，需提前筛选有效分组。

内容的提问来源于stack exchange，提问作者Juan Carlos