随着芯片和测序水平的发展，使得我们研究所有基因在整个基因组里的表达情况成为了可能。合理地利用和解释这些数据，能够帮助我们探索相关的生物过程和人类疾病的机制。目前已经有一些软件或方法，可以将具有相似表达模式的基因或者样本进行聚类，但是都有自身的限制。NMF包基于非负矩阵分解(non-negative matrix factorization，以下简称NMF)方法，提取基因表达矩阵内数据的生物相关系数，通过对基因和样本进行组织，抓住数据的内部结构特征，从而对样本进行分组，目前在疾病分型方面受到广泛应用。我前面已经介绍过了NMF的基本原理【NMF(非负矩阵分解)的算法原理】，这里我介绍R语言实现NMF。下面是一篇今年刚发的一篇纯生信的分析文章，用的就是NMF这个方法来对肿瘤进行分型。影响因子为4.8。
文章大家自己去下载：Development and Clinical Validation of a Novel 4-Gene Prognostic Signature Predicting Survival in Colorectal Cancer。
这篇文章做的是结直肠癌。用了2个数据集，一个来自TCGA-COAD，一个GEO数据库：GSE39582，我们就以这篇文章来介绍NMF。
1.NMF包介绍
首先得安装包NMF，这个包属于生物导体包，普通安装方法就行。
代码语言：javascript复制install.packages("NMF") # 安装包的命令
library(NMF) # 加NMF包用法：
代码语言：javascript复制nmf(x, rank,
method, seed = nmf.getOption("default.seed"),
rng = NULL, nrun = if (length(rank) > 1) 30 else 1,
model = NULL, .options = list(),
.pbackend = nmf.getOption("pbackend"),
.callback = NULL, ...)接下来我们看下nmf函数的主要参数：
x：就是我们的表达矩阵； rank：因式分解秩的说明。它通常是一个单一的数值，但也可能是其他类型的值(例如矩阵)，为其实现特定的方法。你可以理解成分几群。
method:就是NMF算法的选择，总共有11种，可通过nmfAlgorithm()函数查看，默认brunet。
关键算法的如下：
nrun：范围内每个值的迭代次数；
model：如果有已经构建好的模型，那就是直接把模型写在这里，优化等级计算。构建模型的函数是nmfModel（rank,c(features,samples)）或者是nmfModel（rank,data,W,H）。 .options:可以设置是否保留每次的运算结果：keep.all=T。例：.options=list(keep.all=TRUE);
2.加载RNAseq数据
下面是该文章的TCGA数据处理方式。
TCGA-COAD的RNAseq数据和临床数据可以从UCSC下载，可参考文章【UCSC数据库下载TCGA数据需要注意的细节】，也可以直接从TCGA数据库下载，自己整理，我之前也把处理过的数据处理过【TCGA数据库33个Project的RNA-Seq转录组数据为你整理打包好了】。
好了，开始...
代码语言：javascript复制rm(list = ls())
options(stringsAsFactors = F)
setwd("K:/MedBioInfoCloud/NMF")###---------加载从UCSC下载的RNA-Seq-FPKM数据
ucsc.coad.fpkm <- read.table("TCGA-COAD.htseq_fpkm.tsv",header = T,check.names = F)
ucsc.coad.fpkm$Ensembl_ID <- substr(ucsc.coad.fpkm$Ensembl_ID,1,15)
load("K:/BioInfoFiles/hsaEnsembl2Symbol.Rdata")
coad.fpkm <- merge(hsaEnsembl2Symbol[,1:2],ucsc.coad.fpkm,by.x ="Ensembl",by.y = "Ensembl_ID")
coad.fpkm <- coad.fpkm[!duplicated(coad.fpkm$gene_name),]
rownames(coad.fpkm) <- coad.fpkm$gene_name
coad.fpkm <- coad.fpkm[,-c(1:2)]
我这里用的是从UCSC下载的数据，从UCSC下载的数据需要自己注释。hsaEnsembl2Symbol这个数据框是我自己注释过的人Ensembl和Symbol之间的对应关系文件。自己注释，参考文章【生信中各种ID转换】。
3.读入临床数据
这里的临床数据我是从UCSC下载的【UCSC数据库下载TCGA数据需要注意的细节】，你可以自己去下载，我也会上传到文章【TCGA数据库33个Project的RNA-Seq转录组数据为你整理打包好了】中对应的网盘中。
代码语言：javascript复制###---------读入临床数据
clin <- read.table("TCGA-COAD.survival.tsv",header = T)
coad.clin <- clin[grep("-01A$",clin$sample),]##删除正常样本
coad.clin <- coad.clin[coad.clin$OS.time>=30,]
coad.clin <- coad.clin[!duplicated(coad.clin$sample),]#去重4.读入基因文件
上面得到的数据是TCGA-COAD所有肿瘤所有基因的FPKM数据，这篇文章做的是能量代谢相关通路的基因，基因信息如下图，自己可以去Reactome数据库搜集。
当然，我这里也准备好了，直接读入数据。
代码语言：javascript复制protein <- read.csv(file = "enengy_protein.csv",header = T)由于搜集数据的时间不同，我搜集的基因是593个，文章是594个，差不多。
代码语言：javascript复制> dim(protein)
[1] 593 3
5.融合数据
将数据与表达矩阵融合。其实也就是只要能量代谢相关基因的表达数据。
代码语言：javascript复制###---------融合数据
gene <- intersect(protein$SYMBOL,rownames(coad.fpkm))
sample <- intersect(coad.clin$sample,colnames(coad.fpkm))
enengyTurExp <- coad.fpkm[gene,sample]
dim(enengyTurExp)
##Genes with FPKM of 0 in half of the samples were removed.
library(tidyr)
genetable <- lapply(as.data.frame(t(enengyTurExp)),function(r){
sum(r==0)}) %>% as.data.frame()
enengyTurExp <- enengyTurExp[colnames(genetable)[genetable[1,] < length(colnames(enengyTurExp))/2],]
sum(is.na(enengyTurExp))##看看有没有空值先看一下表达矩阵的维度
代码语言：javascript复制> dim(enengyTurExp)
[1] 540 421540个基因、421个样本被纳入分析，文章中是371样本。怎么作者少了那么多，我也不知道。
6.NMF进行肿瘤分型
这里我们可以看文章的处理方式。方法用的是brunet，迭代次数=50，ranks为2到10。
因为我们不知道分几个亚群合适，所以设置ranks为2到10，以找到合适的ranks。
代码语言：javascript复制library(NMF)
coad.log2fpkm.enengy <- enengyTurExp##
ranks <- 2:10
estim.coad <- nmf(coad.log2fpkm.enengy,ranks, nrun=50)
duplicated(colnames(coad.log2fpkm.enengy))选择合适的rank。
代码语言：javascript复制#Estimation of the rank: Quality measures computed from 10 runs for each value of r.
plot(estim.coad)
具体怎么选择rank，官方文档是这样建议的：
(Brunet2004) proposed to take the first value of r for which the cophenetic coefficient starts decreasing, (Hutchins2008) suggested to choose the first value where the RSS curve presents aninflection point, and (Frigyesi2008) considered the smallest value at which the decrease in the RSS is lower than the decrease of the RSS obtained from random data.
文章用了三个指标：cophenetic, dispersion 和silhouette。判断最佳rank值的准则就是，cophenetic值随K变化的最大变动的前点，上面结果中cophenetic值在rank为4-5时是第一个变化最大的拐点，所以选择最佳rank值为4。文章中也是4。
代码语言：javascript复制#再次NMF,rank=4
seed = 2020820
nmf.rank4 <- nmf(coad.log2fpkm.enengy,
rank = 4,
nrun=50,
seed = seed,
method = "brunet")绘制分群热图。
代码语言：javascript复制#设置颜色
jco <- c("#2874C5","#EABF00","#C6524A","#868686")
index <- extractFeatures(nmf.rank4,"max")
sig.order <- unlist(index)
NMF.Exp.rank4 <- coad.log2fpkm.enengy[sig.order,]
NMF.Exp.rank4 <- na.omit(NMF.Exp.rank4) #sig.order有时候会有缺失值
group <- predict(nmf.rank4) # 提出亚型
table(group)
consensusmap(nmf.rank4,
labRow = NA,
labCol = NA,
annCol = data.frame("cluster"=group[colnames(NMF.Exp.rank4)]),
annColors = list(cluster=c("1"=jco[1],"2"=jco[2],"3"=jco[3],"4"=jco[4])))
从结果来看，分4个亚群还是可以的。
文章中只给出了consensus matrix这个图（如下）。
得到分群后，就可以进行下游分分析了，可以参考之前TCGA数据库的相关文章【TCGA】 。
