0.NMF
NMF是非负矩阵分解，也是常用的划分亚组（亚型）的聚类方法之一咯。
代码主要参考了帮助文档和曾老板的博客：http://www.bio-info-trainee.com/7046.html
1.输入数据
这里使用的还是KIRC的fpkm数据和临床信息，下载自xena。
rm(list=ls())
dat = read.table('TCGA-KIRC.htseq_fpkm.tsv.gz',row.names = 1,header = T,check.names = F)
# 样本筛选，只要tumor
library(tinyarray)
exp = dat[,make_tcga_group(dat)=="tumor"]
# 生存信息
meta = read.table("TCGA-KIRC.survival.tsv",header = T,row.names = 1)
k2 = meta$OS.time>=30;table(k2)
## k2
## FALSE TRUE
## 20 959
meta = meta[k2,]
patient = intersect(colnames(exp),rownames(meta))
exp = exp[,patient]
meta = meta[patient,]
exp = as.matrix(exp)
identical(rownames(meta),colnames(exp))
## [1] TRUE
colnames(meta)[c(1,3)] = c("event","time")
2.用于聚类的基因
tableS3.txt是这篇文章的附表三，里面是免疫相关的基因，有2006个。
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8377979/#MOESM3
为了缩小一下基因数量，我选了免疫相关的基因，并做了批量单因素cox，是tinyarray里的函数。
# 基因筛选
exp = exp[apply(exp, 1, function(x)sum(x==0) < 0.25 *ncol(exp)),]
exp = trans_exp(exp)
tmp = read.table("tables3.txt",header = T)[,1]
exp = exp[rownames(exp)%in%tmp,]
#批量单因素cox
cs = rownames(surv_cox(exp,meta,pvalue_cutoff = 0.001))
length(cs)
## [1] 313
nmf.input <- exp[cs,]
dim(nmf.input)
## [1] 313 518
这两步是为了把基因数量缩小。不是必须这样做哦，单因素cox的p值阈值比较小是因为我不想用太多的基因，基因数量和阈值都没有标准答案哦。
3.完成NMF聚类
也是需要先确定分成几个亚组比较合适。最明确的一个标准就是cophenetic 曲线下降最大的前点。这个操作运行耗时较长，所以我设置了存在即跳过。
library(NMF)
ranks <- 2:6
f = "rank.rdata"
seed <- 2019
if(!file.exists(f)){
result = nmf(nmf.input,ranks,seed =seed)
save(result,file = f)
}
load(f)
plot(result)
选择聚成三类，因为3-4就是最大的下降（第一张图）
rank = 3
result2 <- nmf(nmf.input,
rank = rank,
seed = seed)
index <- extractFeatures(result2,"max") # 能给出选了哪些关键的基因，可以用这些基因再次聚类
nmf.input2 = nmf.input[unlist(index),]
result2 <- nmf(nmf.input2,
rank = rank,
seed = seed)
group <- predict(result2) # 提出亚型
table(group)
## group
## 1 2 3
## 124 54 340
4.热图
这个是NMF R包自带的函数，是把每一个亚组画在了一列，组间的差异非常直观。
basismap(result2,
cexCol = 1.5,
cexRow = 1,
annColors=list(c("#2874C5","#EABF00","#C6524A","#868686")))
如果是常规热图的话，也能看出一定的规律，不如上面明显
dp = nmf.input2[,order(group)]
draw_heatmap(dp,sort(group),
color_an = c("#2874C5","#EABF00","#C6524A","#868686"),
annotation_legend = T,
cluster_cols = F,
show_rownames = T)
5.KM-plot
library(tidyverse)
table(group)
## group
## 1 2 3
## 124 54 340
identical(names(group),rownames(meta))
## [1] TRUE
meta$group = group
library(survival)
library(survminer)
sfit <- survfit(Surv(time, event) ~ group,
data = meta)
ggsurvplot(sfit,pval = T,palette = "jco")
6.PCA和t-SNE可视化
draw_pca(nmf.input,group)
library(Rtsne)
tsne_out = Rtsne(t(nmf.input),perplexity = 30)
pdat = data.frame(tsne_out$Y,factor(group))
colnames(pdat) = c("Y1","Y2","group")
head(pdat)
## Y1 Y2 group
## TCGA-B2-4101-01A -5.444515 3.113312 1
## TCGA-BP-4342-01A 11.238216 1.933829 3
## TCGA-B0-4691-01A 15.687236 6.811072 3
## TCGA-BP-4167-01A 12.930816 3.143940 2
## TCGA-B8-4620-01A 10.009050 5.663727 3
## TCGA-BP-4769-01A -10.508895 13.678469 1
library(ggplot2)
library(paletteer)
ggplot(pdat,aes(Y1,Y2))+
geom_point(aes(Y1,Y2,fill = group),shape = 21,color = "black")+
stat_ellipse(aes(color = group,fill = group),
geom = "polygon",
alpha = 0.3,
linetype = 2)+
scale_color_paletteer_d("RColorBrewer::Set3")+
scale_fill_paletteer_d("RColorBrewer::Set3")+
theme_classic()+
theme(legend.position = "top")
看起来好像比较简单，中间做的时候发现了两个问题：
1.有时extractFeatures提取关键基因，会有某一两个亚组对应的基因是NA，大概是因为这个亚组确实没有什么拿得出手的基因吧。
2.聚类挺好，但KM-plot或者PCA、t-SNE的样本聚类图不能很好的分开。这也是正常的，前面若干步骤选出的用于聚类的基因不同，结果也会大相径庭。
不是所有的问题都有标准答案,想办法在合理范畴调整前面的基因和参数，是完全可以的。
小尾巴
看起来好像比较简单，中间做的时候发现了两个问题：
1.有时extractFeatures提取关键基因，会有某一两个亚组对应的基因是NA，大概是因为这个亚组确实没有什么拿得出手的基因吧。
2.聚类挺好，但KM-plot或者PCA、t-SNE的样本聚类图不能很好的分开。这也是正常的，前面若干步骤选出的用于聚类的基因不同，结果也会大相径庭。
不是所有的问题都有标准答案,调整前面的基因和参数，让结果变得更好，是完全可以的。
