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

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

本文简答的大概介绍一下文献常用的一致性聚类（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-06AMT-CO2         13.45460         12.86170         13.42982         12.99413MT-CO3         12.67183         11.81769         13.10443         12.62491MT-ND4         13.21440         12.37920         13.10761         12.44728MT-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:10seed <- 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     1write.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])))

这样就完成了2种常见分子分型方法方法，后续就可以做各种文献中的结合分析。

1）输入数据的基因可以是某个家族的基因，某个通路的基因，某个预后模型中的基因，hub基因等

2）得到分子分型后，可以对不同亚型的临床特征，病理分期，生存状态，免疫特征（RNAseq|免疫浸润也杀疯了，cibersoert？xCELL？ESTIMATE？你常用哪一个）等进行比较分析

3）可以进行差异分析，得到差异基因后可以批量进行单因素生存分析R|生存分析-结果整理

4）分型可以做生存分析以及KM可视化R|生存分析 - KM曲线，必须拥有姓名和颜值

◆ ◆ ◆ ◆ ◆

更多精心内容详见：精心整理（含图PLUS版）|R语言生信分析，可视化（R统计，ggplot2绘图，生信图形可视化汇总）

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