基因芯片（Affymetrix）分析5：聚类分析

聚类又称非监督分类，是一种探索性数据分析（Exploratory Data Analysis, EDA）方法，其目的是把有限数量的事物划分到若干“自然”的类别中。芯片结果的聚类分析很复杂，因为聚类本身就是个复杂的问题，聚类后的结果解析则更复杂。

在芯片数据分析中为什么要进行基因的聚类分析？主要是为了了解未知基因的功能。通过聚类分析我们把基因按某些特性（表达规律）分成若干类别，聚合成同一类的基因在表达上有相似规律，它们很可能具有某些类似的功能。基于这一假设，我们可以通过了解某一类别中已知功能的基因进一步分析其他未知基因的功能。

聚类分析方法（算法）大体上可分为层次聚类（hierarchical clustering）和分配聚类（partitional clustering）两种，其中层次聚类在芯片分析中用得比较多。经典的层次聚类也有两种：凝聚式层次聚类法（agglomerative hierarchical clustering）和分裂式层次聚类（divisive hierarchical clustering）。凝聚式层次聚类法先把N个基因当成N个独立的簇（或称群、类），通过逐个合并最相似的项最终归为一个簇；而分裂式层次聚类正好相反。

R语言发展到今天，研究者为聚类分析软件包的开发做出了大量的贡献（衷心地感谢他们），使用者有很多选择。读者可以从CRAN的Task view页面做进一步的了解。

1 计算距离：

先读入我们前面保存的分析数据：

data.fc <- read.csv("results.lim.all.csv", header=T, row.names=1)
head(data.fc, 2)

##               T1h  T24h   T7d
## 254818_at  0.3409 6.024 6.215
## 245998_at -0.1367 3.676 2.778

R中计算距离的函数为dist( )。但是使用什么指标来计算距离呢？一般通过计算相关系数作为基因间距离，R函数是cor( )。使用fold change计算距离也可以做出图来，其结果显然和用相关系数做出的图不一样：

par(mfrow=c(2,1))
par(mar=c(0,4,2,0))
dist.fc <- dist(data.fc, method="euclidean")
plot(hclust(dist.fc), hang=-1, labels=FALSE, main="使用Fold-change计算距离", sub="")
dist.fc <- as.dist(1 - cor(t(data.fc)))
plot(hclust(dist.fc), hang=-1, labels=FALSE, main="使用相关系数作为距离", sub="")

cor函数计算相关系数是计算矩阵不同列的相关系数，上面代码中函数 t() 对矩阵行列进行了转置，这样计算出的是基因间的相关系数矩阵。设置cor函数的method参数可以计算出pearson, kendall 或 spearman相关系数矩阵，默认计算pearson相关系数。as.dist函数把该矩阵转成了距离（向量）。如果直接使用dist函数计算距离，method可以设置距离计算的类型。

2 聚类

上面使用的hclust函数是R基本包中的凝聚式层次聚类函数。分裂式层次聚类运算耗时很长，做高通量数据分析一般不采用。

凝聚式层次聚类方法首先把N个基因当成N个独立的簇，每合并出一个新簇就为它计算相似系数。那么这个新簇的系数怎么计算？比如Ci和Cj合并后产生Cij，Cij在相似矩阵里面取什么值？这不是生物学研究者该考虑的问题，我们只用现成的。R的hclust函数提供了7中可用方法，"ward", "single", "complete", "average", "mcquitty", "median" 和 "centroid"。

我们看看其中的6种结果：

par(mfrow=c(3,2))
par(mar=c(0,4,2,0))
plot(hclust(dist.fc, method<