setwd("G:\\肾细胞癌\\肾细胞癌")
cli=read.table('STAD_tpm_for_cibersort.txt', header=T, sep="\t", check.names=F, row.names=1)      #读取临床文件
exprSet_acrg=cli
# 构建GPRscore
# 选择基因
gene=read.table('uniCox_lasso_gpr.txt',header = T)
gene=gene$gene
data=exprSet_acrg[gene,]

# 读取系数
coef=read.csv('GPR_coef.csv',header = T,check.names = F)
coef=coef$Coef..boot_SD
gpr_score=c()
for (i in 1:ncol(data)) {
  score=sum(as.numeric(data[,i])*coef)
  gpr_score=c(gpr_score,score)
}

data=as.data.frame(t(data))
data$gpr_score=gpr_score

#读取生存数据
suv=read.table('time.txt',row.names = 1,header = T,check.names = F)
cli=dplyr::select(suv,'futime','fustat')

##数据合并并输出结果
sameSample=intersect(row.names(data),row.names(cli))
data=data[sameSample,]
cli=cli[sameSample,]

rt=cbind(cli,data)

library(survival)
library(survminer)

### 中位值划分
Type=ifelse(data[,'gpr_score']<= median(rt$gpr_score), "Low", "High")
data=rt
data$group=Type
data$group=factor(data$group, levels=c("Low", "High"))
data_gpr=data
setwd("G:\\肾细胞癌\\肾细胞癌")
# CIBERSORT-Results-GEO.txt请根据预习内容自己评估好！
immune=read.table('CIBERSORT-Results.txt',sep = '\t',header = T,check.names = F,row.names = 1)
immune=immune[,-c(23:25)]


immune=as.data.frame(t(immune))

data=as.data.frame(t(immune))

data=data[,c('B cells naive','Dendritic cells activated', 'Monocytes')]
data=as.data.frame(t(data))

coef=read.csv('TME_coef.csv',header = T,check.names = F)
coef=coef$Coef..boot_SD

TME_score=c()
for (i in 1:ncol(data)) {
  score=sum(as.numeric(data[,i])*coef)
  TME_score=c(TME_score,score)
}

data=as.data.frame(t(data))
data$TME_score=TME_score

#读取生存数据
#读取生存数据
suv=read.table('time.txt',row.names = 1,header = T,check.names = F)
cli=dplyr::select(suv,'futime','fustat')

##数据合并并输出结果
sameSample=intersect(row.names(data),row.names(cli))
data=data[sameSample,]
cli=cli[sameSample,]


rt=cbind(cli,data)
rt$futime=rt$futime

library(survival)
library(survminer)

### 中位值划分，注意此处TME_score负值越大，说明TME浸润越多！！！
### 所以我们认为TME_score越大，其实浸润越low
Type=ifelse(rt[,'TME_score'] <= median(rt$TME_score), "High", "Low")
data=rt
data$group=Type
data$group=factor(data$group, levels=c("Low", "High"))
data_immune=data



data_gpr$GPR_group=ifelse(data_gpr$group=='High','GPR_high','GPR_low')
data_gpr$TME_group=ifelse(data_immune$group=='High','TME_high','TME_low')

data=data_gpr
data$group=paste0(data$GPR_group,'+',data$TME_group)
data$group=stringr::str_replace(data$group,pattern = 'GPR_high\\+TME_high',replacement = 'Mixed')
data$group=stringr::str_replace(data$group,pattern = 'GPR_low\\+TME_low',replacement ='Mixed')

data$group

diff=survdiff(Surv(futime, fustat) ~ group, data = data)
length=length(levels(factor(data[,"group"])))
pValue=1-pchisq(diff$chisq, df=length-1)
if(pValue<0.001){
  pValue="p<0.001"
}else{
  pValue=paste0("p=",sprintf("%.03f",pValue))
}
fit <- survfit(Surv(futime, fustat) ~ group, data = data)

bioCol <- c("#FFA500", "#8B4513", "#1E90FF")

p <- ggsurvplot(fit, 
                data = data,
                conf.int = FALSE,
                pval = pValue,
                pval.size = 6,
                legend.title = 'GPR-TME score',
                legend.labs = levels(factor(data[,"group"])),
                legend = "right",  # 将图例放置在右上角
                font.legend = 12,
                xlab = "Time(Months)",
                palette = bioCol,
                surv.median.line = "hv",
                risk.table = TRUE,
                cumevents = FALSE,
                risk.table.height = 0.25)

p

#8B9E5E、#C3B091、#6495ED
p    
library(scales)

# 设置颜色和透明度
library(scales)

# 设置颜色和透明度
library(RColorBrewer)

# 使用颜色盲友好的调色板
bioCol <- c("#FFA500", "#8B4513", "#1E90FF")

p <- ggsurvplot(fit, 
                data = data,
                conf.int = FALSE,
                pval = pValue,
                pval.size = 6,
                legend.title = 'GPR-TME score',
                legend.labs = levels(factor(data[,"group"])),
                legend = c(0.88, 0.9),
                font.legend = 12,
                xlab = "Time(Months)",
                palette = bioCol,
                surv.median.line = "hv",
                risk.table = TRUE,
                cumevents = FALSE,
                risk.table.height = 0.25)

p



