专栏首页生信宝典基因共表达聚类分析及可视化

基因共表达聚类分析及可视化

共表达基因的寻找是转录组分析的一个部分,样品多可以使用WGCNA,样品少可直接通过聚类分析如K-meansK-medoids (比K-means更稳定)或Hcluster或设定pearson correlation阈值来选择共表达基因。

下面将实战演示K-meansK-medoids聚类操作和常见问题:如何聚类分析,如何确定合适的cluster数目,如何绘制共表达密度图、线图、热图、网络图等。

获得模拟数据集

MixSim是用来评估聚类算法效率生成模拟数据集的一个R包。

library(MixSim)
# 获得5个中心点,8维属性的数据模型    
data = MixSim(MaxOmega=0,  K=5,  p=8,  ecc=0.5,  int=c(10, 100))
# 根据模型获得1000次观察的数据集
A <- simdataset(n=1000,  Pi=data$Pi, Mu=data$Mu, S=data$S, n.out=0)
data <- A$X
# 数据标准化
data <- t(apply(data, 1, scale))
# 定义行列名字
rownames(data) <- paste("Gene", 1:1000, sep="_")
colnames(data) <- letters[1:8]
head(data)
##                  a          b         c        d         e          f
## Gene_1 -0.04735251 -0.7147304 0.3836436 1.322786 0.9718988 -0.5468517
## Gene_2  0.09276733 -0.8066507 0.5476909 1.351780 0.8679073 -0.6019107
## Gene_3 -0.08751894 -0.6461075 0.4371506 1.522767 0.8031865 -0.6904081
## Gene_4  0.11065988 -0.7327674 0.4550544 1.379773 0.9304277 -0.5532253
## Gene_5 -0.02722127 -0.7833089 0.6700604 1.448916 0.7128284 -0.6266295
## Gene_6  0.15119823 -0.7468292 0.4859932 1.351159 0.9179421 -0.5625206
##                g         h
## Gene_1 0.4127370 -1.782130
## Gene_2 0.2852284 -1.736813
## Gene_3 0.3420581 -1.681128
## Gene_4 0.1808146 -1.770737
## Gene_5 0.2936467 -1.688292
## Gene_6 0.1821925 -1.779136

K-means

K-means称为K-均值聚类;k-means聚类的基本思想是根据预先设定的分类数目,在样本空间随机选择相应数目的点做为起始聚类中心点;然后将空间中到每个起始中心点距离最近的点作为一个集合,完成第一次聚类;获得第一次聚类集合所有点的平均值做为新的中心点,进行第二次聚类;直到得到的聚类中心点不再变化或达到尝试的上限,则完成了聚类过程。

聚类模拟如下图:

聚类过程需要考虑下面3点:

1.需要确定聚出的类的数目。可通过遍历多个不同的聚类数计算其类内平方和的变化,并绘制线图,一般选择类内平方和降低开始趋于平缓的聚类数作为较优聚类数, 又称elbow算法。下图中拐点很明显,5

tested_cluster <- 12
wss <- (nrow(data)-1) * sum(apply(data, 2, var))
for (i in 2:tested_cluster) {
    wss[i] <- kmeans(data, centers=i,iter.max=100,  nstart=25)$tot.withinss
}
plot(1:tested_cluster, wss, type="b", xlab="Number of Clusters", ylab="Within groups sum of squares")

2.K-means聚类起始点为随机选取,容易获得局部最优,需重复计算多次,选择最优结果。

library(cluster)
library(fpc)
# iter.max: 最大迭代次数
# nstart: 选择的随机集的数目
# centers: 上一步推测出的最优类数目
center = 5
fit <- kmeans(data, centers=center, iter.max=100, nstart=25)
withinss <- fit$tot.withinss
print(paste("Get withinss for the first run", withinss))
## [1] "Get withinss for the first run 44.381088289378"
try_count = 10
for (i in 1:try_count) {
  tmpfit <- kmeans(data, centers=center, iter.max=100, nstart=25)
  tmpwithinss <- tmpfit$tot.withinss
  print(paste(("The additional "), i, 'run, withinss', tmpwithinss))
  if (tmpwithinss < withinss){
    withins <- tmpwithinss
    fit <- tmpfit
  }
}
## [1] "The additional  1 run, withinss 44.381088289378"
## [1] "The additional  2 run, withinss 44.381088289378"
## [1] "The additional  3 run, withinss 44.381088289378"
## [1] "The additional  4 run, withinss 44.381088289378"
## [1] "The additional  5 run, withinss 44.381088289378"
## [1] "The additional  6 run, withinss 44.381088289378"
## [1] "The additional  7 run, withinss 44.381088289378"
## [1] "The additional  8 run, withinss 44.381088289378"
## [1] "The additional  9 run, withinss 44.381088289378"
## [1] "The additional  10 run, withinss 44.381088289378"
fit_cluster = fit$cluster

简单绘制下聚类效果

clusplot(data, fit_cluster, shade=T, labels=5, lines=0, color=T,
lty=4, main='K-means clusters')

3.预处理:聚类变量值有数量级上的差异时,一般通过标准化处理消除变量的数量级差异。聚类变量之间不应该有较强的线性相关关系。(最开始模拟数据集获取时已考虑)

K-medoids聚类

K-means算法执行过程,首先需要随机选择起始聚类中心点,后续则是根据聚类结点算出平均值作为下次迭代的聚类中心点,迭代过程中计算出的中心点可能在观察数据中,也可能不在。如果选择的中心点是离群点 (outlier)的话,后续的计算就都被带偏了。而K-medoids在迭代过程中选择的中心点是类内观察到的数据中到其它点的距离最小的点,一定在观察点内。两者的差别类似于平均值中值的差别,中值更为稳健。

围绕中心点划分(Partitioning Around Medoids,PAM)的方法是比较常用的(cluster::pam),与K-means相比,它有几个特征: 1. 接受差异矩阵作为参数;2. 最小化差异度而不是欧氏距离平方和,结果更稳健;3. 引入silhouette plot评估分类结果,并可据此选择最优的分类数目; 4. fpc::pamk函数则可以自动选择最优分类数目,并根据数据集大小选择使用pam还是clara (方法类似于pam,但可以处理更大的数据集) 。

fit_pam <- pamk(data, krange=2:10, critout=T)

不同的分类书计算出的silhouette值如下,越趋近于1说明分出的类越好。

## 2  clusters  0.5288058 
## 3  clusters  0.6915659 
## 4  clusters  0.8415226 
## 5  clusters  0.8661989 
## 6  clusters  0.7415207 
## 7  clusters  0.5862313 
## 8  clusters  0.4196284 
## 9  clusters  0.2518583 
## 10  clusters  0.116984

获取分类的数目

fit_pam$nc
## [1] 5
layout(matrix(c(1, 2), 1, 2)) 
plot(fit_pam$pamobject)
layout(matrix(1)) #改回每页一张图

获取分类信息

fit_cluster <- fit_pam$pamobject$clustering

数据提取和可视化

pam的输出结果为例 (上面两种方法的输出结果都已处理为了同一格式,后面的代码通用)。

1.获取每类数值的平均值,利用线图一步画图法获得各个类的趋势特征。

cluster_mean <- aggregate(data, by=list(fit_cluster), FUN=mean)
write.table(t(cluster_mean), file="ehbio.pam.cluster.mean.xls", sep='\t',col.names=F, row.names=T, quote=F)

2.获取按照分类排序的数据,绘制热图 (点击查看)

dataWithClu <- cbind(ID=rownames(data), data, fit_cluster)
dataWithClu <- as.data.frame(dataWithClu)
dataWithClu <- dataWithClu[order(dataWithClu$fit_cluster),]
write.table(dataWithClu, file="ehbio.pam.cluster.xls", 
        sep="\t", row.names=F, col.names=T, quote=F)

3.获取每类数据,绘制多线图和密度图

cluster3 <- data[fit_cluster==3,]
head(cluster3)
##                   a         b         c          d         e        f
## Gene_413 -1.2718728 0.6957162 0.9963399 -0.1895966 0.1786798 1.225407
## Gene_414 -1.1705230 0.6765085 0.8689340 -0.2155533 0.4176178 1.251818
## Gene_415 -0.9545339 0.5635188 0.8437158 -0.1360588 0.2084771 1.316728
## Gene_416 -1.0888687 0.8269888 0.7590209 -0.3090701 0.4478664 1.275057
## Gene_417 -1.1230295 0.8282559 0.9112640 -0.2524612 0.3966905 1.104951
## Gene_418 -1.1291253 0.9574904 0.8405449 -0.1200131 0.1964983 1.155913
##                    g         h
## Gene_413 -0.11021165 -1.524462
## Gene_414 -0.22266636 -1.606135
## Gene_415 -0.03562464 -1.806222
## Gene_416 -0.32047268 -1.590522
## Gene_417 -0.21171923 -1.653951
## Gene_418 -0.27792177 -1.623386
cluster3 <- t(cluster3)

多线图,绘制见线图绘制

data_rownames <- rownames(cluster3)
data_mean <- data.frame(id=data_rownames, data_mean=rowMeans(cluster3))
data_mean
##   id  data_mean
## a  a -1.1708878
## b  b  0.7811888
## c  c  0.8443212
## d  d -0.2008149
## e  e  0.2382650
## f  f  1.2601860
## g  g -0.1612029
## h  h -1.5910553
library(reshape2)
library(ggplot2)
data_melt <- melt(cluster3)
colnames(data_melt) <- c("id", "Gene", "Expr")
ggplot(data_melt, aes(id, Expr)) + stat_density_2d(aes(alpha=..level.., group=1)) + stat_smooth(data=data_mean, mapping=aes(x=id, y=data_mean, colour="red", group=1), se=F) + theme(legend.position='none')
## `geom_smooth()` using method = 'loess'

等高线的颜色越深表示对应的Y轴的点越密,对平均值的贡献越大;颜色浅的点表示分布均匀。不代表点的多少。等高线的变化趋势与平均值曲线一致。

网络图就不多说了,表达矩阵直接导入Cytoscape,跟着视频教程走就好了。

参考

  1. GIF来自 https://commons.wikimedia.org/wiki/File:Kmeans_animation.gif
  2. PAM更多描述见http://shiyanjun.cn/archives/1398.html

本文分享自微信公众号 - 生信宝典(Bio_data),作者:陈同

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2017-12-01

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 数据可视化完美指南-R-python

    Yan Holtz 和Conor Healys两个人关系很好,一起在业余时间开发了这个网站。基于R和Python做的源代码,这里我们不仅可以得到大量优秀的源代码...

    生信宝典
  • R语言学习 - 热图美化

    热图美化 上一期的绘图命令中,最后一行的操作抹去了之前设定的横轴标记的旋转,最后出来的图比较难看。 上次我们是这么写的 p <- p + xlab("samp...

    生信宝典
  • 你知道R中的赋值符号箭头(<-)和等号(=)的区别吗?

    作为一门高级语言,R语言拥有独特的语法,比如今天说道的赋值符号。在其他语言里,赋值符合通常用一个等号(=)表示,而在R语言里,承担这个任务的可以是箭头(<-)符...

    生信宝典
  • R语言构建层次分析模型不看一下吗~

    AHP (Analytic Hierarchy Process)层次分析法是美国运筹学家Saaty教授于二十世纪80年代提出的一种实用的多方案或多目标的决策方法...

    数据小磨坊
  • 利用BootStrap Table插件实现自己的弹出框分页。

    参考链接1:     官网:http://bootstrap-table.wenzhixin.net.cn/zh-cn/home/         开始使用:h...

    别先生
  • [机智的机器在学习] 利用TensorFlow实现多元线性回归分类器

    从今天的推文开始,我打算把经典的机器学习算法,都用tf实现一遍。这样一来可以熟悉一下机器学习算法,二来可以对tf有比较好的掌握,如果你是新手,那就跟着我的节奏,...

    用户1622570
  • 微信小程序 计算器

    text-shadow:0 1px 1px rgba(255,255,255,.3);

    用户5760343
  • R-forestplot包| HR结果绘制森林图

    上一篇简单的介绍了COX生存分析结果绘制森林图Forest plot(森林图) | Cox生存分析可视化,本文将介绍根据数据集合的基本信息以及点估计值(置信区间...

    西游东行
  • ActiveReports 报表应用教程 (11)---交互式报表之文档目录

    通过文档目录,用户可以非常清晰的查看报表数据结构,并能方便地跳转到指定的章节,最终还可以将报表导出为PDF等格式的文件。本文以2012年各月产品销售分类汇总报表...

    葡萄城控件
  • hand first python 选读(3)

    分析数据,有的是以 .连接,有些以 -连接,还有的是 :。需作分类讨论转化为小数。

    一粒小麦

扫码关注云+社区

领取腾讯云代金券