
## 载入获得的CIBERSORT
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))

#删掉正常样品
group=sapply(strsplit(colnames(immune),"\\-"),"[",4)
group=sapply(strsplit(group,""),"[",1)
group=gsub("2","1",group)
immune=immune[,group==0]

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

#读取生存数据
suv=read.table('TCGA-KIRC.survival.tsv',row.names = 1,header = T,check.names = F)
cli=dplyr::select(suv,'OS.time','OS')
colnames(cli)=c("futime", "fustat")

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


## K-M生存分析免疫细胞,找预后相关的
rt=cbind(cli,data)
rt$futime=rt$futime/30
## 删去全为0,或基本为0的免疫细胞!!!!!!!
immune1=immune[-c(2,5,20,21),]
rt1=rt[,-c(4,7,22,23)]
immune_p=c()
immune_figure=list()
library(survival)
library(survminer)

#install.packages('cowplot')
library(cowplot)
dir.create('k_m')
for (i in rownames(immune1)) {
  res.cut=surv_cutpoint(rt1, time="futime", event="fustat", variables=i)
  cutoff=as.numeric(res.cut$cutpoint[1])
  print(cutoff)
  Type=ifelse(data[,i]<= cutoff, "Low", "High")
  data=rt1
  data$group=Type
  data$group=factor(data$group, levels=c("Low", "High"))
  diff=survdiff(Surv(futime, fustat) ~ group, data = data)
  length=length(levels(factor(data[,"group"])))
  pValue=1-pchisq(diff$chisq, df=length-1)
  immune_p=c(immune_p,pValue)
  if(pValue<0.05){
    fit <- survfit(Surv(futime, fustat) ~ group, data = data)
    bioCol=c("#0073C2","#EFC000","#6E568C","#7CC767","#223D6C","#D20A13","#FFD121","#088247","#11AA4D")
    bioCol=bioCol[1:length]
    
    p=ggsurvplot(fit, 
                 data=data,
                 conf.int=F,
                 pval=pValue,
                 pval.size=6,
                 legend.title=i,
                 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=F,
                 cumevents=F,
                 risk.table.height=.25)
    ggsave2(filename = paste0('./k_m/',i,'.pdf'),width = 4,height = 4)
  }
}

## 找出三种免疫细胞


library(survival)
library(survminer)
cox=coxph(Surv(futime, fustat) ~  `B cells naive`+`Dendritic cells activated`+ `Monocytes`,data = rt)
ggforest(cox)

############
#install.packages('boot')
library(boot)
rsq <- function(formula, data, indices) { 
  d <- data[indices,] 
  fit <- coxph(formula, data=d) 
  return(fit$coefficients) 
} 

set.seed(123456)
## 1000次重抽样多因素
boot_results <- boot(data=rt, statistic=rsq, 
                     R=1000, formula=Surv(futime, fustat) ~ `B cells naive`+`Dendritic cells activated`+ `Monocytes`)

cox$coefficients

print(boot_results)

## 获取参数
coef=boot_results$t0
sd=as.data.frame(boot_results$t)
sd=apply(sd, 2, sd)

## 定义coef/sd为新的参数
ratio=coef/sd

TME=data.frame('Coef'=coef,'boot_SD'=sd,'Coef\\/boot_SD'=ratio)
write.csv(TME,file= 'TME_coef.csv',quote = F)







###########TME_score

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))

#删掉正常样品
group=sapply(strsplit(colnames(immune),"\\-"),"[",4)
group=sapply(strsplit(group,""),"[",1)
group=gsub("2","1",group)
immune=immune[,group==0]

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('TCGA-KIRC.survival.tsv',row.names = 1,header = T,check.names = F)
cli=dplyr::select(suv,'OS.time','OS')
colnames(cli)=c("futime", "fustat")

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


## K-M生存分析免疫细胞
rt=cbind(cli,data)
rt$futime=rt$futime/30

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"))
diff=survdiff(Surv(futime, fustat) ~ group, data = data)
length=length(levels(factor(data[,"group"])))
pValue=1-pchisq(diff$chisq, df=length-1)

fit <- survfit(Surv(futime, fustat) ~ group, data = data)
bioCol=c("#0073C2","#EFC000","#6E568C","#7CC767","#223D6C","#D20A13","#FFD121","#088247","#11AA4D")
bioCol=bioCol[1:length]
bioCol <- c("#FFA500", "#1E90FF")
p=ggsurvplot(fit, 
             data=data,
             conf.int=F,
             pval=pValue,
             pval.size=6,
             legend.title='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=T,
             cumevents=F,
             risk.table.height=.25)

p    

data_immune=data
save(data_immune,file ='TMEscore_and_group.Rdata')


## 相关性 #########################################
load('KIRC_tpm.Rdata')
gene=read.table('uniCox_lasso_gpr.txt',header = T)
gene=gene$gene
data=exprSet_tcga_mRNA[gene,]

#删掉正常样品
group=sapply(strsplit(colnames(data),"\\-"),"[",4)
group=sapply(strsplit(group,""),"[",1)
group=gsub("2","1",group)
data=data[,group==0]

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


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))

#删掉正常样品
group=sapply(strsplit(colnames(immune),"\\-"),"[",4)
group=sapply(strsplit(group,""),"[",1)
group=gsub("2","1",group)
immune=immune[,group==0]
immune=as.data.frame(t(immune))

identical(rownames(immune),rownames(data))

ss=intersect(rownames(immune),rownames(data))
immune=immune[ss,]
colnames(immune)


immune=immune[,c('B cells naive','Dendritic cells activated', 'Monocytes')]
data=data[ss,]

rt=cbind(data,immune)


#install.packages('corrplot')
#install.packages('ggcorrplot')

library(corrplot)
M <- cor(rt)
res1 <- cor.mtest(rt, conf.level = .95) 
library(ggcorrplot)
ggcorrplot(
  M,
  hc.order = F,
  #type = "lower",
  outline.color = NA,
  ggtheme = ggplot2::theme_gray,
  colors = c("blue", "white", "red")
)


### 联合预后
library(survival)
library(survminer)
load('GPRscore_and_group.Rdata')
load('TMEscore_and_group.Rdata')


data$GPR_group=ifelse(data$group=='High','GPR_high','GPR_low')
data$TME_group=ifelse(data_immune$group=='High','TME_high','TME_low')
data$group=paste0(data$GPR_group,'+',data$TME_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("#0073C2","#EFC000","#6E568C","#7CC767","#223D6C","#D20A13","#FFD121","#088247","#11AA4D")
bioCol=bioCol[1:length]
bioCol <- c("#FFA500", "#1E90FF")
p=ggsurvplot(fit, 
             data=data,
             conf.int=F,
             pval=pValue,
             pval.size=6,
             legend.title='GPR-TME score',
             legend.labs=levels(factor(data[,"group"])),
             legend = c(0.8, 0.8),
             font.legend=12,
             xlab="Time(Months)",
             palette = bioCol,
             surv.median.line = "hv",
             risk.table=T,
             cumevents=F,
             risk.table.height=.25)
p    
table(data$group)
rt_roc=data[data$group %in% c('GPR_high+TME_low','GPR_low+TME_high'),]

rt_roc$class=ifelse(rt_roc$group=='GPR_high+TME_low',1,0)
pROC::plot.roc(rt_roc$fustat,rt_roc$class,print.auc=T)

#install.packages('timeROC')
library(timeROC)
library(survival)
library(survivalROC)

data=rt_roc

time_roc_res <- timeROC(
  T = data$futime,
  delta = data$fustat,
  marker = data$class,
  cause = 1,
  weighting="marginal",
  times = c(3 * 12, 5 * 12,7*12),
  ROC = TRUE,
  iid = TRUE
)


time_ROC_df <- data.frame(
  TP_3year = time_roc_res$TP[, 1],
  FP_3year = time_roc_res$FP[, 1],
  TP_5year = time_roc_res$TP[, 2],
  FP_5year = time_roc_res$FP[, 2],
  TP_7year = time_roc_res$TP[, 3],
  FP_7year = time_roc_res$FP[, 3]
)


bioCol=c("#0073C2","#EFC000","#6E568C","#7CC767","#223D6C","#D20A13","#FFD121","#088247","#11AA4D")

library(ggplot2)
ggplot(data = time_ROC_df) +
  geom_line(aes(x = FP_3year, y = TP_3year), size = 1.3, color = "#0073C2") +
  geom_line(aes(x = FP_5year, y = TP_5year), size = 1.3, color = "#EFC000") +
  geom_line(aes(x = FP_7year, y = TP_7year), size = 1.3, color = "#6E568C") +
  geom_abline(slope = 1, intercept = 0, color = "grey", size = 1, linetype = 2) +
  theme_bw() +
  annotate("text",
           x = 0.7, y = 0.25, size = 3.2,
           label = paste0("AUC at 3 years = ", sprintf("%.3f", time_roc_res$AUC[[1]])), color = "#0073C2"
  ) +
  annotate("text",
           x = 0.7, y = 0.15, size = 3.2,
           label = paste0("AUC at 5 years = ", sprintf("%.3f", time_roc_res$AUC[[2]])), color = "#EFC000"
  ) +
  annotate("text",
           x = 0.7, y = 0.05, size = 3.2,
           label = paste0("AUC at 7 years = ", sprintf("%.3f", time_roc_res$AUC[[3]])), color = "#6E568C"
  ) +
  labs(x = "False positive rate", y = "True positive rate") 



