R語言繪製聚類樹示例

R語言繪製聚類樹示例
層次聚類(hierarchical clustering)常見兩種形式,“自底向上”的聚合策略(層次聚合)或“自頂向下”的分拆策略(層次分劃),結果一般以聚類樹表示,它表示將物件或聚類群連線在一起的層次結構。在聚類樹中,分支的高度代表了距離的遠近。
對於節點周圍分支的方向,大多數層次聚類方法中都可以任意調整順序;少數方法如TWINSPAN,物件的排列順序和其分類特徵密切相關,分支方向不可隨意調整。
在前文簡介層次聚合分類時,已經在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')

接下來,展示一些可能用到的聚類樹調整方案。
注:下文所展示的方法僅為樹狀圖本身的調整。其它組合型別的樣式,如聚類樹+柱形圖、聚類樹+熱圖、聚類樹+排序圖等,將放在後續的教程中繪製。
直接在plot()作圖時新增引數調整
基本的引數調整已在層次聚合分類時提到,以下是繼續延伸的內容。
#將樣本高度保持在同一水平,以下兩種方法都可以

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包中的系統發育樹風格
聚類樹和系統發育樹都是樹形圖,因此也可以透過系統發育樹的樣式作圖,例如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包的視覺化
例如,給分支標記顏色。
#使用 sparcl 包給分支標記顏色

library(sparcl)
ColorDendrogram(upgma, y = clusMember, labels = names(clusMember), branchlength = 0.1)

A2R包的視覺化
一個彩色聚類樹示例,指令碼也可以下載下來自定義編輯。
# 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'))

參考連結
http://rstudio-pubs-static.s3.amazonaws.com/1876_df0bf890dd54461f98719b461d987c3d.html
連結

相關文章