
rm(list=ls())
library(tidyverse)
library(limma)
library(ggpubr)
library(egg)
library(survival)
library(survminer)

coulor2 = c("#00AFBB", "#FC4E07")
coulor3 = c("#00AFBB", "#E7B800", "#FC4E07")
coulor4 = c("#00AFBB", "#E7B800", "#7bc043", "#FC4E07")


## immune profile
exp = read.table("Meta_TCGA.exp.txt", header = T, sep = "\t", check.names = F)
exp <- tibble::column_to_rownames(exp, "ID")
exp <- avereps(exp)

checkpoint = read.csv("immunomodulator.csv")
genes = intersect(rownames(exp), checkpoint$gene)
dat = exp[genes,]

## clinical
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")

surv = dplyr::filter(data, Cohort %in% c("TCGA_GBM", "TCGA_LGG"))
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"))

col <- c("ID", "GRSscore", "TMEscore", "GRSTME")
cli <- surv[, col]

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

colnames(dat)
res_aov = data.frame()
for (i in colnames(dat)[2:97]) {
  aov = aov(dat[,i] ~ dat$GRSTME, data=dat)
  p = round(summary(aov)[[1]][["Pr(>F)"]][1], 6)
  res_aov = rbind(res_aov,
                  cbind(ID=i, p=p))
}

gene = read.csv("immunomodulator.csv")
res_aov = dplyr::full_join(gene, res_aov, by="ID")

aov = res_aov
table(aov$class)
mol = aov[aov$class=="Immunostimulator",]
mol = dplyr::arrange(mol, p) 
gene = mol$ID[1:5]

gene1 = c("CD27", "CD40", "CXCR4", "IL6", "CD28")
col1 = c("GRSTME", gene1)
dat2 = dat[,col1]

data = dat2 %>% pivot_longer(cols = 2:6, names_to = "gene", values_to = "expression")
p9 = ggboxplot(data,
               x="GRSTME",
               y= "expression",
               fill = "GRSTME",
               #palette = coulor3,
               #add = "jitter",
               outlier.shape = NA,
               shape = "GRSTME") +
  xlab("") + ylab("Relative expression (log2)") +
  theme(axis.title.x = element_blank(), axis.text.x=element_blank()) +
  facet_wrap(~gene, nrow = 1)

my_com = list(c("GRSHigh/TMELow", "Mixed"), 
              c("Mixed", "GRSLow/TMEHigh"), 
              c("GRSHigh/TMELow", "GRSLow/TMEHigh"))
p9 = p9 + stat_compare_means(comparisons = my_com, label = "p.signif")
p9



table(aov$class)
mol = aov[aov$class=="Immunoinhibitor",]
mol = dplyr::arrange(mol, p) 
gene = mol$ID[1:5]

gene1 = c("IDO1", "CTLA4", "LAG3", "PDCD1", "TGFB1")
col = c("GRSTME", gene1)
dat2 = dat[,col]

data = dat2 %>% pivot_longer(cols = 2:6, names_to = "gene", values_to = "expression")
p10 = ggboxplot(data,
               x="GRSTME",
               y= "expression",
               fill = "GRSTME",
               outlier.shape = NA,
               #palette = coulor3,
               #add = "jitter"
               ) +
  xlab("") + ylab("Relative expression (log2)") +
  theme(axis.title.x = element_blank(), axis.text.x=element_blank()) +
  facet_wrap(~gene, nrow = 1)

my_com = list(c("GRSHigh/TMELow", "Mixed"), 
              c("Mixed", "GRSLow/TMEHigh"), 
              c("GRSHigh/TMELow", "GRSLow/TMEHigh"))
p10 = p10 + stat_compare_means(comparisons = my_com, label = "p.signif")
p10


##
library(patchwork)
pa = p1|p2|p3|p4
pb = plot_spacer()|p7|(p8$plot)
pc = p9|p10
pd = pa/pb/pc + plot_annotation(tag_levels = 'A')
pd = p5|p6


  plot_spacer()|p3|p4$plot + plot_annotation(tag_levels = 'A') 
pb


