之前的笔记：
聚类介绍：点这里
层次聚类分析案例（一）
层次聚类分析案例（二）

案例三：基因聚类

获取全基因组表达数据的能力是一项计算复杂度非常高的任务。由于人脑的局限性，是无法解决这个问题。但是，通过将基因分类进数量较少的类别后再进行分析，就能将基因数据加工到更易理解的水平。

聚类的目标是将一组基因进行划分，使相似的基因落入同一个簇，同时不相似的基因落入不同的簇。这里需要考虑的关键问题是如何定义相似性，以及处理已分类基因。这里我们使用两种基因类型的感光性来探索基因聚类问题。

准备工作

为了进行层次聚类，我们使用从实验鼠身上采集的数据集。

第1步：收集和描述数据

该任务使用名为GSE4051_data和GSE4051_design的数据集。该数据集以标准格式存储在名为GSE4051_data.csv和GSE4051_design.csv的CSV格式的文件中。数据获取路径： 在这里

GSE4051_data数据集包含29949行数据和39个变量。数值型变量如下：

GSE4051_design数据集包含39行数据和4个变量。数值型变量是：sidNum
非数值型变量是：sidChar；devStage；gType；

具体实施步骤以下为实现细节。

第2步：探索数据

RColorBrewer包是一个R包，可从http://colorbrewer2.org获取，它提供地图和其他图形的彩色模板。

pvclust包用来实现非确定性的层次聚类分析。在层次聚类中，每个簇通过多尺度有放回抽样计算p值。一个簇的p值在0～1之间。p值有两种类型：近似无偏（approximately unbiased，AU）和有放回概率（bootstrap probability，BP）值。AU p值通过多尺度有放回采样方法计算，经典的有放回采样方法用来计算BP p值。AU p值相比BP p值存在优效性偏见。

xtable包可以生成LaTeX格式的表格。使用xtable可以将特定的R对象转换成xtables。这些xtables能够以LaTeX或HTML的格式输出。

plyr包被用来进行分置合并（split-apply-combine，SAC）过程。它将一个大的问题切分成易处理的小块，在每个小块上进行操作，然后将所有小块合并起来。

载入以下包：

library(RColorBrewer)
library(cluster)
library(pvclust)
library(xtable)
library(plyr)

让我们探索并理解变量间的关系。从导入名为GSE4051_data.csv的CSV文件开始。我们将该文件数据存储到GSE4051_data数据框中：

GSE4051_data = read.csv("ClusteringAnalysis/Practical-Machine-Learning-Cookbook/Chapter03/Data/GSE4051_data.csv",header = T)

接下来，输出GSE4051_data数据框的信息。str（）函数返回GSE4051_data的结构信息。它简略显示了GSE4051_data数据框的内部结构。max.level指明了为了显示网状结构的最大等级。

str(GSE4051_data, max.level = 0)

结果如下：

下面，我们导入名为GSE4051_design.csv的CSV文件，将其数据保存到GSE4051_design数据框中：

GSE4051_design = read.csv("ClusteringAnalysis/Practical-Machine-Learning-Cookbook/Chapter03/Data/GSE4051_design.csv",header = T)

输出GSE4051_design数据框的内部结构。

str(GSE4051_design)

结果如下：

第3步：转换数据

为了便于后续的可视化阶段，需要对每一行数据进行拉伸操作。这是由于在目前的要求下，不同基因表达之间存在绝对值的差距，因此需要对每一行数据进行拉伸。

中心化变量和创建z值是两个常见的数据分析方法。scale函数中心化并拉伸数值型矩阵的列。

变换矩阵。传入GSE4051_data数据框用t（）函数进行数据框变换。

trans_GSE4051_data <- t(scale(t(GSE4051_data)))

接下来，我们输出GSE4051_data数据框的信息。通过设置give.attr=FALSE，次级结构的属性不会被显示。

str(trans_GSE4051_data,max.level=0, give.attr = FALSE)

结果如下：

num [1:29949, 1:39] 0.0838 0.1758 0.7797 -0.3196 0.8358 ...

round（）函数用于舍入到最接近的整数。语法形式只有1种：Y = round(X)，这里的X可以是数，向量，矩阵，输出对应。

head（）函数返回一个向量、矩阵、表、数据框或函数的头部。GSE4051_data和trans_GSE4051_data数据框被当作对象传入。rowMeans（）函数计算每列的平均值。data.frame（）函数创建数据框耦合变量集合，并且共享许多指标的性质：

round(data.frame(avgBefore = rowMeans(head(GSE4051_data)),
avgAfter = rowMeans(head(trans_GSE4051_data)),
varBefore = apply(head(GSE4051_data),1,var),
varAfter = apply(head(trans_GSE4051_data),1,var)),2)

结果如下：

第4步：训练模型

接下来是训练模型。第一步是计算距离矩阵。dist（）函数用来计算并返回距离矩阵，可以使用特定的距离度量方法来计算数据矩阵中各行间的距离。这里可使用的距离度量方法有欧式距离、最大距离、曼哈顿距离、堪培拉距离、二进制距离，或闵可夫斯基距离。这里使用欧式距离。欧式距离计算两个向量间的距离公式为sqrt（sum（（x_i-y_i）^2））。转换后的trans_GSE4051_data数据框被用来计算距离。结果存储在pair_dist_GSE4051_data数据框中。

pair_dist_GSE4051_data <- dist(t(trans_GSE4051_data),method = "euclidean")

接下来，使用interaction（）函数计算并返回gType、devStage变量间相互作用的无序因子。无序因子的结果连同GSE4051_design数据框一同被传入with（）函数。该函数计算产生一个新的因子代表gType、devStage变量的相互作用：

GSE4051_design$group <- with(GSE4051_design,interaction(gType,devStage))

summary（）函数用来生成GSE4051_design$group数据框的结果总结：

summary(GSE4051_design$group)

结果如下：

下面，使用多种不同的联合类型计算层次聚类。

使用hclust（）函数对n个不同对象进行聚类分析。第一个阶段，每个对象被指派给自己的簇。算法在每个阶段迭代聚合两个最相似的簇。持续该过程直到只剩一个单独的簇。hclust（）函数要求我们以距离矩阵的形式提供数据。pair_dist_GSE4051_data数据框被传入。

在第一个例子中使用single聚类方法：

pr.hc.single <- hclust(pair_dist_GSE4051_data,method = "single")
pr.hc.single的调用结果是现实使用的聚集方法、距离计算方法和对象数量：
pr.hc.single

结果如下：

在第二个例子中使用complete聚集方法。

pr.hc.complete <- hclust(pair_dist_GSE4051_data,method = "complete")

调用pr.hc.complete的结果是显示所使用的聚集方法、距离计算方法和对象数量：

pr.hc.complete

结果如下：

在第三个例子中使用average聚类方法：

pr.hc.average <- hclust(pair_dist_GSE4051_data,method = "average")
pr.hc.average

调用pr.hc.complete的结果是显示所使用的聚集方法、距离计算方法和对象数量：
结果如下：

在第四个例子中使用ward聚类方法：

pr.hc.ward <- hclust(pair_dist_GSE4051_data,method = "ward.D2")
pr.hc.ward

pr.hc.ward的调用结果是显示所使用的聚集方法、距离计算方法和对象数量：
结果如下：

plot（）函数是绘制R对象的通用函数。

第一次调用plot（）函数，传递pr.hc.single数据框作为输入对象：

plot(pr.hc.single,labels=FALSE, main="Single Linkage Representation", xlab="")

结果如下：

第二次调用plot（）函数，传入pr.hc.complete数据框作为输入对象：

plot(pr.hc.complete,labels=FALSE, main="Complete Linkage Representation", xlab="")

结果如下：

第三次调用plot（）函数，传入pr.hc.average数据框作为输入对象：

plot(pr.hc.average, labels=FALSE, main="Arverage Linkage Representation", xlab="")

结果如下：

第四次调用plot（）函数，传入pr.hc.ward数据框作为输入对象：

plot(pr.hc.ward, labels=FALSE, main="Ward Linkage Representation", xlab="")

结果如下：

第5步：绘制模型

plot（）函数是绘制R对象的通用函数。这里，plot（）函数用来绘制系统树图。
rect.hclust（）函数强调不同的簇，并在系统树图的枝干上绘制长方形。系统树图首先在某个等级上被剪切，之后在选定的枝干上绘制长方形。
RColorBrewer使用从http://colorbrewer2.org获得的包来选择绘制R图像的颜色模板。
颜色分为三组：

• 时序：低数据——浅色；高数据——深色。
• 分歧：中间数据——浅色；低和高范围数据——相反的深色。
• 定性的：设计颜色以强调不同簇之间的最大视觉差。

最重要的一个RColorBrewer函数是brewer.pal（）。通过向该函数传入颜色的数量和配色的名字，可以从display.brewer.all（）函数中选择一个配色方案。
在第一个例子中，pr.hc.single作为一个对象传入plot（）函数：

plot(pr.hc.single, labels = GSE4051_design$group, cex = 0.6, main = "Single Hierarchical Cluster - 10 clusters")
rect.hclust(pr.hc.single,k=10)

结果如下：

下面创建热度图，使用single聚集方法。heatmap（）函数默认使用euclidean聚集方法：

op <- par(mar=c(1,4,4,1))
par(op)
jGraysFun <- colorRampPalette(brewer.pal(n=9, "Blues"))
gTypeCols <- brewer.pal(9,"Spectral")[c(4,7)]

heatmap(as.matrix(trans_GSE4051_data),Rowv = NA, col = jGraysFun(256),
    hclustfun = function(x) hclust(x, method = 'single'),
    scale = "none", labCol = GSE4051_design$group, labRow = NA, margins = c(8,1),
    ColSideColors = gTypeCols[unclass(GSE4051_design$gType)])
legend("topright",legend=levels(GSE4051_design$gType),col = gTypeCols, lty = 1, lwd = 5, cex = 0.5)

结果如下：

image.png

在第二例子中，pr.hc.complete作为对象传入plot（）函数：

plot(pr.hc.complete, labels = GSE4051_design$group, cex = 0.6, main = "Complete Hierarchical Cluster - 10 clusters")
rect.hclust(pr.hc.complete,k=10)

结果如下：

下面使用complete聚集方法创建热度图：

par(op)
jGraysFun <- colorRampPalette(brewer.pal(n=9, "Greens"))
gTypeCols <- brewer.pal(11,"PRGn")[c(4,7)]

heatmap(as.matrix(trans_GSE4051_data),Rowv = NA, col = jGraysFun(256),
        hclustfun = function(x) hclust(x, method = 'complete'),
        scale = "none", labCol = GSE4051_design$group, labRow = NA, margins = c(8,1),
        ColSideColors = gTypeCols[unclass(GSE4051_design$gType)])

legend("topright",legend=levels(GSE4051_design$gType),col = gTypeCols,
       lty = 1, lwd = 5, cex = 0.5)

结果如下：

在第三个例子中，pr.hc.average作为对象传入plot（）函数：

plot(pr.hc.average, labels = GSE4051_design$group, cex = 0.6, main = "Average Hierarchical Cluster - 10 clusters")
rect.hclust(pr.hc.average, k=10)

结果如下：

下面创建average聚集方法的热度图：

jGraysFun <- colorRampPalette(brewer.pal(n=9, "Oranges"))
gTypeCols <- brewer.pal(9,"Oranges")[c(4,7)]

heatmap(as.matrix(trans_GSE4051_data),Rowv = NA, col = jGraysFun(256),
        hclustfun = function(x) hclust(x, method = 'average'),
        scale = "none", labCol = GSE4051_design$group, labRow = NA, margins = c(8,1),
        ColSideColors = gTypeCols[unclass(GSE4051_design$gType)])

legend("topright",legend=levels(GSE4051_design$gType),col = gTypeCols,
       lty = 1, lwd = 5, cex = 0.5)

结果如下：

在第四个例子中，pr.hc.ward作为对象传入plot（）函数：

plot(pr.hc.ward, labels = GSE4051_design$group,
      cex = 0.6, main = "Ward Hierarchical Cluster - 10 clusters")
rect.hclust(pr.hc.ward,k=10)

结果如下：

下面绘制ward聚集方法的热度图：

jGraysFun <- colorRampPalette(brewer.pal(n=9, "Reds"))
gTypeCols <- brewer.pal(9,"Reds")[c(4,7)]

heatmap(as.matrix(trans_GSE4051_data),Rowv = NA, col = jGraysFun(256),
        hclustfun = function(x) hclust(x, method = 'ward.D2'),
        scale = "none", labCol = GSE4051_design$group, labRow = NA, margins = c(8,1),
        ColSideColors = gTypeCols[unclass(GSE4051_design$gType)])

legend("topright",legend=levels(GSE4051_design$gType),col = gTypeCols,
       lty = 1, lwd = 5, cex = 0.5)

结果如下：