load('GPR_TME_combined_group.Rdata')


library(survival)
library(survminer)
diff=survdiff(Surv(futime, fustat) ~ group, data = data)
length=length(levels(factor(data[,"group"])))
pValue=1-pchisq(diff$chisq, df=length-1)
if(pValue<0.001){
  pValue="p<0.001"
}else{
  pValue=paste0("p=",sprintf("%.03f",pValue))
}
data$group=factor(data$group,levels = c('GPR_high+TME_low','Mixed','GPR_low+TME_high'))
fit <- survfit(Surv(futime, fustat) ~ group, data = data)
bioCol <- c("#FFA500", "#8B4513", "#1E90FF")

p=ggsurvplot(fit, 
             data=data,
             conf.int=F,
             pval=pValue,
             pval.size=6,
             legend.title='GPR-TME score',
             legend.labs=levels(factor(data[,"group"])),
             legend = c(0.88, 0.9),
             font.legend=12,
             xlab="Time(Months)",
             palette = bioCol,
             surv.median.line = "hv",
             risk.table=T,
             cumevents=F,
             risk.table.height=.25)
p    




########################独立预后分析#################################
load('GPR_TME_combined_group.Rdata')
cli=read.table('clinical2.txt', header=T, sep="\t", check.names=F, row.names=1)      #读取临床文件
library(stringr)
data$id=stringr::str_sub(rownames(data),1,12)
data<- data[!duplicated(data$id), ]
rownames(data)<- data$id
data$id<-NULL
risk=data
## 注意此处，我们把三个亚型分成3层
risk$group=ifelse(risk$group=='GPR_low+TME_high',0,ifelse(risk$group=='Mixed',1,2))

risk=risk[,c('futime','fustat','group')]

#数据合并
sameSample=intersect(row.names(cli),row.names(risk))
risk=risk[sameSample,]
cli=cli[sameSample,]
rt=cbind(futime=risk[,1], fustat=risk[,2], cli, GPR_TME_classifier=risk[,3])

#单因素独立预后分析
uniTab=data.frame()
for(i in colnames(rt[,3:ncol(rt)])){
  cox <- coxph(Surv(futime, fustat) ~ rt[,i], data = rt)
  coxSummary = summary(cox)
  uniTab=rbind(uniTab,
               cbind(id=i,
                     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.table(uniTab,file='independent_uni.txt',sep="\t",row.names=F,quote=F)
#读取输入文件
rt <- read.table('independent_uni.txt', header=T, sep="\t", check.names=F, row.names=1)
gene <- rownames(rt)
hr <- sprintf("%.3f",rt$"HR")
hrLow  <- sprintf("%.3f",rt$"HR.95L")
hrHigh <- sprintf("%.3f",rt$"HR.95H")
Hazard.ratio <- paste0(hr,"(",hrLow,"-",hrHigh,")")
pVal <- ifelse(rt$pvalue<0.001, "<0.001", sprintf("%.3f", rt$pvalue))

#输出图形
n <- nrow(rt)
nRow <- n+1
ylim <- c(1,nRow)
layout(matrix(c(1,2),nc=2),width=c(3,2.5))

#绘制森林图左边的临床信息
xlim = c(0,3)
par(mar=c(4,2.5,2,1))
plot(1,xlim=xlim,ylim=ylim,type="n",axes=F,xlab="",ylab="")
text.cex=0.8
text(0,n:1,gene,adj=0,cex=text.cex)
text(1.5-0.5*0.2,n:1,pVal,adj=1,cex=text.cex);text(1.5-0.5*0.2,n+1,'pvalue',cex=text.cex,font=2,adj=1)
text(3.1,n:1,Hazard.ratio,adj=1,cex=text.cex);text(3.1,n+1,'Hazard ratio',cex=text.cex,font=2,adj=1)

#绘制右边的森林图
par(mar=c(4,1,2,1),mgp=c(2,0.5,0))
xlim = c(0,max(as.numeric(hrLow),as.numeric(hrHigh)))
plot(1,xlim=xlim,ylim=ylim,type="n",axes=F,ylab="",xaxs="i",xlab="Hazard ratio")
arrows(as.numeric(hrLow),n:1,as.numeric(hrHigh),n:1,angle=90,code=3,length=0.05,col="darkblue",lwd=2.5)
abline(v=1,col="black",lty=2,lwd=2)
boxcolor = ifelse(as.numeric(hr) > 1, '#0072b5', '#0072b5')
points(as.numeric(hr), n:1, pch = 15, col = boxcolor, cex=1.5)
axis(1)


dev.off()


#多因素独立预后分析

load('GPR_TME_combined_group.Rdata')
cli=read.table('clinical2.txt', header=T, sep="\t", check.names=F, row.names=1)      #读取临床文件
library(stringr)
data$id=stringr::str_sub(rownames(data),1,12)
data<- data[!duplicated(data$id), ]
rownames(data)<- data$id
data$id<-NULL


risk=data
risk$group=ifelse(risk$group=='GPR_low+TME_high',0,ifelse(risk$group=='Mixed',1,2))

risk=risk[,c('futime','fustat','group')]

#数据合并
sameSample=intersect(row.names(cli),row.names(risk))
risk=risk[sameSample,]
cli=cli[sameSample,]
rt=cbind(futime=risk[,1], fustat=risk[,2], cli, GPR_TME_classifier=risk[,3])

multiCox=coxph(Surv(futime, fustat) ~ ., data = rt)
multiCoxSum=summary(multiCox)
multiTab=data.frame()
ggforest(multiCox)

##############临床表型亚组生存分析##########
#install.packages("survival")
#install.packages("survminer")

load('GPR_TME_combined_group.Rdata')
cli=read.table('clinical.txt', header=T, sep="\t", check.names=F, row.names=1)      #读取临床文件
library(stringr)
data$id=stringr::str_sub(rownames(data),1,12)
data<- data[!duplicated(data$id), ]
rownames(data)<- data$id
data$id<-NULL

risk=data

#数据合并
sameSample=intersect(row.names(cli), row.names(risk))
risk=risk[sameSample,]
cli=cli[sameSample,]
data=cbind(futime=risk[,1],fustat=risk[,2],cli,risk=risk[,'group'])
data$risk
data$risk=factor(data$risk,levels = c('GPR_high+TME_low',"Mixed","GPR_low+TME_high"))
#对临床信息进行循环
for(i in colnames(data[,3:(ncol(data)-1)])){
  rt=data[,c("futime","fustat",i,"risk")]
  rt=rt[(rt[,i]!="unknow"),]
  colnames(rt)=c("futime","fustat","clinical","risk")
  tab=table(rt[,"clinical"])
  tab=tab[tab!=0]
  #对每个临床信息里面的分类进行循环
  for(j in names(tab)){
    rt1=rt[(rt[,"clinical"]==j),]
    tab1=table(rt1[,"risk"])
    tab1=tab1[tab1!=0]
    labels=names(tab1)
    ## 注意三层就=3
    if(length(labels)==3){
      titleName=j
      if((i=="age") | (i=="Age") | (i=="AGE")){
        titleName=paste0("age",j)
      }
      diff=survdiff(Surv(futime, fustat) ~risk,data = rt1)
      pValue=1-pchisq(diff$chisq,df=1)
      if(pValue<0.001){
        pValue="p<0.001"
      }else{
        pValue=paste0("p=",sprintf("%.03f",pValue))
      }
      fit <- survfit(Surv(futime, fustat) ~ risk, data = rt1)
      
      #绘制生存曲线
      surPlot=ggsurvplot(fit, 
                         data=rt1,
                         conf.int=F,
                         pval=pValue,
                         pval.size=6,
                         title=paste0("Patients with ",titleName),
                         legend.title="Risk",
                         legend.labs=labels,
                         font.legend=12,
                         xlab="Time(months)",
                         palette =c("#EFC000","#0073C2","#6E568C"),
                         risk.table=F,
                         risk.table.title="",
                         risk.table.col = "strata",
                         risk.table.height=.25)
      #输出图片
      j=gsub(">=","ge",j);j=gsub("<=","le",j);j=gsub(">","gt",j);j=gsub("<","lt",j)
      pdf(file=paste0("survival.",i,"_",j,".pdf"),onefile = FALSE,
          width = 8,        #图片的宽度
          height =6)        #图片的高度
      print(surPlot)
      dev.off()
    }
  }
}


############# 验证集###################
# 载入数据，这是一个我整理好的胃癌芯片数据
# 如果其他癌，请使用day0的预习内容芯片处理整理好GEO验证集！！！！！！！！！
load('ACRG.Rdata')

# 构建GPRscore
# 选择基因
gene=read.table('uniCox_lasso_gpr.txt',header = T)
gene=gene$gene
data=exprSet_acrg[gene,]

# 读取系数
coef=read.csv('GPR_coef.csv',header = T,check.names = F)
coef=coef$Coef..boot_SD
gpr_score=c()
for (i in 1:ncol(data)) {
  score=sum(as.numeric(data[,i])*coef)
  gpr_score=c(gpr_score,score)
}

data=as.data.frame(t(data))
data$gpr_score=gpr_score

#读取生存数据
suv=read.table('time.txt',row.names = 1,header = T,check.names = F)
cli=dplyr::select(suv,'futime','fustat')

##数据合并并输出结果
sameSample=intersect(row.names(data),row.names(cli))
data=data[sameSample,]
cli=cli[sameSample,]

rt=cbind(cli,data)

library(survival)
library(survminer)

### 中位值划分
Type=ifelse(data[,'gpr_score']<= median(rt$gpr_score), "Low", "High")
data=rt
data$group=Type
data$group=factor(data$group, levels=c("Low", "High"))
data_gpr=data

# CIBERSORT-Results-GEO.txt请根据预习内容自己评估好！
immune=read.table('CIBERSORT-Results-GEO.txt',sep = '\t',header = T,check.names = F,row.names = 1)
immune=immune[,-c(23:25)]


immune=as.data.frame(t(immune))

data=as.data.frame(t(immune))

data=data[,c('T cells CD8','T cells CD4 memory activated', 'Dendritic cells activated')]
data=as.data.frame(t(data))

coef=read.csv('TME_coef.csv',header = T,check.names = F)
coef=coef$Coef..boot_SD

TME_score=c()
for (i in 1:ncol(data)) {
  score=sum(as.numeric(data[,i])*coef)
  TME_score=c(TME_score,score)
}

data=as.data.frame(t(data))
data$TME_score=TME_score

#读取生存数据
#读取生存数据
suv=read.table('time.txt',row.names = 1,header = T,check.names = F)
cli=dplyr::select(suv,'futime','fustat')

##数据合并并输出结果
sameSample=intersect(row.names(data),row.names(cli))
data=data[sameSample,]
cli=cli[sameSample,]


rt=cbind(cli,data)
rt$futime=rt$futime/30

library(survival)
library(survminer)

### 中位值划分，注意此处TME_score负值越大，说明TME浸润越多！！！
### 所以我们认为TME_score越大，其实浸润越low
Type=ifelse(rt[,'TME_score'] <= median(rt$TME_score), "High", "Low")
data=rt
data$group=Type
data$group=factor(data$group, levels=c("Low", "High"))
data_immune=data



data_gpr$GPR_group=ifelse(data_gpr$group=='High','GPR_high','GPR_low')
data_gpr$TME_group=ifelse(data_immune$group=='High','TME_high','TME_low')

data=data_gpr
data$group=paste0(data$GPR_group,'+',data$TME_group)
data$group=stringr::str_replace(data$group,pattern = 'GPR_high\\+TME_high',replacement = 'Mixed')
data$group=stringr::str_replace(data$group,pattern = 'GPR_low\\+TME_low',replacement ='Mixed')

data$group

diff=survdiff(Surv(futime, fustat) ~ group, data = data)
length=length(levels(factor(data[,"group"])))
pValue=1-pchisq(diff$chisq, df=length-1)
if(pValue<0.001){
  pValue="p<0.001"
}else{
  pValue=paste0("p=",sprintf("%.03f",pValue))
}
fit <- survfit(Surv(futime, fustat) ~ group, data = data)
bioCol=c("#0073C2","#EFC000","#6E568C","#7CC767","#223D6C","#D20A13","#FFD121","#088247","#11AA4D")
bioCol=bioCol[1:length]

p=ggsurvplot(fit, 
             data=data,
             conf.int=F,
             pval=pValue,
             pval.size=6,
             legend.title='GPR-TME score',
             legend.labs=levels(factor(data[,"group"])),
             legend = c(0.7, 0.9),
             font.legend=12,
             xlab="Time(Months)",
             palette = bioCol,
             surv.median.line = "hv",
             risk.table=T,
             cumevents=F,
             risk.table.height=.25)
p    
