层次聚类结果的比较和评估及R操作
下文的示例数据和R代码的获取链接：
https://pan.baidu.com/s/17Vcjhf5n6IGiyODn1eloMA
示例数据的层次聚类分析
文件“otu_table.txt”为某16S测序所得的细菌OTU丰度表格，包含了15个环境样本（分为a、b、c三组，即对应于3种不同的环境，每组各5个样本）中细菌物种的丰度组成。
接下来对该示例数据进行层次聚类划分样本集合，观测细菌群落组成的整体相似性或差异。
##聚类分析#读入原始物种丰度数据
dat_otu <- read.delim('otu_table.txt', row.names = 1)
dat_otu <- t(dat_otu)
#计算样本距离，以 Bray-curtis 距离和弦距离为例
library(vegan)
dis_bray <- vegdist(dat_otu, method = 'bray')
dis_ch <- vegdist(decostand(dat_otu, 'normalize'), method = 'euc')
#使用 Bray-curtis 距离的 UPGMA 聚类
clust_average <- hclust(dis_bray, method = 'average')
clust_average
summary(clust_average)
#使用弦距离的 Ward 最小方差聚类
clust_ward <- hclust(dis_ch, method = 'ward.D')
clust_ward
summary(clust_ward)
#聚类树
par(mfrow = c(1, 2))
plot(clust_average, main = 'UPGMA聚类\nBray-curtis距离', sub = '', xlab = 'Sample', ylab = 'Height')
plot(clust_ward, main = 'Ward最小方差聚类\n弦距离', sub = '', xlab = 'Sample', ylab = 'Height')
Bray-curtis距离和弦距离（弦距离具有欧式距离属性）都是物种多度分析中常用的距离测度，并选择了两种方法进行聚类分析（以便后续比较说明）。
通过绘制聚类树，可从中直观地辨别样本分布及各样本中细菌群落组成的整体差异。summary()用于查看并提取聚类树更多的细节信息，例如各样本在聚类树中的排列顺序以及高度属性等。
由于不确定那种结果更为理想，因此接下来通过一些方法评估并选择最合适的聚类模式。
评估聚类结果表征原始距离的精度
以下内容简介几种常用的评估方法。
同表型相关
一个聚类树内两个对象之间同表型距离（Cophenetic distance）是两个对象在同一组分类水平内的距离。任意两个对象，在聚类树上从一个对象向上走，到达与另一个对象交汇点向下走，势必会到达第二个对象，交汇点所在的层次水平即是两个对象同表型距离。同表型矩阵是表达所有对象的同表型距离的矩阵。同表型相关（Cophenetic correlation）是原始的距离矩阵和同表型距离矩阵之间的Pearson相关系数。具有最高的同表型相关系数的聚类方法可视为原始矩阵最好的聚类模式。
注：同表型矩阵由原始距离矩阵推导而来，两组距离矩阵不独立，不能检验同表型相关的显著性。同表型相关强烈依赖独立于数据的聚类方法的选择。
##同表型相关，比较上述 UPGMA、ward.D 结果#cophenetic() 获取聚类树中两两对象间的同表型矩阵，详情 ?cophenetic
clust_average.coph <- cophenetic(clust_average)
clust_ward.coph <- cophenetic(clust_ward)
#同表型距离矩阵和原始距离矩阵的同表型相关系数（即同表型相关）可使用 cor() 计算，
pearson_average <- cor(dis_bray, clust_average.coph, method = 'pearson')
pearson_ward <- cor(dis_ch, clust_ward.coph, method = 'pearson')
Gower距离
Gower距离等于原始距离与同表型距离之间差值的平方和，同样可用于比较聚类结果。一般来讲，具有最小Gower距离的聚类方法也可视为最具代表原始距离的聚类模式。
注：由于不同方法的评判标准不同，对于相同的数据来讲，Gower距离比较结果可能与上述同表型相关的比较结果并不完全一致。
#Gower 距离，比较上述 UPGMA、ward.D 结果
gow_average <- sum((dis_bray - clust_average.coph)^2)
gow_ward <- sum((dis_ch - clust_ward.coph)^2)
结合此处Gower距离评估结果与上文同表型相关评估结果，发现Ward最小方差聚类似乎并没有很好地反映原始距离测度。
Shepard图评估
还可借助Shepard图观测聚类树内对象的同表型聚类与原始对象距离的相互关系，以选择合适的聚类模式。同样的方法在非度量多维尺度（NMDS）分析中也有提及，在该文中通过使用Shepard图比较NMDS排序图内对象的距离与原始对象距离去评估NMDS结果。
#Shepard 图评估同表型距离与原始距离的相关性
par(mfrow = c(1, 2))
plot(dis_bray, clust_average.coph, pch = 19, col = 'blue', abline(0, 1, col = 'blue'), xlab = 'Bray-curtis距离', ylab = '同表型距离',
main = paste('UPGMA聚类\n同表型相关（Pearson）：', round(pearson_average, 3), '\nGower距离：', round(gow_average, 3)))
lines(lowess(dis_bray, clust_average.coph), col = 'red')
plot(dis_ch, clust_ward.coph, pch = 19, col = 'blue', abline(0, 1, col = 'blue'), xlab = '弦距离', ylab = '同表型距离',
main = paste('Ward最小方差聚类\n同表型相关（Pearson）：', round(pearson_ward, 3), '\nGower距离：', round(gow_ward, 3)))
lines(lowess(dis_ch, clust_ward.coph), col = 'red')
蓝色直线和红色曲线分别为常规线性拟合与lowess平滑拟合线，两种方法比较后，可以看到使用Bray-curtis距离的UPGMA方法是更理想的聚类模式。
当然并不是说使用弦距离不适合，更多原因为Ward最小方差聚类难以重现弦距离的原始结构。可以尝试使用UPGMA在弦距离的基础上进行聚类，比较其与Bray-curtis距离的UPGMA在聚类结果上的区别。
聚类树的生物学意义
当然，能充分体现生物学意义才是最关键的。
上述方法仅是建立在数学优度上而言的，评估的最优结果可能不足以描述所观测到的生物学现象；类似地，评估排名靠后的聚类模式可能反而在描述所观测到的生物学现象具有代表性。因此，这些方法通常作为辅助手段来使用，实际情况中结合目标问题综合看待。
下文的方法也是如此。
寻找可解读的聚类簇
为了解读和比较聚类的结果，通常需要寻找可解读的聚类簇，即寻找某个或某些划分水平来解读聚类的结果。划分水平的选择可以通过主观判断，也可以通过满足一些标准来确定，例如预先设定聚类簇的数量。
一般情况下，做完上一步的评估找到最合适的聚类结果后基本就满足需求了。毕竟在大多数条件下，我们很清楚自己的样本来源，进行聚类分析是想观测不同来源的样本是否能够被明显地划分为不同的集合，组间差异是否能够被区分，组内差异波动性是否较大，聚类树的分支结构是否像预期的那样具有某种分布或者趋势等。以本示例的数据为例，通过样本分组信息表已知15个样本分属于3个环境分组，找到的最适聚类模式恰好能将这15个样本明显地归为3大分支（UPGMA聚类树结果见上文），对应了其已知的3个分组，且相互之间未见模糊的重叠，那么就可以认为已经达到预期的分组目的。
但有的时候，我们还想继续往下探索聚类结果，例如期望在更高（或更低）层次继续划分类群，对聚类树进行裁剪，以寻找已知分组样本间的潜在层次关联；或者，手中拥有未知来源属性的一批数据，例如事先并不清楚分组信息的样本，期望通过聚类分析将它们合理地归类等；有时可能仅仅是想做个“测试”，评判已知归类分支的合理性等。再回到上文示例，虽然通过UPGMA方法恰好将15个样本聚为3大类，代表了3种不同环境来源，但若还想通过其它的一些方法做个评估，看这些方法是否也有效支持将聚类树划分为3大簇是最合适的。那么可以继续下文的内容。
融合水平值图
聚类树的融合水平值（fusion level value）是聚类树中两个分支融合处的相异性的数值，通过融合水平值的变化图有助于我们判断划分聚类簇的水平。
##融合水平值图，评判最优聚类簇数量
#样本数量
sample_num <- nrow(dat_otu)
#以 UPGMA 结果为例绘制融合水平值图
plot(clust_average$height, sample_num:2, type = 'S', col = 'blue', xlab = '节点高度（height）', ylab = '聚类簇数量（k）', main = 'UPGMA平均聚合聚类')
text(clust_average$height, sample_num:2, sample_num:2, col = 'red')
融合水平值图显示，当聚类簇的数量选择为3或4时，中间会产生一个明显的跳跃，因此建议最合适的聚类簇选择为3。
此时可以使用cutree()，根据预期的聚类簇数量对聚类树中的对象进行合理分类。融合水平值图显示建议保留3大簇（这与示例数据一开始已知的分组数量一致），因此默认将聚类树裁为3簇。
#裁剪聚类树，详情 ?cutreecut_average <- cutree(clust_average, k = 3)
#也可通过列联表比较多种聚类模式下的分类结果，例如将上文 Ward 聚类考虑进来
cut_ward <- cutree(clust_ward, k = 3)
table(cut_average, cut_ward)
两种聚类结果一致显示，15个样本的无监督模式下的分类与其已知的先验分组吻合。如果不考虑聚类方法在表征原始距离测度的精度上的问题，只论聚类簇是否正确识别了来自不同环境的群落，上文所得的两种聚类结果都是合适的。
轮廓宽度图
轮廓宽度（Silhouette width）是描述一个对象与其所属聚类簇归属程度的测度，是一个对象同同一组内其它对象的平均距离与该对象同最临近聚类簇内所有对象的平均距离的比较。
轮廓宽度值的范围(-1,1)，数值越大则聚类越好，负值则意味着该对象有可能被错分到当前聚类簇内。
继续使用轮廓宽度图帮助划分UPGMA结果的最优聚类簇数量。
#使用循环依次获取各数量下的聚类族，使用 cluster 包 silhouette() 获取各数量聚类簇下的轮廓宽度值
library(cluster)
silh_average <- rep(0, sample_num)
for (k in 2:(sample_num - 1)) silh_average[k] <- summary(silhouette(cutree(clust_average, k = k), dis_bray))$avg.width
#最佳（最大）轮廓宽度值
max(silh_average)
silh_average_max <- which.max(silh_average)
#绘制轮廓宽度图
plot(1:sample_num, silh_average, type = 'h', col = 'blue', xlab = '聚类簇数量（k）', ylab = '平均轮廓宽度',
main = paste('UPGMA-轮廓宽度图\n最优聚类簇数 k =', silh_average_max, '\n最佳轮廓宽度值：', round(max(silh_average), 3)))
points(silh_average_max, max(silh_average), type = 'h', col = 'red')
首先使用循环，依次将聚类树（UPGMA树）从划分为2个簇开始，直到第n-1个簇结束（n为对象总数），计算不同数量所得簇下的轮廓系数，并从中获取最大的值，其所对应的划分的聚类簇数量即为评估出的最适聚类簇数量。
结果清晰地展示了从划分的第2至第n-1个聚类簇下的轮廓宽度值，本示例UPGMA树最佳聚类簇数量为3。
同上文，此时可继续使用cutree()，根据预期的聚类簇数量对聚类树中的对象进行合理分类后观测分类结果，以进行最终确定。
距离矩阵和代表分组的二元矩阵的比较（以Mantel相关为例）
这种方法是通过计算原始距离与代表不同分类水平的二元矩阵（从聚类树计算获得）之间的相关性，选择最高的相关系数所对应的分类水平作为最优分组方案。Mantel相关是评估两个距离矩阵之间相关性的一种方法，相当于距离矩阵之间的Pearson相关系数。
注：代表聚类树的二元矩阵由原始距离矩阵推导而来，两组距离矩阵不独立，不能检验这种相关系数的显著性。
以下示例计算原始距离矩阵与聚类树中代表分组的二元矩阵的Mantel相关，通过相关系数的大小确定最佳聚类簇数量。
##原始距离矩阵和代表分组的二元矩阵的比较，以 UPGMA 结果为例，以 Mantel 相关为例
#使用循环依次获取各数量下的聚类族，使用 cluster 包 daisy() 获取各数量聚类簇下的二元矩阵
mantel_average <- rep(0, sample_num)
for (k in 2:(sample_num - 1)) {
clust_average_cut_k <- cutree(clust_average, k)
distgr_average <- daisy(as.data.frame(as.factor(clust_average_cut_k)), 'gower')
mantel_average[k] <- cor(dis_bray, distgr_average, method = 'pearson')
}
#最高相关系数
max(mantel_average)
mantel_average_max <- which.max(mantel_average)
#绘制轮廓宽度图
plot(1:sample_num, mantel_average, type = 'h', col = 'blue', xlab = '聚类簇数量（k）', ylab = '平均轮廓宽度',
main = paste('UPGMA-Mantel相关\n最优聚类簇数 k =', mantel_average_max, '\n最高Mantel相关系数：', round(max(mantel_average), 3)))
points(mantel_average_max, max(mantel_average), type = 'h', col = 'red')
同样使用循环，依次将UPGMA树从划分为2个簇开始，直到第n-1个簇结束（n为对象总数），获取不同数量所得簇下的二元矩阵，并计算其与原始距离矩阵（Bray-curtis距离）的Pearson相关系数，并从中获取最大的值，其所对应的划分的聚类簇数量即为评估出的最适聚类簇数量。
基于Mantel相关的评估结果显示，UPGMA树的最佳聚类簇数量均为3。
同上文，此时可继续使用cutree()，根据预期的聚类簇数量对聚类树中的对象进行合理分类后观测分类结果，以进行最终确定。
分组轮廓图
融合水平值图评估结果、轮廓宽度图评估结果、原始距离矩阵和代表分组的二元矩阵的比较评估结果均显示，建议保留的最佳簇数量（聚类分组数量）为3，这与示例数据一开始已知的分组数量一致；且裁剪聚类树观测后也得知所划分的3个分组中所包含的对象（样本）名称也与已知的一致。最终可以放心地使用UPGMA结果。
当然，若对于事先未知明确分组的样本来讲，或者在已知分组信息样本的基础上继续探索样本间可能存在的更高（或更低）层次的潜在联系时，若不同的评估方案中所给出的最佳聚类簇数量也一致，则也可据此划分聚类簇。
但是，若不同的评估方法之间存在分歧时，该如何做出判断呢？此时轮廓图可以提供一种参考，以判断不同数量的选择分组中哪些更为合理。
#轮廓图评估聚类簇数量的合理性（上述示例为 3），以 UPGMA 为例clust_average_cut <- cutree(clust_average, k = 3)
silh_final <- sortSilhouette(silhouette(clust_average_cut, dis_bray))
rownames(silh_final) <- rownames(as.matrix(dis_bray))[attr(silh_final, 'iOrd')]
plot(silh_final, main = 'UPGMA（Bray-curtis距离）-轮廓宽度图', col = clust_average_cut + 1)
若各对象的轮廓宽度值中未出现负值，则表明聚类簇数量的选择合理。
同样地，由于轮廓宽度值越大越好，此时可通过比较多种可能的结果查看，选择相对最合适的那个。
#不妨再看一下 Ward 聚类的轮廓图评估clust_ward_cut <- cutree(clust_ward, k = 3)
silh_final <- sortSilhouette(silhouette(clust_ward_cut, dis_ch))
rownames(silh_final) <- rownames(as.matrix(dis_ch))[attr(silh_final, 'iOrd')]
plot(silh_final, main = 'Ward聚类（弦距离）-轮廓宽度图', col = clust_ward_cut + 1)
此时可查看UPGMA和Ward聚类的聚类树，比较二者的区别，例如c5样本的分区。
标记聚类簇
在plot()绘制聚类树时，不妨在其中标记分组，看自动识别的簇与已知先验分组的一致性。
#UPGMA 结果为例
plot(clust_average, main = 'UPGMA\nBray-curtis distance', sub = '', xlab = 'Sample', ylab = 'Height')
rect.hclust(clust_average, k = 3, border = c('red', 'blue', 'green'))
heatmap(as.matrix(dis_bray), Rowv = TRUE , Colv = TRUE,
col = colorRampPalette(c('orange3', 'wheat'))(30),
RowSideColors = rep(c('blue', 'red', 'green3'), c(5, 5, 5)),
ColSideColors = rep(c('blue', 'red', 'green3'), c(5, 5, 5)))
参考资料