
rm(list = ls())
library(tidyverse)

## plot
library(tidyverse)
library(corrplot)
library(ggcor)
library(ggnewscale)

##
hrs = read.csv("0105 Glycolysis.lasson.csv")
tme = read.csv("00202 TME.lasson.csv")
tme = tme[,-c(2,3,4)]

data = dplyr::full_join(hrs, tme, by = "ID")
colnames(data)

##  hypoxia survival plot
library(survival)
library(survminer)
library(ggthemes)

surv = dplyr::filter(data, Cohort %in% c("TCGA_GBM", "TCGA_LGG"))
surv$GRS <- ifelse(surv$GRSscore > median(surv$GRSscore),'High', 'Low') 
surv$GRS <- factor(surv$GRS, levels = c("Low", "High"))

# survival fit
fit <- survfit(Surv(OS.time, OS) ~ GRS, data = surv)
p1 <- ggsurvplot(fit,
                 #title = "TCGA.OV",
                 pval = TRUE, 
                 #conf.int = TRUE,
                 risk.table = T,
                 #risk.table.pos='in',
                 risk.table.col = "strata",
                 xlab = "Time (months)",
                 ylab = "Overall Survival",
                 #ylab = "Disease Free Survival (%)",
                 linetype = 1,
                 #surv.median.line = "hv", 
                 ggtheme = theme_article(),
                 legend.labs = c("Low", "High"),
                 palette = c("#00BFC4", "#F8766D"),
                 legend = c(0.78,0.90),
                 legend.title = "GRS score",
                 xlim = c(0,150),
                 #break.x.by = 24,
                 surv.scale = "percent")
p1


##  TME survival plot
# surv <- dplyr::filter(data, GSE=="TCGA.OV")
surv$TME <- ifelse(surv$TMEscore > median(surv$TMEscore),'Low', 'High') 
surv$TME <- factor(surv$TME, levels = c("Low", "High"))

# survival fit
fit <- survfit(Surv(OS.time, OS) ~ TME, data = surv)
p2 <- ggsurvplot(fit,
                 #title = "All HGSOC patients",
                 pval = TRUE, 
                 #conf.int = TRUE,
                 risk.table = T,
                 #risk.table.pos='in',
                 risk.table.col = "strata",
                 xlab = "Time (months)",
                 ylab = "Overall Survival",
                 #ylab = "Disease Free Survival (%)",
                 linetype = 1,
                 #surv.median.line = "hv", 
                 ggtheme = theme_article(),
                 legend.labs = c("Low", "High"),
                 palette = c("#00BFC4", "#F8766D"),
                 legend = c(0.78,0.90),
                 legend.title = "TME score",
                 xlim = c(0,150),
                 #break.x.by = 24,
                 surv.scale = "percent")
p2


## classier
# surv <- dplyr::filter(data, GSE=="TCGA.OV")
surv$GRS <- ifelse(surv$GRSscore > median(surv$GRSscore),'GRSHigh', 'GRSLow')
surv$TME <- ifelse(surv$TMEscore > median(surv$TMEscore),'TMELow', 'TMEHigh') 

surv$group = paste(surv$GRS, surv$TME, sep = "/")
table(surv$group)

fit <- survfit(Surv(OS.time, OS) ~ group, data = surv)
p3 <- ggsurvplot(fit,
                 title = "",
                 pval = TRUE, 
                 #conf.int = TRUE,
                 risk.table = F, 
                 risk.table.col = "strata",
                 xlab = "Time (months)",
                 ylab = "Overall Survival",
                 #ylab = "Disease Free Survival (%)",
                 linetype = "strata",
                 #surv.median.line = "hv", 
                 ggtheme = theme_article(),
                 #legend.labs = c("Low TME", "High TME"),
                 palette = c("#C77CFF", "#F8766D", "#00BFC4", "#7CAE00"),
                 legend = c(0.45,0.22),
                 legend.title = "",
                 xlim = c(0,150),
                 #break.x.by = 24,
                 surv.scale = "percent")
p3



# surv <- dplyr::filter(data, GSE=="TCGA.OV")
surv$GRS <- ifelse(surv$GRSscore > median(surv$GRSscore),'GRSHigh', 'GRSLow')
surv$TME <- ifelse(surv$TMEscore > median(surv$TMEscore),'TMELow', 'TMEHigh') 

surv$group = paste(surv$GRS, surv$TME, sep = "/")
table(surv$group)

surv$GRSTME = forcats::fct_collapse(surv$group, Mixed = c("GRSHigh/TMEHigh", "GRSLow/TMELow"))
surv$GRSTME = factor(surv$GRSTME, 
                     levels = c("GRSLow/TMEHigh", "Mixed", "GRSHigh/TMELow"),
                     labels = c("GRSLow/TMEHigh", "Mixed", "GRSHigh/TMELow"))

fit <- survfit(Surv(OS.time, OS) ~ GRSTME, data = surv)
p4 <- ggsurvplot(fit,
                 title = "",
                 pval = TRUE, 
                 #conf.int = TRUE,
                 risk.table = F, 
                 risk.table.col = "strata",
                 xlab = "Time (months)",
                 ylab = "Overall Survival",
                 #ylab = "Disease Free Survival (%)",
                 linetype = "strata",
                 #surv.median.line = "hv", 
                 ggtheme = theme_article(),
                 #legend.labs = c("Low TME", "High TME"),
                 palette = c("#00BFC4", "#7CAE00", "#F8766D"),
                 legend = c(0.45,0.22),
                 legend.title = "",
                 xlim = c(0,150),
                 #break.x.by = 24,
                 surv.scale = "percent")
p4




## time-ROC
set.seed(2022)

surv = dplyr::filter(data, Cohort %in% c("TCGA_GBM", "TCGA_LGG"))
cox <- coxph(Surv(OS.time, OS) ~ TMEscore + GRSscore, data = surv)
surv$GRSTMEscore=predict(cox, type="risk", newdata=surv)

library(timeROC)
result <-with(surv, timeROC(T=OS.time,
                            delta=OS,
                            marker=GRSTMEscore,
                            cause=1,
                            times=c(12,36,60),
                            iid = TRUE))
plot_dat = data.frame(fpr = as.numeric(result$FP),
                      tpr = as.numeric(result$TP),
                      time = rep(as.factor(c(12,36,60)), each = nrow(result$TP)))

library(ggplot2)
p6 <- ggplot() + 
  geom_line(data = plot_dat, aes(x = fpr, y = tpr, color = time), size = 1) + 
  scale_color_manual(name = NULL,values = c("#00BFC4", "#7CAE00", "#F8766D"),
                     labels = paste0("AUC of ",c(1,3,5)," year: ",
                                     format(round(result$AUC,3),nsmall = 3)))+
  geom_line(aes(x=c(0,1), y=c(0,1)), color = "black")+
  theme_article()+
  theme(axis.text = element_text(color = "black"),
        legend.position = c(0.78,0.125))+
  scale_x_continuous(expand = c(0.005,0.005))+
  scale_y_continuous(expand = c(0.005,0.005))+
  labs(x = "1 - Specificity",
       y = "Sensitivity")
p6


## pca
library(tidyverse)
library(ggplot2)
library(ggsci)
library(Cairo)

colnames(surv)
surv$TMEscore = -1*surv$TMEscore
pca = prcomp(surv[,c("GRSscore", "TMEscore")])

pca.res = pca[["x"]] %>%
  as.data.frame()

pca.res = cbind(pca.res, surv[,"GRSTME"])
colnames(pca.res)[3] = 'group'

pca.var = pca$sdev^2 %>%
  as.data.frame()
pca.var$var = round(pca.var$. / sum(pca.var) * 100, 2) # 计算各主成分所占百分比
pca.var$pc = colnames(pca.res)[1:(ncol(pca.res)-1)]

p5 = ggplot(pca.res, aes(PC1, PC2, color = group, shape = group))+ 
  # 选择X轴Y轴并映射颜色和形状
  geom_point(size = 1.5)+ # 画散点图并设置大小
  geom_hline(yintercept = 0, linetype="dashed") + # 添加横线
  geom_vline(xintercept = 0, linetype="dashed") + # 添加竖线
  #scale_color_manual(c("#00BFC4", "#7CAE00", "#F8766D"))
  
  theme_article() + # 加上边框
  #stat_ellipse(level = 0.95)+ # 添加置信椭圆
  # 自动提取主成分解释度进行绘图
  labs(x = paste('PC1(', pca.var$var[1],'%)', sep = ''),
       y = paste('PC2(', pca.var$var[2],'%)', sep = '')) +
  theme(legend.position = c(0.85,0.85)) # 设置图例位置，此处为相对位置

p5

### Figure 
###################
library(patchwork)
library(ggplotify)
p = (p1$plot|p2$plot|p3)/(p4$plot|p5$plot|p6)
pa = plot_spacer() | p + plot_layout(ncol = 1, heights = c(3,4.5)) + plot_annotation("A")
p

p = (p1$plot|p2$plot|p3$plot)/(p4$plot|p5|p6) 


