通过不同的聚类算法，依据不同的基因或指标给样本聚类，是数据挖掘文章里面经常出现的操作咯。接下来要写写不同聚类方法的介绍和代码实现啦。
1.输入数据
目前从网上找到的代码大多是帮助文档的示例数据ALL，是个log过后的表达芯片数据，作者取了mad最大的5000个基因，并且把每个基因的表达量减去了中位数，使数据分布范围在0上下。
芯片数据，就照这个样子处理。那RNA-seq数据呢，是应该用log2后的FPKM、TPM这些，是直接使用还是也需要减一下中位数呢。我查了R包的帮助文档、使用手册、github、文章，都没有提出明确的输入数据要求啊，已发表的文章里面也大多没有给出具体的说明。其实也不只是可以用基因表达量，已发表的文章里也有拿蛋白表达数据去聚类。 那…就自由发挥咯，如果以后有了更新，我会在这里补充出来。
帮助文档的示例代码如下
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数据，来做做这个聚类，顺便把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这个包确实简单，一个函数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”，其他参数的意思可以看帮助文档。(这里提到的大佬就是下面PAC方法的作者)
此时工作目录下的test文件夹里已经有一个consensus.pdf和一个ici.pdf，里面各有多张图了。
结果就不全部贴出来了，结合下面的筛选标准贴一部分。
3.筛选最佳聚类数
也就是公式里的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
针不戳，和上面的指标结果基本统一，最佳聚类数是2。但是呐，这个方法也会失灵，有的时候你maxK设置几，最后结果就返回几。这时就要对照一下前面的指标结果啦,它最佳聚类数返回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")
