比较聚类
在上期文章层次聚类与聚类树中，不同对象之间的关系可以通过聚类树展现出来，通过聚类树我们可以观察哪些对象比较相似，哪些对象距离较远，从而对所有对象的关系有一个整体的把握。然而，这时候我们并没有获得一个明显的聚类簇划分，也即不知道对象可以划分为几类、谁和谁归为一类，以及这个聚类结果是不是合理，这可以通过比较聚类来实现。
同表型相关
同表型距离（cophenetic distance）是指聚类树上两个对象共有节点对应的距离，可以理解为聚类距离，对聚类算法依赖较大，我们可以比较同表型距离与原始距离的相关关系来评价聚类的效果，方法如下所示：
代码语言：javascript复制#读取群落数据并计算Bray-Curtis距离矩阵
data=read.table(file="sample.subsample.otu_table.txt", header=T, check.names=FALSE)
rownames(data)=data[,1]
data=data[,-1]
library(vegan)
data=decostand(data, MARGIN=2, "total")
otu=t(data)
otu_dist=vegdist(otu, method="bray", diag=TRUE, upper=TRUE, p=2)
#进行聚类分析
hclust=hclust(otu_dist, method="average")
#同表型相关
coph=cophenetic(hclust)
#计算相关系数并作图
rcor=cor(otu_dist, coph)
plot(otu_dist, coph, xlab="Bray distance", ylab="Cophenetic distance", asp=1,
xlim=c(0, 1), ylim=c(0, 1), main=c("UPGMA Clustering", paste("Cophenetic correlation", round(rcor,3))))
abline(0, 1)结果如下图所示：
可以看到聚类后的距离与原始Bray-Curtis距离相关性较好，聚类结果较为可信。
融合水平值
为了更好地比较和解读聚类结果，需要确定可解读的聚类簇数目，也即需要对聚类树层次进行修剪（聚类树最高层次聚类簇数目就是样品数）来确定有效的聚类簇数目。聚类树的融合水平值（fusion level value）是聚类树中两个分支融合处相异性的数值（该节点高度聚类簇的数目），可以绘制融合水平值变化图来确定聚类树的修剪水平，方法如下所示：
代码语言：javascript复制#总结聚类结果，查看节点数目（height）
summary(hclust)
plot(hclust$height, nrow(otu):2, type="S", xlab="h(node height)", ylab="k(number of clusters)", col="grey", lwd=2, main="Fusion Levels - UPGMA Tree")
text(hclust$height, nrow(otu):2, nrow(otu):2, col="red", cex=0.8)作图结果如下所示：
这里需要说明的是，在聚类树中节点数等于样品数减一，hclust$height里面即为节点对应的高度值（即距离）如下所示为19个节点对应的高度：
高度最大时第一个节点聚类簇数目为2，之后每增加一个节点聚类簇数目加一，高度最小（距离最小）时聚类簇数目即为样品数。因此上图就展示了不同高度水平下聚类簇数目的变化。从右往左看，随着高度的降低，聚类簇数目增加，图像呈现阶梯状。一般来说曲线越缓，“台阶”越宽，也即增加一个聚类簇间隔的距离大，其聚类约有意义。从上图来看聚类簇由2到3间隔了较宽的距离，而之后聚类簇数目快速增加，因此有效聚类簇数目应为2。
轮廓宽度
轮廓宽度（silhouette width）是描述一个对象与所属聚类簇归属程度的测度，是一个对象同同一组内其他对象的平均距离与该对象同最邻近聚类簇内对象平均距离的比较。轮廓宽度越大，对象聚类越好，轮廓宽度为负时可能存在错分情况。接下来我们计算轮廓宽度值，并确定最佳轮廓宽度值时聚类簇数目：
代码语言：javascript复制#计算轮廓宽度
library(cluster)
asw=numeric(nrow(otu))
for (i in 2:(length(asw)-1)) {
sil=silhouette(cutree(hclust, k=i), otu_dist)
asw[i]=summary(sil)$avg.width
}
#确定最佳轮廓宽度值聚类簇数目
k.best=which.max(asw)
#绘制轮廓宽度图
plot(1:length(asw), asw, type="h", main="Silhouette of UPGMA Tree", lwd=2,
xlab="k(number of clusters", ylab="average silhouette")绘图结果如下所示：
可以看出k小于5时随着聚类簇数目增加轮廓宽度下降明显，其中k=2时对应的轮廓宽度最大。
绘制聚类树
经过上面的分析，最佳聚类簇数目为2，接下来修剪聚类树，并标识不同的聚类簇，方法如下所示：
代码语言：javascript复制#根据前面分析结果确定最佳聚类簇数目，并绘制聚类树
#根据距离矩阵的样品顺序对聚类树做相应旋转，使样品排列尽可能接近原来顺序
orhclust=reorder(hclust, otu_dist)
library(dendextend)
library(RColorBrewer)
tree=as.dendrogram(orhclust)
clusMember=cutree(tree, 2)
labelColors=brewer.pal(n=2, name="Set1")
colLab=function(n) {
if (is.leaf(n)) {
a=attributes(n)
labCol=labelColors[clusMember[which(names(clusMember)==a$label)]]
attr(n, "nodePar")=c(a$nodePar, lab.col=labCol)
}
n
}
clusDendro=dendrapply(tree, colLab)
plot(clusDendro, main="UPGMA Tree", type="rectangle", horiz=TRUE)作图结果如下所示：
可以看到，这个聚类结果与样品原本的分组是一致，进一步证明了结果的合理性，因此在样品原本没有分组的情况下可以尝试通过以上步骤进行分类。
END
