
rm(list = ls())
library(GSEABase)
library(GSVA)
library(survival)
library(survminer)
library(My.stepwise)
library(egg)
library(tidyverse)
library(edgeR)

## ssGSEA
# exp <- read.csv("TCGA.OV_exp.csv", check.names = F)

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

geneSets <- getGmt("h.all.v2023.1.Hs.symbols.gmt")

ssgsea <- gsva(
  expr = exp,
  gset.idx.list = geneSets,
  mx.diff = F, # 数据为正态分布则T，双峰则F
  kcdf = "Gaussian", # CPM, RPKM, TPM数据就用默认值"Gaussian"， read count数据则为"Poisson"，
  parallel.sz = 10, # 并行线程数目
  method = "ssgsea",
  abs.ranking = TRUE)

## normalize
normalize=function(x){
  return((x-min(x))/(max(x)-min(x)))}

ssgsea = normalize(ssgsea)
ssgsea = as.data.frame(t(ssgsea))
colnames(ssgsea) = str_remove(colnames(ssgsea), "HALLMARK_")

## 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]

## data merge
ssgsea = tibble::rownames_to_column(ssgsea, "ID")
dat = dplyr::inner_join(cli, ssgsea, by = "ID")

dat$group = dat$GRSTME

## Anov
options(scipen = 200)
colnames(dat)
res_aov <- data.frame()
for (i in colnames(dat)[5:54]) {
  aov <- aov(dat[, i] ~ dat$group, data = dat)
  p <- round(summary(aov)[[1]][["Pr(>F)"]][1], 4)
  res_aov <- rbind(
    res_aov,
    cbind(ID = i, p = p)
  )
}

res_aov$sig <- ifelse(res_aov$p < 0.05,
                      ifelse(res_aov$p < 0.01,
                             ifelse(res_aov$p < 0.001, "***", "**"), "*"
                      ), ""
)
res_aov$ID1 <- paste(res_aov$ID, res_aov$sig, " ")

res_aov <- res_aov[order(res_aov$p), ]
cli_name <- res_aov[, c("ID", "ID1")]
write.csv(res_aov, "res_aov.csv", row.names = F)

### data merge
ssgsea = tibble::column_to_rownames(ssgsea, "ID")
ssgsea = as.data.frame(t(ssgsea))
ssgsea = tibble::rownames_to_column(ssgsea, "ID")
ssgsea = inner_join(cli_name, ssgsea, by = "ID")

ssgsea = ssgsea[,-1]
ssgsea = tibble::column_to_rownames(ssgsea, "ID1")
ssgsea = as.data.frame(t(ssgsea))
ssgsea = tibble::rownames_to_column(ssgsea, "ID")

cli_exp <- inner_join(cli, ssgsea, by = "ID")
write.csv(cli_exp, "TCGA.OV_cli_ssgsea.csv", row.names = F)

## plot
dat = read.csv("TCGA.OV_cli_ssgsea.csv",check.names = F)
colnames(dat)
dat = dat[,1:36]

dat$Hypoxia = ifelse(dat$riskscore > mean(dat$riskscore),'HypoxiaHigh', 'HypoxiaLow')
dat$TME = ifelse(dat$TMEscore > mean(dat$TMEscore),'TMELow', 'TMEHigh') 

dat$group = paste(dat$Hypoxia, dat$TME, sep = "/")
table(dat$group)

dat$group = forcats::fct_collapse(dat$group, Mixed = c("HypoxiaHigh/TMEHigh", "HypoxiaLow/TMELow"))
table(dat$group)

dat$group = factor(dat$group, 
                   levels = c("HypoxiaLow/TMEHigh", "Mixed","HypoxiaHigh/TMELow"),
                   labels = c("HypoxiaLow/TMEHigh", "Mixed","HypoxiaHigh/TMELow"))



annotation_col <- data.frame(row.names = dat$ID,
                             Hypoxia = dat$Hypoxia,
                             TME = dat$TME,
                             classifier = dat$group)
colnames(dat)
# annotation_col = annotation_col[,-1]

colnames(dat)
exp <- dat[,7:36]
rownames(exp) <- dat$ID
exp_plot <- as.data.frame(t(exp))

mycolors <-
  colorRampPalette(c("#0392cf", "white", "#ee4035"), bias = 1.2)(100)
# ann_colors = list(MRSTME = c(A="#00BFC4", B="#F8766D", C="red"))

p8 <- pheatmap::pheatmap(
  exp_plot,
  col = mycolors,
  annotation_col = annotation_col,
  # annotation_colors =ann_colors,
  scale = "row",
  breaks = seq(-2, 2, length.out = 150),
  show_colnames = F,
  cluster_cols = T,
  annotation_legend = T,
  cluster_rows = F,
  show_rownames = T)
p8

