
### 免疫治疗
## 读取数据
rm(list=ls())
load("IMvigor210_expcli.Rdata")

##
exp = i210exp_tpm
exp = dplyr::select_if(exp, ~!any(is.na(.)))

GeneCoef <- read.csv("0104 hypoxia genecoef.csv")
lassongene <- GeneCoef$Gene
gene <- intersect(exp$ID, lassongene)

meta = exp[lassongene,]
#meta = tibble::column_to_rownames(meta, "ID")
meta = as.data.frame(t(meta))

# 计算riskscore
actCoef = GeneCoef$Coef
FinalGeneExp = meta
myFun = function(x){crossprod(as.numeric(x),actCoef)} 
FinalGeneExp$riskscore = apply(FinalGeneExp,1,myFun)

##
surv = FinalGeneExp
surv$risk = ifelse(surv$riskscore > mean(surv$riskscore),'High','Low') 
surv$risk = factor(surv$risk, levels = c("Low", "High"))

surv = tibble::rownames_to_column(surv, "ID")
write.csv(surv, file = "0105 hypoxia.lasson.csv", row.names = F)

##
## ssGSEA数据前处理
gmtFile = "scRNA_marker.gmt" 
geneSet=getGmt(gmtFile, geneIdType=SymbolIdentifier())

dat = exp
dimnames=list(rownames(dat),colnames(dat))
mat=matrix(as.numeric(as.matrix(dat)),nrow=nrow(dat),dimnames=dimnames)
mat=avereps(mat)
mat=mat[rowMeans(mat)>0,]

## ssgsea
ssgseaScore=gsva(mat, geneSet, method='ssgsea', 
                 kcdf='Gaussian', abs.ranking=TRUE)

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

ssgseaOut = normalize(ssgseaScore)


## LASSO Cox
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)

GeneCoef <- read.csv("0201 TME genecoef.csv")
lassongene <- GeneCoef$Gene
gene <- intersect(row.names(ssgseaOut), lassongene)

# 计算riskscore
actCoef <- GeneCoef$Coef

ssgseaOut = as.data.frame(t(ssgseaOut))
FinalGeneExp <- ssgseaOut[,gene] 
myFun = function(x){crossprod(as.numeric(x),actCoef)} 
FinalGeneExp$riskscore = apply(FinalGeneExp,1,myFun)

##
surv = FinalGeneExp
surv$TME <- ifelse(surv$riskscore > mean(surv$riskscore),'Low','High') 
surv$TME <- factor(surv$TME, levels = c("Low", "High"))

surv = tibble::rownames_to_column(surv, "ID")
write.csv(surv, file = "0105 TME.lasson.csv", row.names = F)














lassongene <- GeneCoef$Gene
gene <- intersect(rownames(i210exp_tpm), lassongene)

meta <- as.data.frame(t(i210exp_tpm))
actCoef <- GeneCoef$Coef

# 计算riskscore
FinalGeneExp <- meta[,lassongene] 
myFun = function(x){crossprod(as.numeric(x),actCoef)} 
FinalGeneExp$riskscore = apply(FinalGeneExp,1,myFun)

# 合并表达数据和临床数据
FinalGeneExp$ID <- rownames(FinalGeneExp)
i210 <- inner_join(i210cli, FinalGeneExp, by="ID")


surv <- i210
cut_value <- surv_cutpoint(surv, time="OS_time", event="OS", variables="riskscore")
cut_value <- as.numeric(cut_value[["cutpoint"]]["cutpoint"])
surv$risk <- ifelse(surv$riskscore >cut_value,'High','Low') 

# 生存函数绘图
library(egg)
fit <- survfit(Surv(OS_time, OS) ~ risk, data = surv)
p6 <- 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")
p6


###
library(dplyr)
surv <- rename(surv, Response="BOR_binary")
dat = count(surv, risk, Response)
dat = dat %>% group_by(risk) %>% 
  summarise(Response = Response,n = n/sum(n))
dat$Response = factor(dat$Response,levels = c("SD/PD","CR/PR"))
dat$risk <- factor(dat$risk, levels = c("Low", "High"))

library(ggplot2)
p7 <- ggplot(data = dat)+
  geom_bar(aes(x = risk,y = n,
               fill = Response),
           stat = "identity")+
  scale_fill_manual(values =  c("#00BFC4","#F8766D"))+
  geom_text(aes(x = risk,y = n,
                label = scales::percent(n)),
            color = "white",
            size = 6,
            position = position_fill(vjust = 0.5))+
  theme_article()+
  theme(legend.position = "top") +
  labs(x=" ", y="Percent weight") + 
  scale_x_discrete(labels=c("Low"="Low GRS", "High"="High GRS"))
p7

###
fct_count(surv$Response)
compaired <- list(c("CR/PR","SD/PD"))
p8 <- ggplot(surv, aes(x=Response, y=riskscore)) +
  geom_boxplot(aes(fill=Response)) +
  theme_article() +
  theme(legend.position = "none",
        axis.text = element_text(colour = "black"),
        plot.title = element_text(hjust = 0.5)) +
  labs(x=" ", y="GRS") +
  stat_compare_means(comparisons = compaired,
                     label = "p.signif") +
  scale_fill_manual(values = c("#00BFC4", "#F8766D")) + 
  scale_x_discrete(labels=c("CR/PR"="Response", "SD/PD"="Nonresponse"))
p8












BiocManager::install(c('sva', 'car', 'genefilter', 'preprocessCore', 'ridge'))
rm(list = ls())
library(pRRophetic)
library(tidyverse)
library(survival)
library(survminer)
library(ggpubr)
library(ggsignif)
library(egg)

## 读取数据
lasson <- read.csv("03 TCGA.OV.HRSTME.csv", row.names = 1)

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

colnames(lasson)
exp <- lasson[, c(4:11)]
exp <- t(exp)
exp = as.matrix(exp)

## 预测化合物活性
# 顺铂 Cisplatin
# 多西紫杉醇 Docetaxel
# 阿霉素 Doxorubicin
# 吉西他滨 Gemcitabine
#  Mitomycin.C
# 紫杉醇 Paclitaxel



# Cisplatin
lasson$Cisplatin <- pRRopheticPredict(exp, "Cisplatin", selection = 1)

# Docetaxel
lasson$Docetaxel <- pRRopheticPredict(exp, "Docetaxel", selection = 1)

# Doxorubicin
lasson$Doxorubicin <- pRRopheticPredict(exp, "Doxorubicin", selection = 1)

# Gemcitabine
lasson$Gemcitabine <- pRRopheticPredict(exp, "Gemcitabine", selection = 1)

# Sorafenib
lasson$Sorafenib <- pRRopheticPredict(exp, "Sorafenib", selection = 1)


## plot
colnames(lasson)

dat <- lasson %>% pivot_longer(cols = 22:26,
                               names_to = "drug",
                               values_to = "IC50")

compaired <- list(c("Low","High"))
p1 <- ggplot(lasson, aes(x=risk, y=Cisplatin)) +
  geom_boxplot(aes(fill=risk)) +
  theme_article() +
  theme(legend.position = "none",
        axis.text = element_text(colour = "black"),
        plot.title = element_text(hjust = 0.5)) +
  labs(x=" ", y="Predicted Sensitivity") +
  ggtitle("Cisplatin") +
  stat_compare_means(comparisons = compaired,
                     label = "p.signif") +
  scale_fill_manual(values = c("#00BFC4", "#F8766D")) + 
  scale_x_discrete(labels=c("Low"="Low GRS", "High"="High GRS"))
p1

compaired <- list(c("Low","High"))
p2 <- ggplot(lasson, aes(x=risk, y=Docetaxel)) +
  geom_boxplot(aes(fill=risk)) +
  theme_article() +
  theme(legend.position = "none",
        axis.text = element_text(colour = "black"),
        plot.title = element_text(hjust = 0.5)) +
  labs(x=" ", y="Predicted Sensitivity") +
  ggtitle("Docetaxel") +
  stat_compare_means(comparisons = compaired,
                     label = "p.signif") +
  scale_fill_manual(values = c("#00BFC4", "#F8766D")) + 
  scale_x_discrete(labels=c("Low"="Low GRS", "High"="High GRS"))
p2


compaired <- list(c("Low","High"))
p3 <- ggplot(lasson, aes(x=risk, y=Doxorubicin)) +
  geom_boxplot(aes(fill=risk)) +
  theme_article() +
  theme(legend.position = "none",
        axis.text = element_text(colour = "black"),
        plot.title = element_text(hjust = 0.5)) +
  labs(x=" ", y="Predicted Sensitivity") +
  ggtitle("Doxorubicin") +
  stat_compare_means(comparisons = compaired,
                     label = "p.signif") +
  scale_fill_manual(values = c("#00BFC4", "#F8766D")) + 
  scale_x_discrete(labels=c("Low"="Low GRS", "High"="High GRS"))
p3


compaired <- list(c("Low","High"))
p4 <- ggplot(lasson, aes(x=risk, y=Gemcitabine)) +
  geom_boxplot(aes(fill=risk)) +
  theme_article() +
  theme(legend.position = "none",
        axis.text = element_text(colour = "black"),
        plot.title = element_text(hjust = 0.5)) +
  labs(x=" ", y="Predicted Sensitivity") +
  ggtitle("Gemcitabine") +
  stat_compare_means(comparisons = compaired,
                     label = "p.signif") +
  scale_fill_manual(values = c("#00BFC4", "#F8766D")) + 
  scale_x_discrete(labels=c("Low"="Low GRS", "High"="High GRS"))
p4


compaired <- list(c("Low","High"))
p5 <- ggplot(lasson, aes(x=risk, y=Sorafenib)) +
  geom_boxplot(aes(fill=risk)) +
  theme_article() +
  theme(legend.position = "none",
        axis.text = element_text(colour = "black"),
        plot.title = element_text(hjust = 0.5)) +
  labs(x=" ", y="Predicted Sensitivity") +
  ggtitle("Sorafenib") +
  stat_compare_means(comparisons = compaired,
                     label = "p.signif") +
  scale_fill_manual(values = c("#00BFC4", "#F8766D")) + 
  scale_x_discrete(labels=c("Low"="Low GRS", "High"="High GRS"))
p5

### 免疫治疗
## 读取数据
GeneCoef <- read.csv("GeneCoef.csv")
load("IMvigor210_expcli.Rdata")

lassongene <- GeneCoef$Gene
gene <- intersect(rownames(i210exp_tpm), lassongene)

meta <- as.data.frame(t(i210exp_tpm))
actCoef <- GeneCoef$Coef

# 计算riskscore
FinalGeneExp <- meta[,lassongene] 
myFun = function(x){crossprod(as.numeric(x),actCoef)} 
FinalGeneExp$riskscore = apply(FinalGeneExp,1,myFun)

# 合并表达数据和临床数据
FinalGeneExp$ID <- rownames(FinalGeneExp)
i210 <- inner_join(i210cli, FinalGeneExp, by="ID")


surv <- i210
cut_value <- surv_cutpoint(surv, time="OS_time", event="OS", variables="riskscore")
cut_value <- as.numeric(cut_value[["cutpoint"]]["cutpoint"])
surv$risk <- ifelse(surv$riskscore >cut_value,'High','Low') 

# 生存函数绘图
library(egg)
fit <- survfit(Surv(OS_time, OS) ~ risk, data = surv)
p6 <- 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")
p6


###
library(dplyr)
surv <- rename(surv, Response="BOR_binary")
dat = count(surv, risk, Response)
dat = dat %>% group_by(risk) %>% 
  summarise(Response = Response,n = n/sum(n))
dat$Response = factor(dat$Response,levels = c("SD/PD","CR/PR"))
dat$risk <- factor(dat$risk, levels = c("Low", "High"))

library(ggplot2)
p7 <- ggplot(data = dat)+
  geom_bar(aes(x = risk,y = n,
               fill = Response),
           stat = "identity")+
  scale_fill_manual(values =  c("#00BFC4","#F8766D"))+
  geom_text(aes(x = risk,y = n,
                label = scales::percent(n)),
            color = "white",
            size = 6,
            position = position_fill(vjust = 0.5))+
  theme_article()+
  theme(legend.position = "top") +
  labs(x=" ", y="Percent weight") + 
  scale_x_discrete(labels=c("Low"="Low GRS", "High"="High GRS"))
p7

###
fct_count(surv$Response)
compaired <- list(c("CR/PR","SD/PD"))
p8 <- ggplot(surv, aes(x=Response, y=riskscore)) +
  geom_boxplot(aes(fill=Response)) +
  theme_article() +
  theme(legend.position = "none",
        axis.text = element_text(colour = "black"),
        plot.title = element_text(hjust = 0.5)) +
  labs(x=" ", y="GRS") +
  stat_compare_means(comparisons = compaired,
                     label = "p.signif") +
  scale_fill_manual(values = c("#00BFC4", "#F8766D")) + 
  scale_x_discrete(labels=c("CR/PR"="Response", "SD/PD"="Nonresponse"))
p8


#
library(patchwork)
(p1+p2+p3+p4+p5+plot_layout(nrow = 1))/(p8+p6$plot+p7) + plot_annotation(tag_levels = 'A')


