
options(scipen=100, digits=4)
rm(list = ls())
library(stringr)
library(ggpubr)
library(tidyverse)
library(survival)
library(survminer)
library(stringi)
library(egg)

##
library(immunedeconv)
exp = read.table("Meta_TCGA.exp.txt", sep = "\t", header = T, check.names = F, row.names = 1)
score = immunedeconv::deconvolute_estimate(exp)

score = as.data.frame(t(score))
score = tibble::rownames_to_column(score, "ID")
gsea <- read.csv("TCGA.ssGSEA.csv")
dat <- inner_join(score, gsea, by="ID")

colnames(dat)
#cut_value <- surv_cutpoint(dat, time="OS_time", event="OS", variables="riskscore")
#cut_value <- as.numeric(cut_value[["cutpoint"]]["cutpoint"])
#dat$risk <- ifelse(dat$Glycolysis >cut_value,'High','Low') 
dat$risk <- ifelse(dat$Glycolysis > median(dat$Glycolysis),'High','Low') 
dat$risk = factor(dat$risk, levels = c("Low","High"))

# plot
colnames(data)
compaired <- list(c("Low","High"))
p1 <- ggplot(dat, aes(x=risk, y=StromalScore)) +
  geom_violin(aes(fill=risk)) +
  theme_article() +
  geom_boxplot(width = 0.2, position = position_dodge(width = 1.0)) +
  theme(legend.position = "none",
        axis.text = element_text(colour = "black")) +
  labs(x="Glycolysis", y="StromalScore") +
  stat_compare_means(comparisons = compaired,
                     label = "p.signif") +
  scale_fill_manual(values = c("#00BFC4", "#F8766D"))

p1

p2 <- ggplot(dat, aes(x=risk, y=ImmuneScore)) +
  geom_violin(aes(fill=risk)) +
  theme_article() +
  geom_boxplot(width = 0.2, position = position_dodge(width = 1.0)) +
  theme(legend.position = "none",
        axis.text = element_text(colour = "black")) +
  labs(x="Glycolysis", y="ImmuneScore") +
  stat_compare_means(comparisons = compaired,
                     label = "p.signif") +
  scale_fill_manual(values = c("#00BFC4", "#F8766D"))

p2

p3 <- ggplot(dat, aes(x=risk, y=TumorPurity)) +
  geom_violin(aes(fill=risk)) +
  theme_article() +
  geom_boxplot(width = 0.2, position = position_dodge(width = 1.0)) +
  theme(legend.position = "none",
        axis.text = element_text(colour = "black")) +
  labs(x="Glycolysis", y="TumorPurity") +
  stat_compare_means(comparisons = compaired,
                     label = "p.signif") +
  scale_fill_manual(values = c("#00BFC4", "#F8766D"))

p3


## TIDE
dat = read.table("TCGA.GBM.logTPM.txt", header = T, sep = '\t', check.names = F, row.names = 1)
TIDE = dat
TIDE <- sweep(TIDE,2, apply(TIDE,2,median,na.rm=T))
TIDE <- sweep(TIDE,1, apply(TIDE,1,median,na.rm=T))
write.table(TIDE,"TIDE_input.txt",sep = "\t",row.names = T,col.names = NA,quote = F)

##
ssGSEA <- read.csv("TCGA.ssGSEA.csv", check.names = F)
TIDE <- read.csv("TIDE.csv", check.names = F)
dat = dplyr::inner_join(ssGSEA, TIDE, by = "ID")

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

dat$Group <- ifelse(dat$Glycolysis >median(dat$Glycolysis),'High','Low') 


###
colnames(dat)
dat1 = dat[,c("MDSC", "CAF", "TAM M2", "Group")]
dat1 = dat1 %>% pivot_longer(cols = 1:3, names_to = "gene", values_to = "expression")
dat1$Group = factor(dat1$Group, levels = c("Low", "High"))

p4 = ggplot(dat1, aes(x=Group, y= expression)) +
  geom_violin(aes(fill=Group)) +
  theme_article() +
  geom_boxplot(width = 0.2, position = position_dodge(width = 1.0)) +
  theme(legend.position = "none", axis.text = element_text(colour = "black")) +
  scale_fill_manual(values = c("#00BFC4", "#F8766D"), labels=c("Low MRS", "High MRS")) +
  labs(x = "Glycolysis", y = "Immunosuppressive cells proportion") +
  facet_wrap(~gene, nrow = 1) +
  stat_compare_means(aes(x=Group, y= expression),
                     method="t.test",
                     symnum.args=list(cutpoints = c(0, 0.001, 0.01, 0.05, 1), 
                                      symbols = c("***", "**", "*", " ")),
                     label = "p.signif")

p4



## data merge
tcga.cell <- read.csv("infiltration_estimation_for_tcga.csv", check.names = F)
colnames(tcga.cell)
tcga.cell <- tcga.cell %>% dplyr::rename(ID = cell_type)
cell <- grep("CIBERSORT", colnames(tcga.cell), value = T)
cell <- cell[1:22]
immune.cell <- tcga.cell %>% dplyr::select("ID", cell)

lihc <- read.csv("TCGA.ssGSEA.csv")
rownames(lihc) <- lihc$ID
lihc$ID <- substr(lihc$ID, 1, 15)
immune.data <- inner_join(immune.cell, lihc, by="ID")

## risk group
dat <- immune.data
colnames(dat)
#cut_value <- surv_cutpoint(dat, time="OS_time", event="OS", variables="Glycolysis")
#cut_value <- as.numeric(cut_value[["cutpoint"]]["cutpoint"])
#dat$risk <- ifelse(dat$riskscore >cut_value,'High','Low')

dat$risk <- ifelse(dat$Glycolysis >median(dat$Glycolysis),'High','Low') 
dat$risk <- factor(dat$risk, levels = c("Low", "High"))

# data transform
colnames(dat)
data1 <- dat %>% pivot_longer(cols = contains("CIBERSORT"), 
                              names_to = "cell", values_to = "exp")
data <- data1 %>% dplyr::select(cell, exp, risk)
str(data)
#data$exp <- as.numeric(data$exp)

# plot 
data$cell <- data$cell %>% str_remove("_CIBERSORT")
median(data$exp)
p5 <- ggplot(data, aes(cell, exp)) + 
  geom_boxplot(aes(fill=risk), outlier.colour = NA) +
  labs(x="", y="Immune cell composition") +
  theme_article() +
  theme(axis.text.x = element_text(angle = 45,hjust = 1, vjust = 1)) +
  scale_fill_manual(values = c("#00BFC4", "#F8766D"), labels=c("Low Glycolysis", "High Glycolysis")) + 
  theme(legend.title=element_blank(), legend.position=c(0.80, 0.80)) + 
  scale_y_continuous(limits = c(0, 0.8))
p5

p5=p5+stat_compare_means(aes(group=risk),
                         method="t.test",
                         symnum.args=list(cutpoints = c(0, 0.001, 0.01, 0.05, 1), 
                                          symbols = c("***", "**", "*", " ")),
                         label = "p.signif") 

p5



###################
###################
library(patchwork)
library(ggplotify)

pa = p1|p2|p3
pb = pa|p4

p = p1/(p2|p3) + plot_annotation("A")
p











