聚類廣泛用於數據分析。去年研究了一下R語言聚類樹的繪圖原理。以芯片分析為例，我們來給一些樣品做聚類分析。聚類的方法有很多種，我們選擇Pearson距離、ward方法。
選擇的樣品有：
"GSM658287.CEL", "GSM658288.CEL", "GSM658289.CEL", "GSM658290.CEL", "GSM658291.CEL", "GSM658292.CEL", "GSM658293.CEL", "GSM658294.CEL", "GSM658295.CEL", "GSM658296.CEL"
R語言代碼實現Pearson聚類：
library(affy) library(bioDist) rawData<-ReadAffy("GSM658287.CEL","GSM658288.CEL", "GSM658289.CEL","GSM658290.CEL", "GSM658291.CEL","GSM658292.CEL", "GSM658293.CEL","GSM658294.CEL", "GSM658295.CEL","GSM658296.CEL") correl <- cor.dist(t(exprs(rawData)),abs=FALSE) switch(tolower("pearson"), "pearson" = {correl <- cor.dist(t(exprs(rawData)),abs=FALSE)}, "spearman" = {correl <- spearman.dist(t(exprs(rawData)),abs=FALSE)}, "euclidean" = {correl <- euc(t(exprs(rawData)))}) clust <- hclust(correl, method = tolower("ward")) plot(clust)
R語言作圖結果：
根據這幾行代碼，我們只知道使用cor.dist、hclust、plot這幾個函數得到的結果，卻看不出這些函數具體做了什么，也不太有人去深究這些問題。
事實上，R語言聚類部分的計算是用Fortran實現的，源碼在https://svn.r-project.org/R/trunk/src/library/stats/src/，把hclust.f復制到本地，用一些工具生成hclust.dll，把hclust.dll和以上10個樣品放在同一個目錄，然后再運行以下的代碼：
library(affy) library(bioDist) dyn.load('hclust.dll') rawData<-ReadAffy("GSM658287.CEL","GSM658288.CEL", "GSM658289.CEL","GSM658290.CEL", "GSM658291.CEL","GSM658292.CEL", "GSM658293.CEL","GSM658294.CEL", "GSM658295.CEL","GSM658296.CEL") correl <- cor.dist(t(exprs(rawData)),abs=FALSE) ## correl是一個下三角矩陣，本文不介紹這里的重要算法 > correl GSM658287.CEL GSM658288.CEL GSM658289.CEL GSM658290.CEL GSM658288.CEL 0.037635566 GSM658289.CEL 0.024346960 0.042032944 GSM658290.CEL 0.024480935 0.040669995 0.025292084 GSM658291.CEL 0.035538210 0.039284603 0.067154783 0.048204704 GSM658292.CEL 0.072405758 0.050517381 0.059166892 0.059722043 GSM658293.CEL 0.060155354 0.060391062 0.041925320 0.043643727 GSM658294.CEL 0.036793132 0.029287344 0.069763710 0.059879668 GSM658295.CEL 0.037397535 0.030773204 0.072159149 0.060667121 GSM658296.CEL 0.068689147 0.031698616 0.068728603 0.065111592 GSM658291.CEL GSM658292.CEL GSM658293.CEL GSM658294.CEL GSM658288.CEL GSM658289.CEL GSM658290.CEL GSM658291.CEL GSM658292.CEL 0.074867692 GSM658293.CEL 0.085559588 0.019655239 GSM658294.CEL 0.023287164 0.059198270 0.071436194 GSM658295.CEL 0.028215326 0.065728329 0.075385956 0.007874206 GSM658296.CEL 0.059225037 0.046602561 0.059663628 0.044584172 GSM658295.CEL GSM658288.CEL GSM658289.CEL GSM658290.CEL GSM658291.CEL GSM658292.CEL GSM658293.CEL GSM658294.CEL GSM658295.CEL GSM658296.CEL 0.048650173 method <- 1 n <- as.integer(attr(correl, "Size")) len <- as.integer(n * (n - 1)/2) members <- rep(1, n) storage.mode(correl) <- "double" hcl <- .Fortran("hclust", n = n, len = len, method = as.integer(method), ia = integer(n), ib = integer(n), crit = double(n), members = as.double(members), nn = integer(n), disnn = double(n), flag = logical(n), diss = correl) hcass <- .Fortran("hcass2", n = n, ia = hcl$ia, ib = hcl$ib, order = integer(n), iia = integer(n), iib = integer(n)) tree <- list(merge = cbind(hcass$iia[1L:(n - 1)], hcass$iib[1L:(n - 1)]), height = hcl$crit[1L:(n - 1)], order = hcass$order, labels = attr(d, "Labels"), method = "ward", dist.method = "cor") 輸出結果： > hcl$crit[1L:(n - 1)] ## 高度 [1] 0.007874206 0.019655239 0.024346960 0.025066360 0.031698616 0.031710258 [7] 0.065868858 0.103249166 0.137220473 > hcass$iia[1L:(n - 1)] [1] -8 -6 -1 -4 -2 -5 5 2 7 > hcass$iib[1L:(n - 1)] [1] -9 -7 -3 3 -10 1 6 4 8 > hcass$order [1] 2 10 5 8 9 6 7 4 1 3
解析：
一、10個樣品原來的順序：
"GSM658287.CEL", ## 第1個 "GSM658288.CEL", ## 第2個 "GSM658289.CEL", ## 第3個 "GSM658290.CEL", ## 第4個 "GSM658291.CEL", ## 第5個 "GSM658292.CEL", ## 第6個 "GSM658293.CEL", ## 第7個 "GSM658294.CEL", ## 第8個 "GSM658295.CEL", ## 第9個 "GSM658296.CEL" ## 第10個
按照hcass$order的順序重新排列，就會得到：
"GSM658288.CEL", ## 第2個 "GSM658296.CEL", ## 第10個 "GSM658291.CEL", ## 第5個 "GSM658294.CEL", ## 第8個 "GSM658295.CEL", ## 第9個 "GSM658292.CEL", ## 第6個 "GSM658293.CEL", ## 第7個 "GSM658290.CEL", ## 第4個 "GSM658287.CEL", ## 第1個 "GSM658289.CEL", ## 第3個
這剛好是聚類圖像里的樣品順序
二、再看看iia和iib：
iia iib -8 -9 -6 -7 -1 -3 -4 3 -2 -10 -5 1 5 6 2 4 7 3
1）第一步的iia是-8，iib是-9，如果iia或者iib的值是負數的話，說明它所代表的樣品是聚類樹最底層的子樹的分支，我們把第8個樣品和第九個樣品連接起來，高度取hcl$crit的第一個值0.007874206，得到：
"GSM658288.CEL" ## 2 "GSM658296.CEL" ## 10 "GSM658291.CEL" ## 5 "GSM658294.CEL" ## 8 --| "GSM658295.CEL" ## 9 --| "GSM658292.CEL" ## 6 "GSM658293.CEL" ## 7 "GSM658290.CEL" ## 4 "GSM658287.CEL" ## 1 "GSM658289.CEL" ## 3
2）根據第2步的-6、-7和第三行的-1、-3，我們把第6個樣品和第7個樣品連接起來，取高度0.019655239，把第1個樣品和第3個樣品連接起來，取高度0.024346960：
"GSM658288.CEL" ## 2 "GSM658296.CEL" ## 10 "GSM658291.CEL" ## 5 "GSM658294.CEL" ## 8 --| "GSM658295.CEL" ## 9 --| "GSM658292.CEL" ## 6 ----| "GSM658293.CEL" ## 7 ----| "GSM658290.CEL" ## 4 "GSM658287.CEL" ## 1 -----| "GSM658289.CEL" ## 3 -----|
3）第4步是-4和3，意思是把第4個樣品和剛剛第三步聚類的結果（也就是第1個樣品和第3個樣品聚類的結果）連接起來，取高度0.025066360：
"GSM658288.CEL" ## 2 "GSM658296.CEL" ##10 "GSM658291.CEL" ## 5 "GSM658294.CEL" ## 8 --| "GSM658295.CEL" ## 9 --| "GSM658292.CEL" ## 6 ----| "GSM658293.CEL" ## 7 ----| "GSM658290.CEL" ## 4 -------| "GSM658287.CEL" ## 1 -----|_| "GSM658289.CEL" ## 3 -----|
4）第五步是-2和-10，把第2個樣品和第10個樣品連接起來，取高度0.031698616：
"GSM658288.CEL" ## 2 --------| "GSM658296.CEL" ##10 --------| "GSM658291.CEL" ## 5 "GSM658294.CEL" ## 8 --| "GSM658295.CEL" ## 9 --| "GSM658292.CEL" ## 6 ----| "GSM658293.CEL" ## 7 ----| "GSM658290.CEL" ## 4 -------| "GSM658287.CEL" ## 1 -----|_| "GSM658289.CEL" ## 3 -----|
5）第六步是-5和1，把樣品5和第一步聚類的結果連接起來，取高度0.031710258：
"GSM658288.CEL" ## 2 --------| "GSM658296.CEL" ##10 --------| "GSM658291.CEL" ## 5 ----------| "GSM658294.CEL" ## 8 --|_______| "GSM658295.CEL" ## 9 --| "GSM658292.CEL" ## 6 ----| "GSM658293.CEL" ## 7 ----| "GSM658290.CEL" ## 4 -------| "GSM658287.CEL" ## 1 -----|_| "GSM658289.CEL" ## 3 -----|
6）第七步是5和6，把第五步和第一步聚類的結果連接起來，取高度0.065868858：
"GSM658288.CEL" ## 2 --------|_____ "GSM658296.CEL" ##10 --------| | "GSM658291.CEL" ## 5 ----------|___| "GSM658294.CEL" ## 8 --|_______| "GSM658295.CEL" ## 9 --| "GSM658292.CEL" ## 6 ----| "GSM658293.CEL" ## 7 ----| "GSM658290.CEL" ## 4 -------| "GSM658287.CEL" ## 1 -----|_| "GSM658289.CEL" ## 3 -----|
7）第八步連接第2步和第4步的結果，取高度0.103249166：
"GSM658288.CEL" ## 2 --------|_____ "GSM658296.CEL" ##10 --------| | "GSM658291.CEL" ## 5 ----------|___| "GSM658294.CEL" ## 8 --|_______| "GSM658295.CEL" ## 9 --| "GSM658292.CEL" ## 6 ----|__________ "GSM658293.CEL" ## 7 ----| | "GSM658290.CEL" ## 4 -------|_______| "GSM658287.CEL" ## 1 -----|_| "GSM658289.CEL" ## 3 -----|
8）第九步連接第7步和第8步的結果，取高度0.137220473：
"GSM658288.CEL" ## 2 --------|_____ "GSM658296.CEL" ##10 --------| |____ "GSM658291.CEL" ## 5 ----------|___| | "GSM658294.CEL" ## 8 --|_______| | "GSM658295.CEL" ## 9 --| | "GSM658292.CEL" ## 6 ----|___________ | "GSM658293.CEL" ## 7 ----| |___| "GSM658290.CEL" ## 4 -------|_______| "GSM658287.CEL" ## 1 -----|_| "GSM658289.CEL" ## 3 -----|
完成