15、R语言聚类树的绘图原理

OStack程序员社区-中国程序员成长平台 › 门户 › 编程› R语言›R语言教程

原作者: [db:作者] 来自: [db:来源] 收藏邀请

　　聚类广泛用于数据分析。去年研究了一下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：

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    -----|