微信：基于R语言利用NMF(非负矩阵分解)替代层次聚类进行肿瘤分型
层次聚类最后是将样本根据表达量进行了分类，可以用于样本的分型，但是不能用于基因分cluster。

#加载文件
rt=read.table("normal.txt",sep="\t",header=T,check.names=F)
#基因去重
library(limma)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data=avereps(data)
data<-as.data.frame(data)
nrow(data)
data=data[rowMeans(data)>0.5,]
#行是基因，列是样本
dim(data) #11890个基因，6个样本

##数据读取完成，选择方差大于1/4的基因
m.vars=apply(data,1,var)
data.upper=data[which(m.vars>quantile(m.vars, probs = seq(0, 1,0.25))[4]),]
dim(data.upper)
write.table(data.upper,file="data_variancetop0.25.txt",sep='\t',quote=F,row.names=T)


#NMF方法用的是brunet，迭代次数=50，ranks为2到10。
library(NMF)
coad.log2fpkm.enengy <- data.upper
ranks <- 2:10
#rank就是分几群，nrun是迭代次数，method默认就是brunet
estim.coad <- nmf(coad.log2fpkm.enengy,ranks, nrun=50,.options=list(keep.all=TRUE))
duplicated(colnames(coad.log2fpkm.enengy))
#选择合适的rank,文章用了三个指标：cophenetic, dispersion 和silhouette
#判断最佳rank值的准则就是，cophenetic值随K变化的最大变动的前点，
#上面结果中cophenetic值在rank为3时是第一个变化最大的拐点，
plot(estim.coad)


#再次NMF,rank=3
seed = 2020820
nmf.rank4 <- nmf(coad.log2fpkm.enengy, 
                 rank = 3, 
                 nrun=50,
                 seed = seed, 
                 method = "brunet")

#绘制分群热图
#设置颜色
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) # 提出亚型

group#这是把样本进行分群
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])))

data