
#if (!requireNamespace("BiocManager", quietly = TRUE))
#    install.packages("BiocManager")
#BiocManager::install("limma", version = "3.8")

library("limma")

gene=geneName="FTL"                                                           #需要提取的基因

library("ggplot2")
library(survival)
library(survminer)

MyPal <- c(brewer.pal(12,"Set3")[c(3,6,8,10,11)],brewer.pal(8,"Set2")[-6])
MyPal <- c("#E41A1C","#377EB8","#4DAF4A","#FF7F00","#7CC767","#223D6C","#D20A13","#FFD121","#088247","#11AA4D")#（"#0571b0"-蓝,"#ca0020"-红）
bioCol=MyPal


#---OS---#-------------------------

##---
filename.OSDFS.exp=paste("The data for Figure S1C and D.txt",sep = "")
rt.OSDFS=read.table(filename.OSDFS.exp,header=T, row.names=1,sep="\t")                   #读取输入文件
rt=rt.OSDFS[,-(3:4)]

res.cut=surv_cutpoint(rt, time = "Follow_days", event = "Follow_status", variables = gene)
cutoff=as.numeric(res.cut$cutpoint[1])
print(cutoff)

Median=median(rt[,gene])
print(Median)

Type1=ifelse(rt[,gene]<=Median, "Low", "High")
Type2=ifelse(rt[,gene]<=cutoff, "Low", "High")
rt$risk=Type1
rt$riskCutoff=Type2

outTaball=rbind(ID=c("futime","fustat","riskscore","riskMedian","riskCutoff"), rt)
colnames(outTaball)=outTaball[1,]
outTaball=outTaball[-1,]
outTaball$futime[outTaball$futime <0 ] = 0
outTaball$futime=as.numeric(outTaball$futime)
outTaball$futime=outTaball$futime/365

#Median
outTabMedian=outTaball[,c("futime","fustat","riskscore","riskMedian")]
outTabMedian=rbind(ID=c("futime","fustat","riskscore","risk"), outTabMedian)

filename.OS.riskMedian<- paste(gene,"-OS.firebrowse.down-Median.txt",sep = "")
write.table(outTabMedian, file=filename.OS.riskMedian, sep="\t", quote=F, col.names=F)



##--Median------

#
data=read.table(filename.OS.riskMedian,header=T, row.names=1,sep="\t")                   #读取输入文件

#10years
data$fustat[data$futime>10 ] <- 0
data$futime[data$futime>10 ] <- 10

cli="OS"
surDatarisk = data

#---

diffrisk=survdiff(Surv(futime, fustat) ~ risk, data = surDatarisk)
lengthrisk=length(levels(factor(surDatarisk[,"risk"])))
pValuerisk=1-pchisq(diffrisk$chisq, df=1)

#if(pValuerisk<0.01){
#  pValuerisk="p<0.01"
#}else{
#  pValuerisk=paste0("p=",sprintf("%.03f",pValuerisk))
#}
pValuerisk=signif(pValuerisk,4)
pValuerisk=format(pValuerisk, scientific = TRUE)
pValuerisk.OS.Median=pValuerisk
fit.risk <- survfit(Surv(futime, fustat) ~ risk, data = surDatarisk)
#print(surv_median(fit))

#????????????
bioCol=MyPal
bioColrisk=bioCol[1:length(levels(factor(surDatarisk[,"risk"])))]
surPlotrisk=ggsurvplot(fit.risk, 
                       data=surDatarisk,
                       conf.int=F,
                       pval=pValuerisk,
                       pval.size=6,
                       legend.labs=levels(factor(surDatarisk[,"risk"])),
                       legend = c(0.8, 0.8),
                       legend.title=paste(gene),
                       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)

width=7
height=5.5
filename.tcga.OS.RiskMedian.pdf <- paste("OS.pdf",sep = "")
pdf(file=filename.tcga.OS.RiskMedian.pdf,onefile = FALSE, width=width, height=height)
print(surPlotrisk)
dev.off()



#---DFS---#-------------------------

filename.OSDFS.exp=paste("The data for Figure S1C and D.txt",sep = "")
rt.OSDFS=read.table(filename.OSDFS.exp,header=T, row.names=1,sep="\t")                   #读取输入文件
rt=rt.OSDFS[,-(1:2)]

res.cut=surv_cutpoint(rt, time = "DFS_days", event = "DFS_status", variables = gene)

cutoff=as.numeric(res.cut$cutpoint[1])
print(cutoff)

Median=median(rt[,gene])
print(Median)

Type1=ifelse(rt[,gene]<=Median, "Low", "High")
Type2=ifelse(rt[,gene]<=cutoff, "Low", "High")
rt$risk=Type1
rt$riskCutoff=Type2

outTaball=rbind(ID=c("futime","fustat","riskscore","riskMedian","riskCutoff"), rt)
colnames(outTaball)=outTaball[1,]
outTaball=outTaball[-1,]
outTaball$futime[outTaball$futime <0 ] = 0
outTaball$futime=as.numeric(outTaball$futime)
outTaball$futime=outTaball$futime/365

#--输出Median的分组
outTabMedian=outTaball[,c("futime","fustat","riskscore","riskMedian")]
outTabMedian=rbind(ID=c("futime","fustat","riskscore","risk"), outTabMedian)

filename.DFS.riskMedian<- paste(geneName,"-DFS.firebrowse.down-Median.txt",sep = "")
write.table(outTabMedian, file=filename.DFS.riskMedian, sep="\t", quote=F, col.names=F)



##--Median的分组

#比较高低风险组生存差异，得到显著性p值TCGA取10年，GEO取5年
data=read.table(filename.DFS.riskMedian,header=T, row.names=1,sep="\t")                   #读取输入文件

#对时间选择10年
data$fustat[data$futime>10 ] <- 0
data$futime[data$futime>10 ] <- 10

cli="DFS"
surDatarisk = data
#---生存分析
diffrisk=survdiff(Surv(futime, fustat) ~ risk, data = surDatarisk)
lengthrisk=length(levels(factor(surDatarisk[,"risk"])))
pValuerisk=1-pchisq(diffrisk$chisq, df=1)

#if(pValuerisk<0.01){
#  pValuerisk="p<0.01"
#}else{
#  pValuerisk=paste0("p=",sprintf("%.03f",pValuerisk))
#}
pValuerisk=signif(pValuerisk,4)
pValuerisk=format(pValuerisk, scientific = TRUE)
pValuerisk.DFS.Median=pValuerisk
fit.risk <- survfit(Surv(futime, fustat) ~ risk, data = surDatarisk)
#print(surv_median(fit))

#????????????
bioCol=MyPal
bioColrisk=bioCol[1:length(levels(factor(surDatarisk[,"risk"])))]
surPlotrisk=ggsurvplot(fit.risk, 
                       data=surDatarisk,
                       conf.int=F,
                       pval=pValuerisk,
                       pval.size=6,
                       legend.labs=levels(factor(surDatarisk[,"risk"])),
                       legend = c(0.8, 0.8),
                       legend.title=paste(geneName),
                       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)

width=7
height=5.5
filename.tcga.DFS.RiskMedian.pdf <- paste("DFS.pdf",sep = "")
pdf(file=filename.tcga.DFS.RiskMedian.pdf,onefile = FALSE, width=width, height=height)
print(surPlotrisk)
dev.off()



