15、R语言聚类树的绘图原理
2015-05-22 15:26
393 查看
聚类广泛用于数据分析。去年研究了一下R语言聚类树的绘图原理。以芯片分析为例，我们来给一些样品做聚类分析。聚类的方法有很多种，我们选择Pearson距离、ward方法。
选择的样品有：
R语言代码实现Pearson聚类：
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个样品放在同一个目录，然后再运行以下的代码：
解析：
一、10个样品原来的顺序：
按照hcass$order的顺序重新排列，就会得到：
这刚好是聚类图像里的样品顺序
二、再看看iia和iib：
1）第一步的iia是-8，iib是-9，如果iia或者iib的值是负数的话，说明它所代表的样品是聚类树最底层的子树的分支，我们把第8个样品和第九个样品连接起来，高度取hcl$crit的第一个值0.007874206，得到：
2）根据第2步的-6、-7和第三行的-1、-3，我们把第6个样品和第7个样品连接起来，取高度0.019655239，把第1个样品和第3个样品连接起来，取高度0.024346960：
3）第4步是-4和3，意思是把第4个样品和刚刚第三步聚类的结果（也就是第1个样品和第3个样品聚类的结果）连接起来，取高度0.025066360：
4）第五步是-2和-10，把第2个样品和第10个样品连接起来，取高度0.031698616：
5）第六步是-5和1，把样品5和第一步聚类的结果连接起来，取高度0.031710258：
6）第七步是5和6，把第五步和第一步聚类的结果连接起来，取高度0.065868858：
7）第八步连接第2步和第4步的结果，取高度0.103249166：
8）第九步连接第7步和第8步的结果，取高度0.137220473：
完成
选择的样品有：
"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 -----|
完成
内容来自用户分享和网络整理，不保证内容的准确性，如有侵权内容，可联系管理员处理 
相关文章推荐
- Android中Canvas绘图之PorterDuffXfermode使用及工作原理详解
- (15)session原理，应用（防止用户非法登录、验证码、关闭浏览器再开启浏览器还能访问之前的session）
- R语言学习第二天----R语言绘图
- Android 中 Canvas 绘图之 PorterDuffXfermode 使用及工作原理详解
- R语言聚类算法之密度聚类(Density-based Methods)
- R语言学习笔记之聚类分析
- opencv for python （15） 图像梯度（Sobel算子、scharr算子与laplacian算子原理及卷积模板）
- 慕课—R语言之数据可视化—学习笔记 3.6ggplot2绘图系统（中）
- [R语言绘图]条状图barplot
- 3D绘图过程及原理简介
- R语言绘图布局
- R语言绘图布局
- R语言低级绘图函数-rect
- Java绘图原理（一） Graphics的各种。。
- 分布式服务框架-原理与实践:15---服务降级[展示API]-学习笔记
- R语言曲线拟合函数（绘图）
- Android - View绘图原理总结
- 机器学习算法原理与实践（二）、meanshift算法图解以及在图像聚类、目标跟踪中的应用
- redis相关原理及面试官由浅到深必问的15大问题（高级）
- R语言聚类算法之系谱聚类(Hierarchical Method)