使用pvclust结合Bray-Curtis相似度计算p值的问题

问题根源 #

你当前代码的核心问题是聚类对象和bootstrap抽样方式不匹配：

• pvclust默认对矩阵的列（变量）进行聚类，而你需要聚类的是行（样本）
• 默认的method.boot = "centroid"是针对变量聚类的抽样策略，不适合样本聚类的稳定性评估，导致p值计算异常（全为0）

正确的pvclust用法 #

要实现样本的Bray-Curtis聚类及bootstrap p值计算，需要调整两个关键参数：cluster="rows"指定聚类行（样本），method.boot="sample"指定对样本进行bootstrap抽样：

# 加载依赖包
library(vegan)
library(pvclust)

# 加载示例数据
data(dune)

# 运行pvclust：聚类样本，Bray-Curtis距离，Ward聚类，样本bootstrap抽样
pv_result <- pvclust(dune,
                     cluster = "rows",          # 指定聚类行（样本）
                     method.hclust = "ward.D",
                     method.dist = function(x) vegan::vegdist(x, "bray"),
                     method.boot = "sample",    # 对样本进行bootstrap抽样
                     n = 1000,
                     parallel = TRUE)

# 可视化结果，标注p值
plot(pv_result)
pvrect(pv_result, alpha = 0.95)  # 高亮显著稳定的聚类节点

替代方案：手动结合hclust与boot包 #

如果pvclust的自定义距离仍有问题，可以用vegan计算距离，hclust构建聚类树，再用boot包手动实现bootstrap节点稳定性检验：

library(vegan)
library(boot)
library(fpc)

# 定义bootstrap函数：返回聚类树的节点结构
boot_clust <- function(data, indices) {
  # 抽样样本
  sampled_data <- data[indices, ]
  # 计算Bray-Curtis距离
  dist_mat <- vegdist(sampled_data, "bray")
  # 构建聚类树
  hc <- hclust(dist_mat, method = "ward.D")
  # 返回聚类树的合并矩阵
  return(hc$merge)
}

# 用clusterboot简化节点支持度计算
cb_result <- clusterboot(dune, clustermethod = hclustCBI,
                         method = "ward.D",
                         distance = vegdist,
                         distargs = list(method = "bray"),
                         B = 1000)

# 可视化聚类树并标注支持度
hc <- hclust(vegdist(dune, "bray"), method = "ward.D")
plot(hc)
nodelabels(cb_result$bootmean, cex = 0.8)

另一个替代方案：使用phyloseq包（针对微生物组数据） #

如果你的数据是微生物组物种丰度表，phyloseq包提供了更便捷的样本聚类及bootstrap支持：

library(phyloseq)
library(vegan)

# 将dune转换为phyloseq对象
ps <- phyloseq(otu_table(dune, taxa_are_rows = FALSE))

# 计算Bray-Curtis距离并构建聚类树
dist_bray <- distance(ps, "bray")
hc <- hclust(dist_bray, method = "ward.D")

# 用phyloseq的bootstrap函数计算节点支持度
boot_support <- bootstrap_trees(ps, FUN = function(x) {
  hclust(distance(x, "bray"), method = "ward.D")
}, R = 1000)

# 可视化带支持度的聚类树
plot_phylo_tree(boot_support, type = "dendrogram", show.node.label = TRUE)

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