A. 层次聚类分析案例(三)
之前的笔记:
聚类介绍: 点这里
层次聚类分析案例(一)
层次聚类分析案例(二)
获取全基因组表达数据的能力是一项计算复杂度非常高的任务。由于人脑的局限性,是无法解决这个问题。但是,通过将基因分类进数量较少的类别后再进行分析,就能将基因数据加工到更易理解的水平。
聚类的目标是将一组基因进行划分,使相似的基因落入同一个簇,同时不相似的基因落入不同的簇。这里需要考虑的关键问题是如何定义相似性,以及处理已分类基因。这里我们使用两种基因类型的感光性来探索基因聚类问题。
准备工作
为了进行层次聚类,我们使用从实验鼠身上采集的数据集。
第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)过程。它将一个大的问题切分成易处理的小块,在每个小块上进行操作,然后将所有小块合并起来。
载入以下包:
让我们探索并理解变量间的关系。从导入名为GSE4051_data.csv的CSV文件开始。我们将该文件数据存储到GSE4051_data数据框中:
接下来,输出GSE4051_data数据框的信息。str()函数返回GSE4051_data的结构信息。它简略显示了GSE4051_data数据框的内部结构。max.level指明了为了显示网状结构的最大等级。
结果如下:
下面,我们导入名为GSE4051_design.csv的CSV文件,将其数据保存到GSE4051_design数据框中:
输出GSE4051_design数据框的内部结构。
结果如下:
第3步:转换数据
为了便于后续的可视化阶段,需要对每一行数据进行拉伸操作。这是由于在目前的要求下,不同基因表达之间存在绝对值的差距,因此需要对每一行数据进行拉伸。
中心化变量和创建z值是两个常见的数据分析方法。scale函数中心化并拉伸数值型矩阵的列。
变换矩阵。传入GSE4051_data数据框用t()函数进行数据框变换。
接下来,我们输出GSE4051_data数据框的信息。通过设置give.attr=FALSE,次级结构的属性不会被显示。
结果如下:
round()函数用于舍入到最接近的整数。语法形式只有1种:Y = round(X),这里的X可以是数,向量,矩阵,输出对应。
head()函数返回一个向量、矩阵、表、数据框或函数的头部。GSE4051_data和trans_GSE4051_data数据框被当作对象传入。rowMeans()函数计算每列的平均值。data.frame()函数创建数据框耦合变量集合,并且共享许多指标的性质:
结果如下:
第4步:训练模型
接下来是训练模型。第一步是计算距离矩阵。dist()函数用来计算并返回距离矩阵,可以使用特定的距离度量方法来计算数据矩阵中各行间的距离。这里可使用的距离度量方法有欧式距离、最大距离、曼哈顿距离、堪培拉距离、二进制距离,或闵可夫斯基距离。这里使用欧式距离。欧式距离计算两个向量间的距离公式为sqrt(sum((x_i-y_i)^2))。转换后的trans_GSE4051_data数据框被用来计算距离。结果存储在pair_dist_GSE4051_data数据框中。
接下来,使用interaction()函数计算并返回gType、devStage变量间相互作用的无序因子。无序因子的结果连同GSE4051_design数据框一同被传入with()函数。该函数计算产生一个新的因子代表gType、devStage变量的相互作用:
summary()函数用来生成GSE4051_design$group数据框的结果总结:
结果如下:
下面,使用多种不同的联合类型计算层次聚类。
使用hclust()函数对n个不同对象进行聚类分析。第一个阶段,每个对象被指派给自己的簇。算法在每个阶段迭代聚合两个最相似的簇。持续该过程直到只剩一个单独的簇。hclust()函数要求我们以距离矩阵的形式提供数据。pair_dist_GSE4051_data数据框被传入。
在第一个例子中使用single聚类方法:
结果如下:
在第二个例子中使用complete聚集方法。
调用pr.hc.complete的结果是显示所使用的聚集方法、距离计算方法和对象数量:
结果如下:
在第三个例子中使用average聚类方法:
调用pr.hc.complete的结果是显示所使用的聚集方法、距离计算方法和对象数量:
结果如下:
在第四个例子中使用ward聚类方法:
pr.hc.ward的调用结果是显示所使用的聚集方法、距离计算方法和对象数量:
结果如下:
plot()函数是绘制R对象的通用函数。
第一次调用plot()函数,传递pr.hc.single数据框作为输入对象:
结果如下:
第二次调用plot()函数,传入pr.hc.complete数据框作为输入对象:
结果如下:
第三次调用plot()函数,传入pr.hc.average数据框作为输入对象:
结果如下:
第四次调用plot()函数,传入pr.hc.ward数据框作为输入对象:
结果如下:
第5步:绘制模型
plot()函数是绘制R对象的通用函数。这里,plot()函数用来绘制系统树图。
rect.hclust()函数强调不同的簇,并在系统树图的枝干上绘制长方形。系统树图首先在某个等级上被剪切,之后在选定的枝干上绘制长方形。
RColorBrewer使用从 http://colorbrewer2.org 获得的包来选择绘制R图像的颜色模板。
颜色分为三组:
最重要的一个RColorBrewer函数是brewer.pal()。通过向该函数传入颜色的数量和配色的名字,可以从display.brewer.all()函数中选择一个配色方案。
在第一个例子中,pr.hc.single作为一个对象传入plot()函数:
结果如下:
下面创建热度图,使用single聚集方法。heatmap()函数默认使用euclidean聚集方法:
结果如下:
在第二例子中,pr.hc.complete作为对象传入plot()函数:
结果如下:
下面使用complete聚集方法创建热度图:
结果如下:
在第三个例子中,pr.hc.average作为对象传入plot()函数:
结果如下:
下面创建average聚集方法的热度图:
结果如下:
在第四个例子中,pr.hc.ward作为对象传入plot()函数:
结果如下:
下面绘制ward聚集方法的热度图:
结果如下:
B. 在进行系统聚类分析时,不同的类间距离计算方法有何区别
聚类分析有两种主要计算方法,分别是凝聚层次聚类(Agglomerative hierarchical method)和K均值聚类(K-Means)。
一、层次聚类
层次聚类又称为系统聚类,首先要定义样本之间的距离关系,距离较近的归为一类,较远的则属于不同的类。可用于定义“距离”的统计量包括了欧氏距离 (euclidean)、马氏距离(manhattan)、 两项距离(binary)、明氏距离(minkowski)。还包括相关系数和夹角余弦。
层次聚类首先将每个样本单独作为一类,然后将不同类之间距离最近的进行合并,合并后重新计算类间距离。这个过程一直持续到将所有样本归为一类为止。在计算类间距离时则有六种不同的方法,分别是最短距离法、最长距离法、类平均法、重心法、中间距离法、离差平方和法。
下面我们用iris数据集来进行聚类分析,在R语言中所用到的函数为hclust。首先提取iris数据中的4个数值变量,然后计算其欧氏距离矩阵。然后将矩阵绘制热图,从图中可以看到颜色越深表示样本间距离越近,大致上可以区分出三到四个区块,其样本之间比较接近。
data=iris[,-5]
dist.e=dist(data,method='euclidean')
heatmap(as.matrix(dist.e),labRow = F, labCol = F)
X
然后使用hclust函数建立聚类模型,结果存在model1变量中,其中ward参数是将类间距离计算方法设置为离差平方和法。使用plot(model1)可以绘制出聚类树图。如果我们希望将类别设为3类,可以使用cutree函数提取每个样本所属的类别。
model1=hclust(dist.e,method='ward')
result=cutree(model1,k=3) 为了显示聚类的效果,我们可以结合多维标度和聚类的结果。先将数据用MDS进行降维,然后以不同的的形状表示原本的分类,用不同的颜色来表示聚类的结果。可以看到setose品种聚类很成功,但有一些virginica品种的花被错误和virginica品种聚类到一起。
C. 层次聚类分析案例(一)
关于聚类分析的介绍,可参见本人之前的笔记: 聚类分析
案例一:世界银行样本数据集
创建世界银行的一个主要目标是对抗和消除贫困。在这个不断发展的世界中,世界银行持续的发展并精细地调整它的政策,已经帮助这个机构逐渐实现了消除贫困的目标。消除贫困的成果以下指标的改进衡量,这些指标包括健康、教育、卫生、基础设施以及其他需要用于改进穷人生活的服务。与此同时,发展成果必须保证以一种环保的、全社会的、经济可持续的方式达成。
准备工作
为了进行层次聚类,我们需要使用从世界银行收集的数据集。
第1步:收集和描述数据
该任务使用名为WBClust2013的数据集。该数据以标准格式存储在名为WBClust2013.csv的CSV格式的文件中。其有80行数据和14个变量。 点我获取数据
第一列Country为非数值型变量,其他列均为数值型变量。
第2步:探索数据
让我们探索数据并理解变量间的关系。我们通过导入名为WBClust2013.csv的CSV文件开始。存储数据到wbclust数据框中:
下一步输出wbclust数据框,head()函数返回wbclust数据框。wbclust数据框作为一个输入参数传入:
结果如下:
第3步:转换数据
中心化变量和创建z值是两个常见的用于归一化数据的数据分析手段。上面提到的数值型变量需要创建z值。scale()函数是一个通用的函数,其默认方法中心化并比例缩放一个数值化矩阵的列。数据框wbclust被传给该比例函数。只有数据框中数值化的变量会被缩放。结果存储在wbnorm数据框中。
结果如下:
所有的数据框都有rownames属性。rownames()函数用来获取或设置矩阵类变量的行名或列名。数据框wbclust以及第一列被传递给rownames()函数。
调用rownames(wbnorm)方法显示第一列的数值。结果如下:
第4步:训练并评估模型效果
下一步是训练模型。首先使用dist()函数计算距离矩阵。使用特定的距离度量方法计算数据矩阵行间的距离。使用的距离度量可以是欧式距离、最大距离、曼哈顿距离、堪培拉距离、二进制距离,或闵可夫斯基距离。这里的距离度量使用欧式距离。使用欧式距离计算两个向量间的距离为sqrt(sum((x_i-y_i)^2))。结果被存储在一个新的数据框dist1中。
下一步是使用Ward方法进行聚类。hclust()函数对一组不同的n个对象进行聚类分析。第一阶段,每个对象被指派给它自己的簇。之后每个阶段,算法迭代聚合两个最相似的簇。这个过程不断持续直到只剩一个簇。hclust()函数要求我们以距离矩阵的形式提供数据。dist1数据框被作为输入传入。默认使用全链接算法。此外还可以使用不同的聚集方法,包括ward.D、ward.D2、single、complete和average。
输入clust1命令可显示所使用的聚集方法,计算距离的方法,以及数据对象的数量。结果如下:
第5步:绘制模型
plot()函数是一个通用的绘制R语言对象的函数。这里plot()函数用来绘制系统树图:
结果如下:
rect.hclust()函数强调不同的簇,并在系统树图的枝干处绘制长方形。系统树图首先在某个等级上被剪切,之后在选定的枝干上绘制长方形。
clust1对象以及需要形成的簇的数量作为输入变量传入函数。
结果如下:
cuts()函数基于期望的簇数量或者切割高度将树中的元素切割到不同的簇中。这里,clust1对象以及需要形成的簇的数量作为输入变量传入函数。
结果如下:
得到每个簇的国家列表:
结果如下: