rm(list = ls())
library(GSVA)
library(tidyverse)
library(GSEABase)
library(GSVAdata)
load("merge.RDATA")

data=outTab


geneSets <- getGmt("uniSigGeneExp.gmt")
score=GSVA::gsva(data,geneSets,method="ssgsea",verbose=T,ssgsea.norm=TRUE)
score=as.data.frame(t(score))
colnames(score)="score"


scoreOut=rbind(id=colnames(score), score)
write.table(scoreOut, file="score.txt", sep="\t", quote=F, col.names=F)


library(survival)
library(survminer)
scoreFile="score.txt"     
cliFile="clinicaldata.txt"          


score=read.table(scoreFile, header=T, sep="\t", check.names=F, row.names=1)
sampleType=gsub("(.*?)\\_.*", "\\1", row.names(score))
score=cbind(score, sampleType)
library(tidyverse)
rownames(score)=str_replace_all(rownames(score),"TCGA_","")
rownames(score)=str_replace_all(rownames(score),"GSE39582_","")
rownames(score)=str_replace_all(rownames(score),"GSE17538_","")

cli=read.table(cliFile, header=T, sep="\t", check.names=F, row.names=1)
cli=cli[,1:2]
cli$futime=cli$futime/365

sameSample=intersect(row.names(score), row.names(cli))
data=cbind(cli[sameSample,], score[sameSample,])

res.cut=surv_cutpoint(data, time="futime", event="fustat", variables=c("score"))
cutoff=as.numeric(res.cut$cutpoint[1])
print(cutoff)
Type=ifelse(data[,"score"]<=cutoff, "Low", "High")
data$group=Type
outTab=rbind(id=colnames(data), data)
write.table(outTab, file="score.group.txt", sep="\t", quote=F, col.names=F)

data$group=factor(data$group, levels=c("Low", "High"))
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("#0066FF","#FF0000","#6E568C","#7CC767","#223D6C","#D20A13","#FFD121","#088247","#11AA4D")
bioCol=bioCol[1:length]
surPlot=ggsurvplot(fit, 
                   data=data,
                   conf.int=F,
                   pval=pValue,
                   pval.size=6,
                   legend.title="score",
                   legend.labs=levels(factor(data[,"group"])),
                   legend = c(0.8, 0.8),
                   font.legend=12,
                   xlab="Time(years)",
                   break.time.by = 1,
                   palette = bioCol,
                   surv.median.line = "hv",
                   risk.table=T,
                   cumevents=F,
                   risk.table.height=.25)


pdf(file="survival.pdf", onefile = FALSE, width=7.5, height=5.5)
print(surPlot)
dev.off()