
rm(list=ls())
## ssGSEA 
library(GSVA)
library(GSEABase)
library(limma)
library(tidyverse)
library(ggsci)

##
## uniCox
ssgsea = read.csv("Cibersort.csv")
colnames(ssgsea)[1] = "ID"

cli <- read.csv("02 all.cli.csv")
cli = cli[,c("ID", "OS", "OS.time", "Cohort")]
#cli$ID <- substr(cli$ID, 1, 15)
meta1 <- dplyr::inner_join(cli, ssgsea, by="ID")
write.csv(meta1, file = "03 Cibersort.csv", row.names = F)

## 单因素Cox分析
library(survminer)
library(survival)
colnames(meta1)

meta = meta1[meta1$OS.time > 1,]

var = fct_unique(meta$Cohort)
var = as.character(var)

result = data.frame()
for (g in var) {
  dat = meta[meta$Cohort %in% g,]
  dat = dat %>% dplyr::select_if(~!any(is.na(.)))
  for(i in colnames(dat[,5:ncol(dat)])){
    cox <- coxph(Surv(OS.time, OS) ~ dat[,i], data = dat)
    coxSummary = summary(cox)
    coxP=coxSummary$coefficients[,"Pr(>|z|)"]
    result=rbind(result,
                 cbind(id=i,
                       Cohort=g,
                       coef=coxSummary[["coefficients"]],
                       HR=coxSummary$conf.int[,"exp(coef)"],
                       HR.95L=coxSummary$conf.int[,"lower .95"],
                       HR.95H=coxSummary$conf.int[,"upper .95"],
                       pvalue=coxSummary$coefficients[,"Pr(>|z|)"])
    )
  }
}

write.csv(result, file="03 Cibersort.uniCox.csv", row.names=F)


## metagen analysis
library(meta)

dat = read.csv("03 Cibersort.uniCox.csv")
dat$lnHR = log(dat$HR)
dat$lnL = log(dat$HR.95L)
dat$lnH = log(dat$HR.95H)
dat$seHR = (dat$lnH-dat$lnL)/(2*1.96)

result = data.frame()
var = unique(dat$id)
for (i in var) {
  dat1 = dplyr::filter(dat, dat$id %in% i)
  pfs=metagen(lnHR, seHR, sm="HR", data=dat1,
              studlab=paste(dat1$Cohort, dat1$id, sep="."))
  result=rbind(result,
               cbind(id=i,
                     HR=round(2.718282^pfs$TE.fixed, 3),
                     HR.95L=round(2.718282^pfs$lower.fixed, 3),
                     HR.95H=round(2.718282^pfs$upper.fixed, 3),
                     pvalue=round(pfs$pval.fixed, 3)))
}

write.csv(result, "03 TME.meta.csv", row.names = F)


## LASSO Cox
rm(list = ls())
library(survminer)
library(survival)
library(stringr)
library(tidyverse)
library(dplyr)
library(plyr)
library(glmnet)
library(pheatmap)
library(ggthemes)
library(patchwork)
library(ggplotify)
library(cowplot)
library(egg)

### 数据前处理
## screening gene
genes = read.csv("03 TME.meta.csv")
genes = genes$id[genes$pvalue < 0.05]
genes = genes$id[genes$HR < 1]

exp = read.csv("03 Cibersort.csv", check.names = F)
exp = dplyr::select(exp, c(c("ID", "OS", "OS.time", "Cohort"), all_of(genes)))
exp = dplyr::select_if(exp, ~!any(is.na(.))) #删除含有NA的列

## training data
table(exp$GSE)

data = exp
#data = dplyr::filter(exp, Cohort %in% c("CGGA_325", "CGGA_693"))
#data = dplyr::filter(exp, Cohort %in% c("TCGA_GBM", "TCGA_LGG"))

## lasso cox分析
x=data.matrix(data[ ,c(5:ncol(data))])
y=data.matrix(Surv(data$OS.time, data$OS))

set.seed(2023)
fit <- glmnet(x, y, family = "cox", alpha = 1)
pdf("lambda.pdf")
plot(fit, xvar = "lambda", label = TRUE)
dev.off()

cvfit <- cv.glmnet(x, y, family="cox", nfolds = 10, alpha=1)
pdf("cvfit.pdf")
plot(cvfit)
dev.off()

## 输出模型相关系数
coef = coef(fit, s = cvfit$lambda.min) 
index = which(coef != 0) 
actCoef = coef[index] 
lassoGene = row.names(coef)[index] 
geneCoef = cbind(Gene=lassoGene,Coef=actCoef) 
geneCoef
write.csv(geneCoef, file = "0201 TME genecoef.csv", row.names = F)

## riskscore
FinalGeneExp = exp[,lassoGene] 
myFun = function(x){crossprod(as.numeric(x),actCoef)} 
FinalGeneExp$TMEscore = apply(FinalGeneExp,1,myFun) 

##
cli = exp[,1:4]
data = cbind(cli, FinalGeneExp)
write.csv(data, file = "00202 TME.lasson.csv", row.names = F)



######## 绘制图形 ########
## KM生存曲线
# 分组
surv = dplyr::filter(data, Cohort %in% c("TCGA_GBM", "TCGA_LGG"))
surv$TME <- ifelse(surv$TMEscore > median(surv$TMEscore),'Low','High') 
#surv$risk <- factor(surv$risk, levels = c("Low", "High"))

# 生存函数绘图
fit <- survfit(Surv(OS.time, OS) ~ TME, data = surv)
p4 <- ggsurvplot(fit,
                 title = "",
                 pval = TRUE, 
                 #conf.int = TRUE,
                 risk.table = TRUE, 
                 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 GRS", "High GRS"),
                 palette = c("#00BFC4", "#F8766D"),
                 legend = c(0.85,0.90),
                 legend.title = "",
                 #xlim = c(0,72),
                 #break.x.by = 24,
                 surv.scale = "percent")
p4


# 绘制时间ROC
set.seed(2022)
library(timeROC)
result <-with(surv, timeROC(T=OS.time,
                           delta=OS,
                           marker=TMEscore,
                           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)
p5 <- ggplot() + 
  geom_line(data = plot_dat, aes(x = fpr, y = tpr, color = time), size = 1) + 
  scale_color_manual(name = NULL,values = c("#00BFC4", "#F8766D", "#e6b707"),
                     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")
p5


# 验证数据集
surv = dplyr::filter(data, Cohort %in% c("CGGA_325", "CGGA_693"))
surv$TME <- ifelse(surv$TMEscore >median(surv$TMEscore),'Low','High') 
#surv$risk <- factor(surv$risk, levels = c("Low", "High"))

# 生存函数绘图
fit <- survfit(Surv(OS.time, OS) ~ TME, data = surv)
p4 <- ggsurvplot(fit,
                 title = "",
                 pval = TRUE, 
                 #conf.int = TRUE,
                 risk.table = TRUE, 
                 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 GRS", "High GRS"),
                 palette = c("#00BFC4", "#F8766D"),
                 legend = c(0.85,0.90),
                 legend.title = "",
                 #xlim = c(0,72),
                 #break.x.by = 24,
                 surv.scale = "percent")
p4

# 绘制时间ROC
set.seed(2022)
library(timeROC)
result <-with(surv, timeROC(T=OS.time,
                            delta=OS,
                            marker=TMEscore,
                            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)
p5 <- ggplot() + 
  geom_line(data = plot_dat, aes(x = fpr, y = tpr, color = time), size = 1) + 
  scale_color_manual(name = NULL,values = c("#00BFC4", "#F8766D", "#e6b707"),
                     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")
p5



