通過不同的聚類算法损话，依據(jù)不同的基因或指標給樣本聚類赛不，是數(shù)據(jù)挖掘文章里面經(jīng)常出現(xiàn)的操作咯窟蓝。接下來要寫寫不同聚類方法的介紹和代碼實現(xiàn)啦摩瞎。
1.輸入數(shù)據(jù)
目前從網(wǎng)上找到的代碼大多是幫助文檔的示例數(shù)據(jù)ALL纹笼，是個log過后的表達芯片數(shù)據(jù)纹份，作者取了mad最大的5000個基因，并且把每個基因的表達量減去了中位數(shù)廷痘，使數(shù)據(jù)分布范圍在0上下蔓涧。
芯片數(shù)據(jù)，就照這個樣子處理笋额。那RNA-seq數(shù)據(jù)呢元暴，是應該用log2后的FPKM、TPM這些鳞陨，是直接使用還是也需要減一下中位數(shù)呢昨寞。我查了R包的幫助文檔、使用手冊厦滤、github援岩、文章，都沒有提出明確的輸入數(shù)據(jù)要求啊掏导，已發(fā)表的文章里面也大多沒有給出具體的說明享怀。其實也不只是可以用基因表達量，已發(fā)表的文章里也有拿蛋白表達數(shù)據(jù)去聚類趟咆。 那…就自由發(fā)揮咯添瓷，如果以后有了更新梅屉，我會在這里補充出來。
幫助文檔的示例代碼如下
rm(list = ls())
library(ALL)
data(ALL)
df <- exprs(ALL)
dim(df)
## [1] 12625 128
df[1:4,1:4]
## 01005 01010 03002 04006
## 1000_at 7.597323 7.479445 7.567593 7.384684
## 1001_at 5.046194 4.932537 4.799294 4.922627
## 1002_f_at 3.900466 4.208155 3.886169 4.206798
## 1003_s_at 5.903856 6.169024 5.860459 6.116890
mads <- apply(df,1,mad)
df <- df[rev(order(mads))[1:5000],]
par(mfrow = c(1,2))
boxplot(df[,1:20],main = "before")
df <- sweep(df,1, apply(df,1,median,na.rm=T))
boxplot(df[,1:20],main = "after")
這里我就采用TCGA的KIRC-fpkm數(shù)據(jù)鳞贷，來做做這個聚類坯汤，順便把KM-plot和PCA、t-SNE樣本聚類做一下搀愧。
TCGA-KIRC.htseq_fpkm.tsv.gz和TCGA-KIRC.survival.tsv是從xena下載的KIRC矩陣和生存信息惰聂。
rm(list = ls())
library(tinyarray)
# 表達矩陣
dat = read.table("TCGA-KIRC.htseq_fpkm.tsv.gz",
header = T,
row.names = 1,
check.names = F)
# 樣本篩選，只要tumor
exp = dat[,make_tcga_group(dat)=="tumor"]
# 基因篩選
exp = exp[apply(exp, 1, function(x)sum(x==0) < 0.5 *ncol(exp)),]
mads = apply(exp, 1, mad)
exp = exp[tail(order(mads),5000),]
# 生存信息
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
df <- sweep(exp,1, apply(exp,1,median,na.rm=T))
2.完成聚類
ConsensusClusterPlus這個包確實簡單咱筛，一個函數(shù)N張圖搓幌，真棒O(∩_∩)O
library(ConsensusClusterPlus)
maxK <- 6 #最多分成幾組
results <- ConsensusClusterPlus(df,
maxK = maxK,
reps = 500,
pItem = 0.8,
pFeature = 1,
clusterAlg = "pam",
seed = 102,
title="test",
innerLinkage="complete",
plot="pdf")
icl = calcICL(results,
title="test",
plot="pdf")
大佬推薦使用innerLinkage=“complete”，不推薦用默認的“average”迅箩，其他參數(shù)的意思可以看幫助文檔溉愁。(這里提到的大佬就是下面PAC方法的作者)
此時工作目錄下的test文件夾里已經(jīng)有一個consensus.pdf和一個ici.pdf，里面各有多張圖了饲趋。
結果就不全部貼出來了拐揭，結合下面的篩選標準貼一部分。
3.篩選最佳聚類數(shù)
也就是公式里的K值篙贸，即最終把這些樣本聚成幾類投队。
關于這個呀。我翻閱了很多資料爵川，確定K值的指標簡直是五花八門的敷鸦，主要標準有：
一致性矩陣熱圖白色塊最干凈，盡量不摻雜藍色
累積分布曲線下降的坡度最平緩
delta area 曲線的肘部點橫坐標
-
聚類一致性直方圖又高又平均
這些指標呢寝贡，又有可能并不指向同一個K值扒披，這時你就開始左右為難啦。
然后就有了大佬提供的PAC標準去篩選圃泡，拒絕選擇困難癥碟案。
Kvec = 2:maxK
x1 = 0.1; x2 = 0.9 # threshold defining the intermediate sub-interval
PAC = rep(NA,length(Kvec))
names(PAC) = paste("K=",Kvec,sep="") # from 2 to maxK
for(i in Kvec){
M = results[[i]]$consensusMatrix
Fn = ecdf(M[lower.tri(M)])
PAC[i-1] = Fn(x2) - Fn(x1)
}#end for i
# The optimal K
optK = Kvec[which.min(PAC)]
optK
## [1] 2
PAC代碼出自：https://www.linkedin.com/pulse/how-use-pac-measure-consensus-clustering-yasin-%C5%9Fenbabao%C4%9Flu
針不戳，和上面的指標結果基本統(tǒng)一颇蜡，最佳聚類數(shù)是2价说。但是吶，這個方法也會失靈风秤，有的時候你maxK設置幾鳖目，最后結果就返回幾。這時就要對照一下前面的指標結果啦,它最佳聚類數(shù)返回10缤弦，你不能真聚成10類吧领迈。。。
4.亞型KM-plot
#聚類結果
library(tidyverse)
table(results[[optK]]$consensusClass)
##
## 1 2
## 283 235
Cluster = results[[optK]]$consensusClass
identical(names(Cluster),rownames(meta))
## [1] TRUE
meta$Cluster = Cluster
library(survival)
library(survminer)
sfit <- survfit(Surv(OS.time, OS) ~ Cluster,
data = meta)
ggsurvplot(sfit,pval = T,palette = "jco")
5.PCA和t-SNE可視化
就是做一下樣本聚類狸捅，看看是否聚在一起衷蜓。
draw_pca(exp,Cluster,addEllipses = F)
library(Rtsne)
tsne_out = Rtsne(t(exp),perplexity = 30)
pdat = data.frame(tsne_out$Y,factor(Cluster))
colnames(pdat) = c("Y1","Y2","group")
head(pdat)
## Y1 Y2 group
## TCGA-B2-4101-01A 8.9694069 6.436323 1
## TCGA-BP-4342-01A 3.3333001 -2.952281 1
## TCGA-B0-4691-01A 4.5284118 -10.801060 1
## TCGA-BP-4167-01A 0.7638295 -10.804211 1
## TCGA-B8-4620-01A 1.3455857 -2.727872 1
## TCGA-BP-4769-01A -21.0383293 1.704198 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")
