在前文簡介層次聚合分類時，已經在R中展示了聚類樹的一些簡單調整方法，本篇繼續作為延伸，展示一些更詳細的視覺化方案。
示例資料和R程式碼的百度盤連結：
https://pan.baidu.com/s/1ysLxEr4kOP8kEAg8HdPYEg
示例資料為來自16S測序所得的15個樣本的細菌OTU丰度表，首先執行層次聚類識別樣本歸類。
#讀取 OTU 丰度表
dat <- read.delim('otu_table.txt', row.names = 1, sep = '\t', head = TRUE, check.names = FALSE)
dat <- t(dat)
#樣本分組檔案
group <- read.delim('group.txt', row.names = 1, sep = '\t', head = TRUE, check.names = FALSE, stringsAsFactors = FALSE)
#計算樣本間距離，以群落分析中常用的 Bray-curtis 距離為例
dis_bray <- vegan::vegdist(dat, method = 'bray')
#層次聚類，以 UPGMA 為例
upgma <- hclust(dis_bray, method = 'average')
upgma
plot(upgma, main = 'UPGMA\n(Bray-curtis distance)', sub = '', xlab = 'Sample', ylab = 'Height')
接下來，展示一些可能用到的聚類樹調整方案。
注：下文所展示的方法僅為樹狀圖本身的調整。其它組合型別的樣式，如聚類樹+柱形圖、聚類樹+熱圖、聚類樹+排序圖等，將放在後續的教程中繪製。
基本的引數調整已在層次聚合分類時提到，以下是繼續延伸的內容。
#將樣本高度保持在同一水平，以下兩種方法都可以
par(mfrow = c(1, 2))
plot(upgma, hang = -1, main = 'UPGMA\n(Bray-curtis distance)', sub = '', xlab = 'Sample', ylab = 'Height')
plot(as.dendrogram(upgma), main = 'UPGMA\n(Bray-curtis distance)', sub = '', xlab = 'Sample', ylab = 'Height')
#三角形的聚類樹
plot(as.dendrogram(upgma), type = 'triangle')
#給樣本標記顏色，可以根據聚類後的簇進行標記，也可以根據先驗分組標記
#這裡按先驗分組標記
clusMember <- group$cluster
names(clusMember) <- group$samples
labelColors <- c('red', 'blue', 'green3')
#標記顏色
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(as.dendrogram(upgma), colLab)
#聚類樹
plot(clusDendro, main = 'UPGMA\n(Bray-curtis distance)', sub = '', xlab = 'Sample', ylab = 'Height')
#配合 R 的其它基礎作圖函式使用，可自定義更改聚類樹主題，例如
#繪製聚類樹主體
op <- par(bg = '#DDE3CA')
plot(upgma, col = '#487AA1', col.main = '#45ADA8', col.lab = '#7C8071',
col.axis = '#F38630', lwd = 3, lty = 3, sub = '', hang = -1, axes = FALSE)
#高度刻度軸
axis(side = 2, at = seq(0, 0.5, 0.1), col = '#F38630', labels = FALSE, lwd = 2)
mtext(seq(0, 0.5, 0.1), side = 2, at = seq(0, 0.5, 0.1), line = 1, col = '#A38630', las = 2)
#用基礎作圖函式繪製進化樹，github 上弄到的，又自己改了下，然後原出處找不到了……
treeplot <- function(tree, grp, grpcol, group_names, …) {
treeline <- function(pos1, pos2, height, col1, col2)
{
meanpos = (pos1[1] + pos2[1]) / 2
segments(y0 = pos1[1] – 0.4, x0 = -pos1[2], y1 = pos1[1] – 0.4, x1 = -height, col = col1,lwd = 2)
segments(y0 = pos1[1] – 0.4, x0 = -height, y1 = meanpos – 0.4, x1 = -height, col = col1,lwd = 2)
segments(y0 = meanpos – 0.4, x0 = -height, y1 = pos2[1] – 0.4, x1 = -height, col = col2,lwd = 2)
segments(y0 = pos2[1] – 0.4, x0 = -height, y1 = pos2[1] – 0.4, x1 = -pos2[2], col = col2,lwd = 2)
}
plot(0, type = 'n', xaxt = 'n', yaxt = 'n', frame.plot = FALSE, xlab = '', ylab = '',
ylim = c(0, length(tree$order)),
xlim = c(-max(tree$height), 0))
legend('topleft',legend = group_names,pch = 15, col = grpcol, bty = 'n', cex = 1.5)
meanpos = matrix(rep(0, 2 * length(tree$order)), ncol = 2)
meancol = rep(0, length(tree$order))
for (step in 1:nrow(tree$merge))
{
if(tree$merge[step, 1] < 0){
pos1 <- c(which(tree$order == -tree$merge[step, 1]), 0)
col1 <- grpcol[as.character(grp[tree$labels[-tree$merge[step, 1]],1])]
}else {
pos1 <- meanpos[tree$merge[step, 1], ]
col1 <- meancol[tree$merge[step, 1]]
}
if(tree$merge[step, 2] < 0){
pos2 <- c(which(tree$order == -tree$merge[step, 2]), 0)
col2 <- grpcol[as.character(grp[tree$labels[-tree$merge[step, 2]],1])]
}else {
pos2 <- meanpos[tree$merge[step, 2], ]
col2 <- meancol[tree$merge[step, 2]]
}
height <- tree$height[step]
treeline(pos1, pos2, height, col1, col2)
meanpos[step, ] <- c((pos1[1] + pos2[1]) / 2, height)
if (col1 == col2){
meancol[step] <- col1
}else {
meancol[step] <- 'grey'
}
}
tree$order
for (y in tree$order) text(x = 0, y = y, labels = rownames(grp)[y], col = grpcol[grp[y,1]])
}
grpcol <- c('red', 'blue', 'green3')
names(grpcol) <- c('1', '2', '3')
treeplot(tree = upgma, grp = group[2], grpcol = grpcol, group_names = c('A', 'B', 'C'))
聚類樹和系統發育樹都是樹形圖，因此也可以透過系統發育樹的樣式作圖，例如ape包中的方法。
library(ape)
#預設風格
plot(as.phylo(upgma))
#“無根樹”的聚類樹樣式
plot(as.phylo(upgma), type = 'unrooted')
#環形的樹
plot(as.phylo(upgma), type = 'fan')
#給樣本標記顏色，同上文，按樣本已知的先驗分組標記
plot(as.phylo(upgma), tip.color = labelColors[clusMember], label.offset = 0.01, cex = 1)
例如，給分支標記顏色。
#使用 sparcl 包給分支標記顏色
library(sparcl)
ColorDendrogram(upgma, y = clusMember, labels = names(clusMember), branchlength = 0.1)
一個彩色聚類樹示例，指令碼也可以下載下來自定義編輯。
# load code of A2R function
source('http://addictedtor.free.fr/packages/A2R/lastVersion/R/code.R')
# colored dendrogram
op <- par(bg = '#EFEFEF')
A2Rplot(upgma, k = 3, boxes = FALSE, col.up = 'gray50', col.down = c('#FF6B6B', '#4ECDC4', '#556270'))
