Tag Archives: 聚类

八 27

根据基因表达量对样品进行分类ConsensusClusterPlus

Posted on by 
bioconductor系列的包都是一样的安装方式：
source("http://bioconductor.org/biocLite.R")
biocLite("ConsensusClusterPlus")

这个包是我见过最简单的包， 加载只有做好输入数据，只需要一句话即可运行，然后默认输出所有结果

读这个包的readme，很容易学会
就是做好一个需要来进行分类的样品的表达量矩阵。或者选择上一篇日志用GEOquery这个包下载的表达量矩阵也可以进行分析
因为这个包是用ALL数据来做测试的，所以可以直接加载这个数据结果，这样就能得到表达矩阵啦
library(ALL)
data(ALL)
d=exprs(ALL)
d[1:5,1:5]
可以看到数据集如下

> d[1:5,1:5]

             01005    01010    03002    04006    04007
1000_at   7.597323 7.479445 7.567593 7.384684 7.905312
1001_at   5.046194 4.932537 4.799294 4.922627 4.844565
1002_f_at 3.900466 4.208155 3.886169 4.206798 3.416923
1003_s_at 5.903856 6.169024 5.860459 6.116890 5.687997
1004_at   5.925260 5.912780 5.893209 6.170245 5.615210
> dim(d)
[1] 12625   128
共128个样品，12625个探针数据
也有文献用RNAs-seq的RPKM值矩阵来做
对上面这个芯片表达数据我们一般会简单的进行normalization ，然后取在各个样品差异很大的那些gene或者探针的数据来进行聚类分析
mads=apply(d,1,mad)
d=d[rev(order(mads))[1:5000],]

d = sweep(d,1, apply(d,1,median,na.rm=T))

#也可以对这个d矩阵用DESeq的normalization 进行归一化，取决于具体情况
library(ConsensusClusterPlus)
#title=tempdir() #这里一般改为自己的目录
title="./" #所有的图片以及数据都会输出到这里的
results = ConsensusClusterPlus(d,maxK=6,reps=50,pItem=0.8,pFeature=1,
 title=title,clusterAlg="hc",distance="pearson",seed=1262118388.71279,plot="png")
这样就OK了，你指定的目录下面会输出大于9个图片

clipboard

大家看看说明书就知道这个包的输出文件是什么了。
很多参数都是需要调整的，一般我们的maxK=6是根据实验原理来调整，如果你的样品应该是要分成6类以上，那么你就要把maxK=6调到一点。
查看结果results[[2]][["consensusClass"] 可以看到各个样品被分到了哪个类别里面去
results[[3]][["consensusClass"]
results[[4]][["consensusClass"] 等等
七 21

R语言用hclust进行聚类分析

Posted on by

聚类的基础就是算出所有元素两两间的距离，我们首先做一些示例数据，如下：

x=runif(10)

y=runif(10)

S=cbind(x,y)                                 #得到2维的数组

rownames(S)=paste("Name",1:10,"")             #赋予名称，便于识别分类

out.dist=dist(S,method="euclidean")           #数值变距离

这个代码运行得到的S是一个矩阵，如下

> S

x         y

Name 1   0.41517985 0.4697017

Name 2   0.35653781 0.1132367

Name 3   0.52253349 0.3680286

Name 4   0.80558684 0.9834687

Name 5   0.04564145 0.8560690

Name 6   0.11044397 0.2988598

Name 7   0.34984447 0.8515141

Name 8   0.28097709 0.1260050

Name 9   0.81771888 0.5976135

Name 10 0.40700158 0.5236567

可以看出里面共有10个点，它们的X,Y坐标均已知，我们有6总方法可以求矩阵

image001

注释：在聚类中求两点的距离有：

1，绝对距离：manhattan

2，欧氏距离：euclidean 默认

3，闵科夫斯基距离：minkowski

4，切比雪夫距离：chebyshev

5，马氏距离：mahalanobis

6，蓝氏距离：canberra

用默认的算法求出距离如下

算出距离后就可以进行聚类啦！

out.hclust=hclust(out.dist,method="complete") #根据距离聚类

注释：聚类也有多种方法：

1，类平均法：average

2，重心法：centroid

3，中间距离法:median

4，最长距离法：complete 默认

5，最短距离法：single

6，离差平方和法：ward

7，密度估计法：density

接下来把聚类的结果图画出来

plclust(out.hclust)                           #对结果画图

image003

rect.hclust(out.hclust,k=3)                   #用矩形画出分为3类的区域

image005

out.id=cutree(out.hclust,k=3)                 #得到分为3类的数值

这里的out.id就是把每个点都分类了的分类数组，1,2,3.