肿瘤分型分析是生信文章中的常客，大致是通过将基因的表达量进行聚类或者非负矩阵分解，发现新的亚型，然后对不同亚型的临床特征，免疫特征等进行比较分析，文章末尾简单的列了一些应用。
本文简答的大概介绍一下文献常用的一致性聚类（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])))
这样就完成了2种常见分子分型方法方法，后续就可以做各种文献中的结合分析。
1）输入数据的基因可以是某个家族的基因，某个通路的基因，某个预后模型中的基因，hub基因等
2）得到分子分型后，可以对不同亚型的临床特征，病理分期，生存状态，免疫特征（RNAseq|免疫浸润也杀疯了，cibersoert？xCELL？ESTIMATE？你常用哪一个）等进行比较分析
3）可以进行差异分析，得到差异基因后可以批量进行单因素生存分析R|生存分析-结果整理
4）分型可以做生存分析以及KM可视化R|生存分析 - KM曲线 ，必须拥有姓名和颜值
