下文的示例資料和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')
透過繪製聚類樹，可從中直觀地辨別樣本分佈及各樣本中細菌群落組成的整體差異。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簇。
#裁剪聚類樹，詳情 ?cutree
cut_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)))
參考資料
