0.NMF
NMF是非負(fù)矩陣分解，也是常用的劃分亞組（亞型）的聚類方法之一咯寻仗。
代碼主要參考了幫助文檔和曾老板的博客：http://www.bio-info-trainee.com/7046.html
1.輸入數(shù)據(jù)
這里使用的還是KIRC的fpkm數(shù)據(jù)和臨床信息匈睁，下載自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是這篇文章的附表三综看，里面是免疫相關(guān)的基因脾拆，有2006個(gè)馒索。
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8377979/#MOESM3
為了縮小一下基因數(shù)量莹妒，我選了免疫相關(guān)的基因，并做了批量單因素cox绰上，是tinyarray里的函數(shù)旨怠。
# 基因篩選
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
這兩步是為了把基因數(shù)量縮小。不是必須這樣做哦渔期，單因素cox的p值閾值比較小是因?yàn)槲也幌胗锰嗟幕蛟讼牛驍?shù)量和閾值都沒有標(biāo)準(zhǔn)答案哦。
3.完成NMF聚類
也是需要先確定分成幾個(gè)亞組比較合適疯趟。最明確的一個(gè)標(biāo)準(zhǔn)就是cophenetic 曲線下降最大的前點(diǎn)拘哨。這個(gè)操作運(yùn)行耗時(shí)較長，所以我設(shè)置了存在即跳過信峻。
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)
選擇聚成三類倦青，因?yàn)?-4就是最大的下降（第一張圖）
rank = 3
result2 <- nmf(nmf.input,
rank = rank,
seed = seed)
index <- extractFeatures(result2,"max") # 能給出選了哪些關(guān)鍵的基因，可以用這些基因再次聚類
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.熱圖
這個(gè)是NMF R包自帶的函數(shù)盹舞，是把每一個(gè)亞組畫在了一列产镐，組間的差異非常直觀。
basismap(result2,
cexCol = 1.5,
cexRow = 1,
annColors=list(c("#2874C5","#EABF00","#C6524A","#868686")))
如果是常規(guī)熱圖的話踢步，也能看出一定的規(guī)律癣亚，不如上面明顯
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")
看起來好像比較簡單，中間做的時(shí)候發(fā)現(xiàn)了兩個(gè)問題：
1.有時(shí)extractFeatures提取關(guān)鍵基因获印，會(huì)有某一兩個(gè)亞組對(duì)應(yīng)的基因是NA述雾，大概是因?yàn)檫@個(gè)亞組確實(shí)沒有什么拿得出手的基因吧。
2.聚類挺好兼丰，但KM-plot或者PCA玻孟、t-SNE的樣本聚類圖不能很好的分開。這也是正常的鳍征，前面若干步驟選出的用于聚類的基因不同黍翎，結(jié)果也會(huì)大相徑庭。
不是所有的問題都有標(biāo)準(zhǔn)答案,想辦法在合理范疇調(diào)整前面的基因和參數(shù)艳丛，是完全可以的匣掸。
小尾巴
看起來好像比較簡單，中間做的時(shí)候發(fā)現(xiàn)了兩個(gè)問題：
1.有時(shí)extractFeatures提取關(guān)鍵基因氮双，會(huì)有某一兩個(gè)亞組對(duì)應(yīng)的基因是NA旺聚，大概是因?yàn)檫@個(gè)亞組確實(shí)沒有什么拿得出手的基因吧。
2.聚類挺好眶蕉，但KM-plot或者PCA、t-SNE的樣本聚類圖不能很好的分開唧躲。這也是正常的造挽，前面若干步驟選出的用于聚類的基因不同碱璃，結(jié)果也會(huì)大相徑庭。
不是所有的問題都有標(biāo)準(zhǔn)答案,調(diào)整前面的基因和參數(shù)饭入，讓結(jié)果變得更好嵌器，是完全可以的。
