
rm(list = ls())

library(tidyverse)
library(ggpubr)
library(stringr)
library(egg)
library(gt)
library(patchwork)
library(forestplot)
library(readxl)
library(ggsignif) 


## forestplot
dat = read.csv("04 hallmarks.meta.csv", header = F)
dat$V6 = log(dat$V6+1)
dat$V7 = log(dat$V7+1)
dat$V8 = log(dat$V8+1)
p1 = forestplot(#labeltext = as.matrix(dat[,c(1:2)]),
                labeltext = as.matrix(dat$V1),
                mean = dat$V6,
                lower = dat$V7,
                upper = dat$V8,
                zero = log(2),
                boxsize = 0.3,
                #lineheight = unit(8, "mm"),
                colgap = unit(3, "mm"),
                lwd.zero = 2,
                lwd.ci = 2,
                col = fpColors(box = "#FC4E07", lines = "#00AFBB", zero = "black"),
                xlab = "HR",
                xlog = TRUE,
                graph.pos = 2,
                ci.vertices = TRUE)
p1


dat = read.csv("04 hallmarks.meta.csv", header = T)
dat$group <- ifelse(dat$HR>1, "h", "l")
dat$group <- factor(dat$group, levels = c("l", "h"))
dat$HR <- log(dat$HR)

p1 <- ggdotchart(dat, x="id", y="HR",
                 color = "group",
                 palette = c("#00BFC4", "#F8766D"),
                 sorting = "descending",
                 add = "segments",
                 rotate = T,
                 dot.size = 2.5,
                 ggtheme = theme_article()) +
  geom_hline(yintercept=0, color="grey", size=0.5) +
  ggtitle("Meta analysis of cancer hallmarks") +
  theme(axis.text = element_text(colour = "black"),
        plot.title = element_text(hjust = 0.5, size = 10, face = "bold"),
        legend.position = 'none') +
  labs(x="", y="Log (HR)")
p1


## Grade boxplot
hallmarks = read.csv("Meta_CGGA.glycolysis.csv")
cli = read.csv("CGGA.cli.csv")

hallmarks = tibble::column_to_rownames(hallmarks, "ID")
hallmarks = as.data.frame(t(hallmarks))
hallmarks = tibble::rownames_to_column(hallmarks, "ID")
dat = dplyr::inner_join(cli, hallmarks, by = "ID")

table(dat$Grade)
dat1 = tidyr::drop_na(dat, Grade)

compaired <- list(c("WHO II","WHO III"), c("WHO II","WHO IV"), c("WHO IV","WHO III"))
p2 = ggplot(dat1, aes(x=Grade, y=Glycolysis)) +
  geom_boxplot(aes(fill=Grade)) +
  theme_article() +
  theme(legend.position = "none") +
  labs(y="Glycolysis ssGSEA score", x="") +
  geom_signif(comparisons = compaired,
              step_increase = 0.3,
              map_signif_level = T,
              test = t.test)
p2

## TCGA survival analysis
library(survival)
library(survminer)
library(egg)
library(ggthemes)

hallmarks = read.csv("Meta_TCGA.glycolysis.csv", check.names = F)
cli = read.csv("TCGA.cli.csv")

hallmarks = tibble::column_to_rownames(hallmarks, "ID")
hallmarks = as.data.frame(t(hallmarks))
hallmarks = tibble::rownames_to_column(hallmarks, "ID")
dat = dplyr::inner_join(cli, hallmarks, by = "ID")
write.csv(dat, "TCGA.ssGSEA.csv", row.names = F)

#table(dat$Cohort)
#dat1 = dplyr::filter(dat, dat$Cohort %in% c("TCGA_GBM", "TCGA_LGG"))
#cut_value <- surv_cutpoint(dat1, time="OS.time", event="OS", variables="Glycolysis")
#cut_value <- as.numeric(cut_value[["cutpoint"]]["cutpoint"])
#dat1$Group <- ifelse(dat1$Glycolysis>cut_value,'High','Low') 
dat$Group <- ifelse(dat$Glycolysis>median(dat$Glycolysis),'High','Low') 

fit <- survfit(Surv(OS.time, OS) ~ Group, data = dat)
p3 <- ggsurvplot(fit,
                 title = "TCGA cohort",
                 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 score", "High score"),
                 #palette = c("#00AFBB", "#FC4E07"),
                 legend = c(0.78,0.90),
                 legend.title = "Glycolysis",
                 xlim = c(0,150),
                 #break.x.by = 24,
                 surv.scale = "percent")
p3

##
table(dat$Cohort)
dat1 = dplyr::filter(dat, dat$Cohort %in% c("TCGA_LGG"))
cut_value <- surv_cutpoint(dat1, time="OS.time", event="OS", variables="Glycolysis")
cut_value <- as.numeric(cut_value[["cutpoint"]]["cutpoint"])
dat1$Group <- ifelse(dat1$Glycolysis>cut_value,'High','Low') 
dat1$Group <- ifelse(dat1$Glycolysis>median(dat1$Glycolysis),'High','Low') 

fit <- survfit(Surv(OS.time, OS) ~ Group, data = dat1)
p3 <- ggsurvplot(fit,
                 title = "TCGA_LGG",
                 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 score", "High score"),
                 palette = c("#00AFBB", "#FC4E07"),
                 legend = c(0.78,0.90),
                 legend.title = "Glycolysis",
                 xlim = c(0,150),
                 #break.x.by = 24,
                 surv.scale = "percent")
p3


##
table(dat$Cohort)
dat1 = dplyr::filter(dat, dat$Cohort %in% c("TCGA_GBM"))
cut_value <- surv_cutpoint(dat1, time="OS.time", event="OS", variables="Glycolysis")
cut_value <- as.numeric(cut_value[["cutpoint"]]["cutpoint"])
dat1$Group <- ifelse(dat1$Glycolysis>cut_value,'High','Low') 
dat1$Group <- ifelse(dat1$Glycolysis>median(dat1$Glycolysis),'High','Low') 

fit <- survfit(Surv(OS.time, OS) ~ Group, data = dat1)
p31 <- ggsurvplot(fit,
                 title = "TCGA_GBM",
                 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 score", "High score"),
                 palette = c("#00AFBB", "#FC4E07"),
                 legend = c(0.78,0.90),
                 legend.title = "Glycolysis",
                 xlim = c(0,150),
                 #break.x.by = 24,
                 surv.scale = "percent")
p31

## CGGA survival analysis

hallmarks = read.csv("Meta_CGGA.glycolysis.csv")
cli = read.csv("CGGA.cli.csv")

hallmarks = tibble::column_to_rownames(hallmarks, "ID")
hallmarks = as.data.frame(t(hallmarks))
hallmarks = tibble::rownames_to_column(hallmarks, "ID")
dat = dplyr::inner_join(cli, hallmarks, by = "ID")

#table(dat$Cohort)
#dat1 = dplyr::filter(dat, dat$Cohort %in% c("CGGA_301", "CGGA_325", "CGGA_693"))

#cut_value <- surv_cutpoint(dat1, time="OS.time", event="OS", variables="Glycolysis")
#cut_value <- as.numeric(cut_value[["cutpoint"]]["cutpoint"])
#dat1$Group <- ifelse(dat1$Glycolysis>cut_value,'High','Low') 

table(dat$cohort)
dat1 = dplyr::filter(dat, dat$cohort=="CGGA_325")
dat1$Group <- ifelse(dat1$Glycolysis>median(dat1$Glycolysis),'High','Low') 

fit <- survfit(Surv(OS.time, OS) ~ Group, data = dat1)
p4 <- ggsurvplot(fit,
                 title = "CGGA_325 cohort",
                 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 score", "High score"),
                 #palette = c("#00AFBB", "#FC4E07"),
                 legend = c(0.78,0.90),
                 legend.title = "Glycolysis",
                 xlim = c(0,150),
                 #break.x.by = 24,
                 surv.scale = "percent") 
p4


table(dat$cohort)
dat1 = dplyr::filter(dat, dat$cohort=="CGGA_693")
dat1$Group <- ifelse(dat1$Glycolysis>median(dat1$Glycolysis),'High','Low') 

fit <- survfit(Surv(OS.time, OS) ~ Group, data = dat1)
p5 <- ggsurvplot(fit,
                 title = "CGGA_693 cohort",
                 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 score", "High score"),
                 #palette = c("#00AFBB", "#FC4E07"),
                 legend = c(0.78,0.90),
                 legend.title = "Glycolysis",
                 xlim = c(0,150),
                 #break.x.by = 24,
                 surv.scale = "percent") 
p5



table(dat$Grade)

fit <- survfit(Surv(OS.time, OS) ~ Grade, data = dat)
p6 <- ggsurvplot(fit,
                 title = "CGGA cohort",
                 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 score", "High score"),
                 #palette = c("#00AFBB", "#FC4E07"),
                 legend = c(0.78,0.90),
                 legend.title = "Grade",
                 xlim = c(0,150),
                 #break.x.by = 24,
                 surv.scale = "percent") 
p6









## Grade survival analysis
table(dat$Cohort)
dat1 = tidyr::drop_na(dat, Grade)

cut_value <- surv_cutpoint(dat1, time="OS.time", event="OS", variables="Glycolysis")
cut_value <- as.numeric(cut_value[["cutpoint"]]["cutpoint"])
dat1$Group <- ifelse(dat1$Glycolysis>cut_value,'High','Low') 
dat1$Group <- ifelse(dat1$Glycolysis>median(dat1$Glycolysis),'High','Low') 

fit <- survfit(Surv(OS.time, OS) ~ Group, data = dat1)
p5 <- ggsurvplot(fit,
                 #title = "TCGA cohort",
                 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 score", "High score"),
                 palette = c("#00AFBB", "#FC4E07"),
                 legend = c(0.78,0.90),
                 legend.title = "Hypoxia",
                 xlim = c(0,150),
                 #break.x.by = 24,
                 surv.scale = "percent",
                 facet.by = "Grade") 
p5


## plot
library(patchwork)
pa = p1|(p2/p3$plot)

pb = p4$plot|p5$plot|p6$plot

pc = pa/pb




