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

library("limma")

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

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---#-------------------------

##--Median-----

filename.OS.riskMedian<- paste("The data for Figure 5 and 6-HNSC.OS.txt",sep = "")

scoreFile=filename.OS.riskMedian     #score

cliFile="The data for Figure 5 and 6-HNSC.clinical.txt"            

score=read.table(scoreFile, header=T, sep="\t",check.names=F, row.names=1)
cli=read.table(cliFile, header=T, sep="\t", check.names=F, row.names=1)

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

data=cbind(futime=score[,"futime"], fustat=score[,"fustat"], cli, group=score[,"risk"])
head(data)

#AJCC
cli.chageName=c("Age","Gender","Alcohol","Smoking","Grade",
                "T",
                "N",
                "M",
                "Stage",
                "Hpv")
#提取临床数据
#trait="Age.2"
for(trait in cli.chageName){
  
  rt=data[,c("futime","fustat",trait,"group")]
  #rt=rt[(rt[,trait]!="unknow"),]
  rt<-na.omit(rt)
  colnames(rt)=c("futime", "fustat", "clinical", "group")
  #对时间选择10年
  rt$fustat[rt$futime>10 ] <- 0
  rt$futime[rt$futime>10 ] <- 10
  
  tab=table(rt[,"clinical"])
  tab=tab[tab!=0]
  
  
  
  #对每个临床信息里面的分类进行循环
  for(j in names(tab)){
    rt1=rt[(rt[,"clinical"]==j),]
    tab1=table(rt1[,"group"])
    tab1=tab1[tab1!=0]
    labels=names(tab1)
    if(length(labels)==2){
      titleName=j
      if((trait=="age") | (trait=="Age") | (trait=="AGE")){
        titleName=paste0("age", j)
      }
      diff=survdiff(Surv(futime, fustat) ~group,data = rt1)
      #设置Pvalue的名字
      pValue.subOSName=paste0("pvalue.",trait,".",j,sep = "")
      pValue.OS=1-pchisq(diff$chisq, df=1)
      
      if(pValue.OS<0.001){
        pval=signif(pValue.OS,4)
        pval=format(pval, scientific = TRUE)
      }else{
        pVal=paste0("p=",sprintf("%.03f", pValue.OS))
      }
      fit <- survfit(Surv(futime, fustat) ~ group, data = rt1)
      
      bioCol=bioCol[1:length(levels(factor(rt1[,"group"])))]
      
      pValue.OS=round(pValue.OS,3)
      
      #绘制生存曲线
      surPlot=ggsurvplot(fit, 
                         data=rt1,
                         conf.int=F,
                         pval=pVal,
                         pval.size=6,
                         title=paste0("Patients with ",titleName),
                         legend.title=gene,
                         legend.labs=labels,
                         font.legend=12,
                         xlab="Time(years)",
                         break.time.by = 1,
                         palette=bioCol,
                         risk.table=TRUE,
                         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(geneName,"-OS.",trait,".",j,".pdf"), onefile = FALSE, width = 5.5, height =5)
      print(surPlot)
      dev.off()
    }
  }
  
}
 


#---DFS分析---#-------------------------

##--Median-------

#

scoreFile="The data for Figure 5 and 6-HNSC.DFS.txt"     #score
score=read.table(scoreFile, header=T, sep="\t",check.names=F, row.names=1)

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

data=cbind(futime=score[,"DFS_futime"], fustat=score[,"DFS_fustat"], cli, group=score[,"risk"])
head(data)
#提取临床数据
#trait="Age.2"
for(trait in cli.chageName){
  
  rt=data[,c("futime","fustat",trait,"group")]
  #rt=rt[(rt[,trait]!="unknow"),]
  rt<-na.omit(rt)
  colnames(rt)=c("futime", "fustat", "clinical", "group")
  #对时间选择10年
  rt$fustat[rt$futime>10 ] <- 0
  rt$futime[rt$futime>10 ] <- 10
  
  tab=table(rt[,"clinical"])
  tab=tab[tab!=0]
 
  #对每个临床信息里面的分类进行循环
  for(j in names(tab)){
    rt1=rt[(rt[,"clinical"]==j),]
    tab1=table(rt1[,"group"])
    tab1=tab1[tab1!=0]
    labels=names(tab1)
    if(length(labels)==2){
      titleName=j
      if((trait=="age") | (trait=="Age") | (trait=="AGE")){
        titleName=paste0("age", j)
      }
      diff=survdiff(Surv(futime, fustat) ~group,data = rt1)
      
      pvalue.subDFSName=paste0("pvalue.",trait,".",j,sep = "")
      
      pvalue.DFS=1-pchisq(diff$chisq, df=1)
      
      if(pvalue.DFS<0.001){
        pval=signif(pvalue.DFS,4)
        pval=format(pval, scientific = TRUE)
      }else{
        pval=paste0("p=",sprintf("%.03f", pvalue.DFS))
      }
      fit <- survfit(Surv(futime, fustat) ~ group, data = rt1)
      
      bioCol=bioCol[1:length(levels(factor(rt1[,"group"])))]
      
      #pvalue.DFS=format(pvalue.DFS, scientific = TRUE)
      pvalue.DFS=round(pvalue.DFS,3)
      
      #绘制生存曲线
      surPlot=ggsurvplot(fit, 
                         data=rt1,
                         conf.int=F,
                         pval=pval,
                         pval.size=6,
                         title=paste0("Patients with ",titleName),
                         legend.title=gene,
                         legend.labs=labels,
                         font.legend=12,
                         xlab="Time(years)",
                         ylab="DFS probability",
                         break.time.by = 1,
                         palette=bioCol,
                         risk.table=TRUE,
                         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(geneName,"-DFS.",trait,".",j,pvalue.DFS,".pdf"), onefile = FALSE, width = 5.5, height =5)
      print(surPlot)
      dev.off()
    }
  }
  
}


