前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >使用NMF代替层次聚类

使用NMF代替层次聚类

作者头像
生信技能树
发布2020-06-09 14:54:09
3K0
发布2020-06-09 14:54:09
举报
文章被收录于专栏:生信技能树

前面我们在教程:使用R包deconstructSigs根据已知的signature进行比例推断,顺利的把508个病人,根据11个signature进行了比例推断,得到的比例矩阵以普通的热图,以及pheatmap包自带的层次聚类如下:

代码是:

代码语言:javascript
复制
rm(list=ls())
options(stringsAsFactors = F) 
load(file = 'mut.wt_from_denovo.Rdata')
a2=mut.wt
## 非负矩阵分解识别驱动性signature并绘图
nmf.input <- t(a2)
nmf.input <- nmf.input[setdiff(rownames(nmf.input),"unknown"),] #去掉unknown
# 需要去除那些没有被计算到的signatures
nmf.input <- nmf.input[rowSums(nmf.input)>0,]
library(pheatmap)
pheatmap(nmf.input)

其中代码里面的escc_denovo_results.Rdata文件,来自于教程:使用R包SomaticSignatures进行denovo的signature推断

可以看到,部分病人的S1,S2比例偏高,很明显,他们应该是一个组,然后剩余的病人可以根据S10进行高低分组,其它signature好像是酱油。

这样的感觉,其实就可以使用NMF算法来实现,尤其是层次聚类并不能很好的把样本进行“泾渭分明”的分组。

第一步:判断最佳分组

就是多次运行 nmf 函数,然后提取cophenetic系数:

代码语言:javascript
复制
library(NMF)
ranks <- 2:10
estim <- lapply(ranks, function(r){
  fit <- nmf(nmf.input, r, nrun = 5, seed = 4, method = "lee") # nrun设置为5以免运行时间过长
  list(fit = fit, consensus = consensus(fit), .opt = "vp",coph = cophcor(fit))
})

names(estim) <- paste('rank', ranks)
sapply(estim, '[[', 'coph')

就可以绘制随rank变化的cophenetic系数图

代码语言:javascript
复制
png("Cophenetic coefficient for seleting optimal nmf rank.png")
par(cex.axis=1.5)
plot(ranks, sapply(estim, '[[', 'coph'), 
     xlab="", ylab="", type="b", 
     col="red", lwd=4,xaxt="n")
axis(side = 1, at=1:10)
title(xlab="number of clusters", ylab="Cophenetic coefficient", cex.lab=1.5)
dev.off()

如下所示:

理论上,根据cophenetic得分选择rank=4,因为它是拐点,但是0元,10小时教学视频直播《跟着百度李彦宏学习肿瘤基因组测序数据分析》 这个文献里面设置为3。这里,我们就尊重文献结果吧。

第二步:筛选signature

前面我们的508个病人,都是11个signature,但是呢,我们的NMF算法运行过后,可以看到有一些signature其实对样本聚类分组并没有意义,所以我们需要挑选一下,代码如下:

代码语言:javascript
复制
rank <- 3
seed <- 2019620
rownames(nmf.input) <- gsub("Signature","Sig",rownames(nmf.input)) # 行名简化
mut.nmf <- nmf(nmf.input, 
               rank = rank, 
               seed = seed, 
               method = "lee") 
index <- extractFeatures(mut.nmf,"max") 
sig.order <- unlist(index)

得到的结果如下;

代码语言:javascript
复制
> index
[[1]]
[1] 2 1

[[2]]
[1]  5  3 11  4

[[3]]
[1] 10  9

attr(,"method")
[1] "max"
> sig.order
[1]  2  1  5  3 11  4 10  9

可以看到,nmf算法得到的3类,都是有各自偏重的signature,11个signature也只有8个需要后续进行分析。如下图,可以看到不同nmf类有各自的偏重的signature。

basismap

第三步:使用挑选出的signature再次NMF

代码语言:javascript
复制
 nmf.input2 <- nmf.input[sig.order,]
library(pheatmap)
pheatmap(nmf.input2,cluster_rows = T,cluster_cols = F)

mut.nmf2 <- nmf(nmf.input2, 
                rank = rank, 
                seed = seed, 
                method = "lee") 
group <- predict(mut.nmf2) # 提出亚型
table(group)

这个时候,病人的分类就是最后的nmf区分成为的3类。泾渭分明,如下:

consensusmap

番外:一些可视化函数

主要是继续参考每个nmf类里面的不同signature的比例,已经不同nmf类的相关性热图

代码语言:javascript
复制
sample.order <- names(group[order(group)])
#设置颜色
jco <- c("#2874C5","#EABF00","#C6524A","#868686")

png(file = "consensusmap.png")
consensusmap(mut.nmf2,
             labRow = NA,
             labCol = NA,
             annCol = data.frame("cluster"=group[colnames(nmf.input)]),
             annColors = list(cluster=c("1"=jco[1],"2"=jco[2],"3"=jco[3],"4"=jco[4])))
dev.off()

png(file = "basismap.png" ) 
# 从这张图可以比较清晰地看到各亚型中的驱动signature(深色),与下面的nmf heatmap对应

basismap(mut.nmf2,
         cexCol = 1,
         cexRow = 0.3,
         annColors=list(c("1"=jco[1],"2"=jco[2],"3"=jco[3],"4"=jco[4])))
dev.off()

aheatmap(as.matrix(nmf.input2[,sample.order]), 
         Rowv=NA, 
         Colv=NA, 
         annCol = data.frame("cluster"=group[sample.order]),
         annColors = list(cluster=c("1"=jco[1],"2"=jco[2],"3"=jco[3],"4"=jco[4])),
         color=c("#EAF0FA","#6081C3","#3454A7"), # 例文的蓝色渐变
         revC=TRUE, 
         cexCol = 0.3, 
         cexRow = 0.3,
         filename = "NMF_heatmap.pdf")

ac=data.frame(group=paste0('NMF_',group))
rownames(ac)=colnames(nmf.input2)
pheatmap(nmf.input2,show_colnames = F,annotation_col = ac)
pheatmap(nmf.input2[,sample.order],
         show_colnames = F,cluster_cols = F,annotation_col = ac)

pheatmap(nmf.input2,show_colnames = F,annotation_col = ac,
         filename = "NMF_heatmap1.pdf")
pheatmap(nmf.input2[,sample.order],
         show_colnames = F,cluster_cols = F,annotation_col = ac,
         filename = "NMF_heatmap2.pdf")
dev.off()

可视化是一条永无止境的路。大家加油哈!

本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2020-06-07,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 生信技能树 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 第一步:判断最佳分组
  • 第二步:筛选signature
  • 第三步:使用挑选出的signature再次NMF
  • 番外:一些可视化函数
相关产品与服务
云直播
云直播(Cloud Streaming Services,CSS)为您提供极速、稳定、专业的云端直播处理服务,根据业务的不同直播场景需求,云直播提供了标准直播、快直播、云导播台三种服务,分别针对大规模实时观看、超低延时直播、便捷云端导播的场景,配合腾讯云视立方·直播 SDK,为您提供一站式的音视频直播解决方案。
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档