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

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

完成

时间： 2024-11-05 12:32:37

15、R语言聚类树的绘图原理的相关文章

R语言学习笔记2——绘图

R语言提供了非常强大的图形绘制功能.下面来看一个例子: > dose <- c(20, 30, 40, 45, 60)> drugA <- c(16, 20, 27, 40, 60)> drugB <- c(15, 18, 25, 31, 40) > plot(dose, drugA, type="b") > plot(dose, drugB, type="b") 该例中,我们引入了R语言中第一个绘图函数plot.pl

R语言聚类（K-Means、层次）

R语言聚类 K-Means 1. 随机生成3个簇点 > c1=cbind(rnorm(20,2,1),rnorm(20,2,1)) > c2=cbind(rnorm(20,3,2),rnorm(20,15,3)) > c3=cbind(rnorm(20,20,2),rnorm(20,20,3)) > v=rbind(c1,c2,c3) 在图中看看这三个簇的分布 > plot(v) 如图, 2. K聚类 clara(x, k),K聚类函数 x是数据集,可以是矩阵或者数据框 k是

R语言数据挖掘 — 决策树直观绘图

R语言数据挖掘 - 决策树直观绘图 1 前言 今天发现一个特别漂亮的决策树绘图方法,特此记录下来,作图工具是R语言,方法特别简单,图形直观美丽大方让我眼界大开. 2 安装包准备 绘制这些漂亮的图需要安装下面的包: library(rpart) library(rattle) library(rpart.plot) library(RColorBrewer) 上面是加载语言,这些包都要 install.packages 安装 3 测试代码 model <- rpart(Species ~ Sepa

R语言基于树的方法：决策树，随机森林，套袋Bagging，增强树

原文链接:http://tecdat.cn/?p=9859 概观 本文是有关  基于树的  回归和分类方法的.用于分割预测变量空间的分割规则可以汇总在树中,因此通常称为  决策树  方法. 树方法简单易懂,但对于解释却非常有用,但就预测准确性而言,它们通常无法与最佳监督学习方法竞争.因此,我们还介绍了装袋,随机森林和增强.这些示例中的每一个都涉及产生多个树,然后将其合并以产生单个共识预测.我们看到,合并大量的树可以大大提高预测准确性,但代价是损失解释能力. 决策树可以应用于回归和分类问题.我们将

R语言：ggplot2精细化绘图——以实用商业化图表绘图为例

本文旨在介绍R语言中ggplot2包的一些精细化操作,主要适用于对R画图有一定了解,需要更精细化作图的人,尤其是那些刚从excel转ggplot2的各位,有比较频繁的作图需求的人.不讨论那些样式非常酷炫的图表,以实用的商业化图表为主.包括以下结构: 1.画图前的准备:自定义ggplot2格式刷 2.画图前的准备:数据塑形利器dplyr / tidyr介绍 3.常用的商业用图: 1)简单柱形图+文本(单一变量) 2)分面柱形图(facet_wrap/facet_grid) 3)簇型柱形图(posi

R语言中文社区历史文章整理（类型篇）

R语言中文社区历史文章整理(类型篇) R包: R语言交互式绘制杭州市地图:leafletCN包简介 clickpaste包介绍 igraph包快速上手 jiebaR,从入门到喜欢 Catterplots包,让你绘制不一样的图 今天再来谈谈REmap包 ggplot2你需要知道的都在这... R访问数据库管理系统(通过RODBC包和RMySQL包两种方式) NLP--自然语言处理(三)text2vec包 Rattle:数据挖掘的界面化操作 借助caret包实现特征选择的工作 R语言的高质量图形渲染

R语言基础编程技巧汇编 - 17

1.       timestamp函数输出当前时间 timestamp() ##------ Sun Apr 05 20:54:06 2015 ------## 该函数可以输入当前的系统时间,可用于耗时很长的程序定时输出当前时间,用于判断程序是否正常运行:也可用于调试,判断哪一段代码效率较低. 2.       多个比较的boxplot图 a=c(1,2,3,4,5,2,1,2,4,2,5,6) b=c("A","A","B","B&

主成分分析（PCA）原理及R语言实现

原理: 主成分分析 - stanford 主成分分析法 - 智库 主成分分析(Principal Component Analysis)原理 主成分分析及R语言案例 - 文库 主成分分析法的原理应用及计算步骤 - 文库 主成分分析之R篇 [机器学习算法实现]主成分分析(PCA)——基于python+numpy scikit-learn中PCA的使用方法 Python 主成分分析PCA 机器学习实战-PCA主成分分析.降维(好) 关于主成分分析的五个问题 主成分分析(PCA)原理详解(推荐) 多变

R语言学习笔记-绘图相关

绘图函数变量解释: 1.绘图函数 plot  高级绘图函数,能自动创建新的绘图窗口 lines,points 低级绘图函数,用于在已有图形上叠加新的图形. lengend 图例,低级绘图函数 下面的例子画cpu1,c2两条曲线: plot(cpu1,type="o",pch=15,lty=1,col="blue") lines(c2,type="o",pch=1,lty=1,col="red") legend("to