结果如图所示:
代码如下:
#引用包
library(limma)
library(tidyverse)
library(ConsensusClusterPlus)
#####导入相关基因##########
gene=c("ACAN",
"ACHE",
"ADAM10",
"ADAM17",
"ADAM9")
########整理输入文件,并根据自己需求更改参数########
####注:导入的输入文件为TCGA的TPM格式的表达数据文件以及整理好的临床数据文件。
expdata=read.table("TCGA_TPM.txt",header = T,sep = "\t")
clinical=read.table("clinical.txt",header = T,row.names = 1,sep = "\t")
clinical$futime=clinical$futime/365
exp_data_T = expdata%>% dplyr::select(str_which(colnames(.), ".01"))
rt=exp_data_T
colnames(rt)=gsub("(.*?)\\.(.*?)\\.(.*?)\\..*", "\\1\\-\\2\\-\\3\\", colnames(rt))#去除.01
sameSample=intersect(rownames(clinical),colnames(rt))#合并共有的样本
clinical_data=clinical[sameSample,][,1:9]#根据你需要的改提取的临床信息
data_ci=t(rt)
data_ci=data_ci[sameSample,]
data_ci=data_ci[,gene]
cli_data=cbind(clinical_data,data_ci)
cli_data=as.data.frame(cli_data)
###colnames(cli_data)=c("futime","fustat",gene)####当只想要生存数据和生存状况数据时启用
data=cli_data[,10:225]#218要根据输入的相关基因更改,10则根据你要的临床数据更改
data=as.data.frame(data)
rownames(data)=rownames(cli_data)
colnames(data)=gene
data=t(data)
####开始进行无监督聚类分析#####
maxK=10 #设置最大的K值,最多可以将样品分成几个组
results=ConsensusClusterPlus(data,
maxK=maxK,#最大聚类数目。
reps=100,#重复次数,默认为100。
pItem=0.8,#每次重取样的样本比例,默认为0.8。
pFeature=1,#每次重取样的特征比例,默认为1。
#title=workDir,
clusterAlg="pam", #聚类的算法
distance="euclidean", #计算距离的方法
seed=2, #设置种子
plot="png") #输出图形的格式
#输出结果
clusterNum=2 #最优的K值
Cluster=results[[clusterNum]][["consensusClass"]]
outTab=cbind(cli_data, Cluster)
outTab[,"Cluster"]=paste0("C", outTab[,"Cluster"])
outTab=rbind(ID=colnames(outTab), outTab)#添加列名到第一列
write.table(outTab, file="cluster.txt", sep="\t", quote=F, col.names=F)