2024-06-17

跟着小花Get无监督聚类！-ConsensusClusterPlus

小花生信果 2023-08-09 19:00:35

转自公众号：生信花
http://mp.weixin.qq.com/s?__biz=MzkzODUxNjcwMA==&mid=2247483989&idx=1&sn=5f55d4242616466918c5c775786d6eb0

点击蓝字关注我们

小伙伴在做生信分析的时候，会经常见到聚类这个名词，却不知道是什么意义。其实平时我们做基因热图时候就用到了聚类的思想。有的小伙伴在做癌症分析时候，也见过给疾病做分子亚型，这个时候无监督聚类就派上用场了。

所以今天小花带大家了解无监督聚类，并学习如何使用它。

首先我们先大致了解一下什么是无监督聚类。就好比一致性聚类方法包括从一组项目中进行次抽样，例如微阵列，并确定特定簇计数（k）的簇。然后，对共识值，两个项目占在同一子样本中发生的次数中有相同的聚类，计算并存储在每个k的对称一致矩阵中。

具体的含义我们不在此一一赘述。

小花在这里用到的是ConsensusClusterPlus这个包对基因进行聚类。

所需要的是下面这种格式的基因表达谱：

#install.packages("BiocManager")#BiocManager::install('ConsensusClusterPlus')#小伙伴没有的自己要下载一下library(tidyverse)library(ConsensusClusterPlus)#导入数据exp <- read.table("text1.txt",sep = "t",row.names = 1,check.names = F,stringsAsFactors = F,header = T)d=as.matrix(exp)gene <- c("MMP9","CCND1","STAT1","AR","AKR1C3","GSTP1","PGR")d <- d[gene,]mads=apply(d,1,mad)d=d[rev(order(mads))[1:7],] #7是基因个数d = sweep(d,1, apply(d,1,median,na.rm=T))title=("JULEI") ##文件夹输出图片的位置set.seed(1) #我发现设不设置种子都一样#这里小花选取自己的示例数据，小伙伴们整理自己想去聚类的基因。

然后我们打开JULEI这个文件夹会发现得到了聚类的结果：

一致性矩阵

一致性累积分布函数图

碎石图

在这里小花就不全部展示了，小伙伴们自己做的话可以慢慢看。

从上述碎石图看出呢，在k=7,k=8左右时候呢，比较平稳，前面K=3，和k=4这两个分组数间相对差值变化越大，说明该分组越有价值。小伙伴需要对症下药，查看自己聚类的结果进行解释和后续的分析。

图中K值呢就代表簇，

接下来我们绘制另一组图：

聚类一致性(cluster-consensus)以及样品一致性(item-consensus)

title=title,clusterAlg=”hc”,distance=”pearson”,seed=1,plot=”pdf”)

results[[2]][[“consensusMatrix”]][1:5,1:5]

results[[2]][[“consensusTree”]]

results[[2]][[“consensusClass”]][1:5]

icl = calcICL(results,title=title,plot=”pdf”) ##画另一组图片

结果如下：

聚类一致性(cluster-consensus)

样品一致性(item-consensus)

小伙伴们把需要那些图片进行展示就可以了。

#上述就是绘制另一组图，然而我们还需要做的是把结果保存下来。

group<-results[[2]][[“consensusClass”]]

group<-as.data.frame(group)

group$group <- factor(group$group,levels=c(1,2))

save(group,file = “group_AY.Rda”)#这里把分组结果保存

然而，又会有小伙伴会问如何展示聚类后的效果呢。

小花在这里教大家一种方式用热图的方式进展展示，

绘制ConsensusClusterPlus后的热图：

这个时候就用到上述咱们保存的分组信息了：

跟着小花的代码往下跑就可以啦！

library(pheatmap)#调用一下热图的包group <- group %>%rownames_to_column("sample")annotation <- group %>% arrange(group) %>% column_to_rownames("sample")a <- group %>% arrange(group) %>% mutate(sample=substring(.$sample,1,12))b <- t(exp_gene) %>%as.data.frame() %>%rownames_to_column("sample") %>%mutate(sample=substring(.$sample,1,12))c <- inner_join(a,b,"sample") %>% .[,-2] %>% column_to_rownames("sample") %>% t(.)pheatmap(c,annotation = annotation,cluster_cols = F,fontsize=5,fontsize_row=5,scale="row",show_colnames=F,fontsize_col=3)bk <- c(seq(-5,-0.1,by=0.01),seq(0,5,by=0.01))pheatmap(c,annotation = annotation,annotation_colors = list(group = c("1" ="#8B1C62","2"= "#EE9572")),cluster_cols = F,fontsize=10,fontsize_row=10,scale="row",show_colnames=F,cluster_row = F,color = c(colorRampPalette(colors = c("#551A8B","white"))(length(bk)/2),colorRampPalette(colors = c("white","#CD3700"))(length(bk)/2)),#legend_breaks=seq(-3,3,1),fontsize_col=3)dev.off()