RNAseq|组学分型-ConsensusClusterPlus（一致性聚类）， NMF（非负矩阵分解）

生信补给站

发布于 2023-08-25 10:11:14

4.3K2

发布于 2023-08-25 10:11:14

文章被收录于专栏：生信补给站

肿瘤分型分析是生信文章中的常客，大致是通过将基因的表达量进行聚类或者非负矩阵分解，发现新的亚型，然后对不同亚型的临床特征，免疫特征等进行比较分析，文章末尾简单的列了一些应用。

本文简答的大概介绍一下文献常用的一致性聚类（ConsensusClusterPlus ）和非负矩阵分解（NMF ）方法。

一载入R包，数据

使用之前得到的RNAseq.SKCM.RData数据集。

library(tidyverse)
library(openxlsx)
#BiocManager::install("ConsensusClusterPlus")
library(ConsensusClusterPlus)
#install.packages("NMF")   # 安装包的命令
library(NMF) # 加NMF包
#使用之前得到的数据
load("RNAseq.SKCM.RData")

#此处展示，选择较小的数据集
table(substr(names(expr),14,16))

expr2 <- expr %>% 
  select(! ends_with("11A"))

expr2 <- as.matrix(expr2)[1:5000,1:100]
expr2[1:4,1:4]

       TCGA-EE-A2GJ-06A TCGA-EE-A2GI-06A TCGA-WE-A8ZM-06A TCGA-DA-A1IA-06A
MT-CO2         13.45460         12.86170         13.42982         12.99413
MT-CO3         12.67183         11.81769         13.10443         12.62491
MT-ND4         13.21440         12.37920         13.10761         12.44728
MT-CO1         13.14420         12.23376         12.93742         12.69075

注意输入数据需要转为矩阵形式，这里的基因可以是某个家族的基因，某个通路的基因，某个预后模型中的基因，hub基因等等，这不就有机的结合起来了。

二一致性聚类（ConsensusClusterPlus）

一致性聚类是一种无监督聚类方法，可以利用ConsensusClusterPlus R包完成分析，表达量矩阵准备好之后，代码很简单，如下

con <- ConsensusClusterPlus(expr, #矩阵形式
                            maxK=10, #最大聚类簇数量
                            reps=100, #抽取的子样本数量
                            pItem=0.8,#抽样样本的比例 
                            pFeature=1,
                            title="resultstrain", #输出文件夹名字
                            clusterAlg="km", #选择聚类算法
                            distance="euclidean", #指定聚类时使用的距离或相关性类型
                            seed=1234, #中子数
                            plot="png", #输出格式 （pdf可能会比较难打开）
                            writeTable=TRUE)

本示例使用的聚类算法是K-means聚类算法，距离是基于欧氏距离（euclidean)，输出格式为png，结果在resultstrain文件夹中。

代码简单但是难点在于选择最优K值，以下仅供参考。 1，Delta area图展示每个K和K-1相比，CDF 曲线下面积的相对变化，值越大表明该k值下的聚类效果相比k-1的聚类效果的优度提升更明显。可以用来帮助决定最佳的K值。

2，一致性累积分布函数 consensus cumulative distribution function，consensus CDF ，图中展示了不同聚类簇数量k下的CDF分布，CDF图可以用来帮助决定最佳的K值

3，一致性矩阵热图矩阵的数值代表同属一个cluster的可能性，取值范围从0到1，颜色从白色到深蓝色，尽量不选择蓝白参杂的K值。

这里可能会倾向选择2 或者 3。（主观，不供参考）

4，每个患者的分型结果在resultstrain （自定义的名字）文件夹中的resultstrain.k=N.consensusClass.csv文件，N为选择的K数字，注意该文件无表头。

三非负矩阵分解（NMF）

除了Consensus Clustering外‍，non-negative matrix factorization (NMF) consensus cluster也是很多文章经常用来分子分型的方式，使用NMF包的nmf函数即可。1，运行NMF

输入表达量矩阵，在初始不清楚rank选择为多少，可以先设置一个范围

ranks <- 2:10
seed <- 1234

result = nmf(expr2,
          ranks,
          method="brunet", 
          nrun=10,
          seed =1234)
plot(result)

如何确定分成几个亚组最合适呢？常用的一个标准就是cophenetic 曲线下降范围最大的前点。由左一图发现4-5下降最大，选择K=4 。确定后，再次进行NMF分析

result2 <- nmf(expr2,
               rank = 4,
               seed = 1234)
index <- extractFeatures(result2,"max") # 提取关键基因
group <- predict(result2) # 提出亚型
table(group)

group
 1  2  3  4 
16 36 28 20

最重要的是保存每个患者的NMF的group结果

NMF_result <- as.data.frame(group)
head(NMF_result,2)
#                group
#TCGA-EE-A2GJ-06A     3
#TCGA-EE-A2GI-06A     1
write.csv(NMF_result,"NMF_result.csv")

2, 可视化 针对NMF结果，一种评估基于指定rank评估聚类稳定性的方法是考虑由多个独立NMF运行结果计算得到的连接矩阵，可以使用consensusmap函数进行绘制。

consensusmap(result2,
             labRow  = NA,
             labCol = NA,
             annCol = data.frame("cluster"=group[colnames(expr2)]))

可以自定义颜色

jco <- c("#2874C5","#EABF00","#C6524A","#868686")
consensusmap(result2,
             labRow  = NA,
             labCol = NA,
             annCol = data.frame("cluster"=group[colnames(expr2)]),
             annColors = list(cluster=c("1"=jco[1],"2"=jco[2],"3"=jco[3],"4"=jco[4])))