rm(list=ls())
options(stringsAsFactors = F)
library("dplyr")
getwd()
Rdata_dir='./Rdata'
Figure_dir='./Figure'
####Reading expression information
miRNA_GA=read.table("CGGA.mRNAseq_325.RSEM-genes.20190701.txt",header = T,sep = '\t',check.names = F)
dim(miRNA_GA)
miRNA_GA[1:4,1:4]
na.omit(miRNA_GA)[1:4,1:4]
dim(na.omit(miRNA_GA))

#####remove duplicate values of expression matrix 
expr=t(miRNA_GA)
write.table(expr,"expr.txt",quote = F,sep = "\t",col.names = F,row.names = T)
z=read.table(file = "expr.txt",sep = "\t",header = T,check.names = F)
expr=z
dim(expr)
d=expr[!duplicated(expr$Gene_Name),]
dim(d)
expr_matirx=d
colnames(expr_matirx)[1]="CGGA_ID"

####Reading clinical information


miRNA_HiSeq=read.table("CGGA.mRNAseq_325.clinical.20190701.txt" ,header = T,sep = '\t')
dim(miRNA_HiSeq)
miRNA_HiSeq[1:4,1:4]
na.omit(miRNA_HiSeq)[1:4,1:4]
dim(na.omit(miRNA_HiSeq))
####Extract the column of clinical information of interest and remove the duplicate barcode
colnames(miRNA_HiSeq)

meta=miRNA_HiSeq

dim(meta)
meta=(meta[complete.cases(meta[,2:10]),])
dim(meta)
table(meta$X_EVENT)
meta=na.omit(meta)
dim(meta)
table(duplicated(meta$CGGA_ID))
unique_sample=unique(meta$sample)           






meta=meta[!duplicated(meta$CGGA_ID),]     
dim(meta)                                       


meta=meta[grep(pattern="WHO II|WHO III",meta[,4]),]    ####Extract the forth column of the meta containing "WHO II or WHO III" characters
meta=meta[grep(pattern="1|0",meta[,9]),]    ####Extract the ninth column of the meta containing "1 or 0" characters
meta=meta[grep(pattern="1|0",meta[,10]),] ####Extract the tenth column of the meta containing "1 or 0" characters
dim(meta)
colnames(meta)
colnames(meta)=c("CGGA_ID","PRS_type","type","grade","gender","age","OS.time",
                 "vital_status","radition","chemotherapy","IDH_status","1p19q_status")
meta_expr=merge(meta,expr_matirx,by="CGGA_ID")  

save(expr_matirx,meta,meta_expr,
     file = file.path(Rdata_dir,'CGGA-LGG-RNAseq-example.Rdata')
) 

###############
rm(list=ls())
options(stringsAsFactors = F)  


Rdata_dir='./Rdata'
Figure_dir='./Figure'       
load( file = 
        file.path(Rdata_dir,'CGGA-LGG-RNAseq-example.Rdata'))

#######Name the column of the clinical data and sorting the data

meta$age_group=ifelse(meta$age>median(meta$age),'older','younger')

meta$grade=gsub("WHO II","II",meta$grade)
meta$grade=gsub("WHO III","III",meta$grade)
meta_expr=merge(meta,expr_matirx,by="CGGA_ID") 
save(meta,meta_expr,expr_matirx,
     file = 
       file.path(Rdata_dir,'CGGA-LGG-RNAseq-example2.Rdata')
)


#######
rm(list=ls())
options(stringsAsFactors = F)



Rdata_dir='./Rdata'
Figure_dir='./Figure'
load( file = 
        file.path(Rdata_dir,'CGGA-LGG-RNAseq-example2.Rdata')
)
write.table(meta,file = "meta.txt",quote = F,sep = "\t",row.names = F)

Ued_exprclinical=meta_expr
z=Ued_exprclinical[,-(2:13)]     ####Delete the 2th to 13th column of meta

z=t(z)
write.table(z,"z.txt",quote = F,sep = "\t",col.names = F,row.names = T)
z=read.table(file = "z.txt",sep = "\t",header = T,check.names = F)





#####multivariable analysis

rm(list=ls())
options(stringsAsFactors = F)
library(survival)
library(survminer)
Rdata_dir='./Rdata/'
Figure_dir='./figures/'
load( file = 
        file.path(Rdata_dir,'CGGA-LGG-RNAseq-example2.Rdata')
)

p=c(colnames(meta_expr))
grep("WEE1",p)
grep("CRTAC1",p)
grep("SEMA4G",p)
multicox_input=meta_expr[,c(1:13,23305,4745,19706)]
save(multicox_input,
     file = 
       file.path(Rdata_dir,'CGGA-159LGG-RNAseq-multicox_input.Rdata')
)
mRNA=multicox_input
write.table(mRNA,file = "mRNA.txt",quote = F,sep = "\t",row.names = F)
miRNA<-read.table(file="mRNA.txt",header=T,sep="\t",row.names = 1,check.names = F,stringsAsFactors = F) 

colnames(miRNA)
fmla1 <- as.formula(Surv(OS.time,vital_status)~WEE1+CRTAC1+SEMA4G)
mycox <- coxph(fmla1,data=miRNA)


summary(mycox)
####Select the selected genes for the rise score
fmla1 <- as.formula(Surv(OS.time,vital_status)~SEMA4G+CRTAC1
                    +WEE1)

mycox <- coxph(fmla1,data=miRNA)
risk_score<-predict(mycox,type="risk",newdata=miRNA)
risk_level<-as.factor(ifelse(risk_score>median(risk_score),"High","Low"))
write.table(cbind(id=rownames(cbind(miRNA[,1:12],risk_score,risk_level)),cbind(miRNA[,1:12],risk_score,risk_level)),"risk_score.txt",sep="\t",quote=F,row.names=F) 



####Positioning gene locations and  survival analysis
library(survival)
library(survminer)

inputmiRNA=multicox_input

risk_score=read.table(file = "risk_score.txt",sep = "\t",header = T,check.names = F)
inputmiRNA$OS.time=inputmiRNA$OS.time/365
risk_score$OS.time=risk_score$OS.time/365

group1<- ifelse(inputmiRNA$WEE1>median(inputmiRNA$WEE1),'high','low')
sfit1<- survfit(Surv(OS.time,vital_status) ~ group1, data = inputmiRNA)

group2<- ifelse(inputmiRNA$CRTAC1>median(inputmiRNA$CRTAC1),'high','low')
sfit2=survfit(Surv(OS.time,vital_status) ~ group2, data = inputmiRNA)

group3<- ifelse(inputmiRNA$SEMA4G>median(inputmiRNA$SEMA4G),'high','low')
sfit3=survfit(Surv(OS.time,vital_status) ~ group3, data = inputmiRNA)




sfit4=survfit(Surv(OS.time,vital_status) ~ risk_level, data = risk_score)

ggsurvplot(sfit4,data = risk_score,
           surv.median.line = "hv",
           legend.title="Three-genes Signature",xlab ="Time(years)",
           conf.int = TRUE,
           pval = TRUE,palette = c("orchid2","dodgerblue2"),
           ggtheme = theme_bw()+theme(plot.title = element_text(hjust = 0.5)) ,title = "CGGA LGGs",legend.labs=c("High risk","Low risk"))


splots <- list()
splots[[1]] <- ggsurvplot(sfit1,data = inputmiRNA,
                          surv.median.line = "hv",
                          legend.title="WEE1",xlab ="Time(years)",
                          conf.int = TRUE,
                          pval = TRUE,palette = c("orchid2","dodgerblue2"),
                          ggtheme = theme_bw()+theme(plot.title = element_text(hjust = 0.5)) ,title = "CGGA-LGGs",legend.labs=c("High","Low"))
splots[[2]] <- ggsurvplot(sfit2,data = inputmiRNA,
                          surv.median.line = "hv",
                          legend.title="CRTAC1",xlab ="Time(years)",
                          conf.int = TRUE,
                          pval = TRUE,palette = c("orchid2","dodgerblue2"),
                          ggtheme = theme_bw()+theme(plot.title = element_text(hjust = 0.5)) ,title = "CGGA-LGGs",legend.labs=c("High","Low"))
splots[[3]] <- ggsurvplot(sfit3,data = inputmiRNA,
                          surv.median.line = "hv",
                          legend.title="SEMA4G",xlab ="Time(years)",
                          conf.int = TRUE,
                          pval = TRUE,palette = c("orchid2","dodgerblue2"),
                          ggtheme = theme_bw()+theme(plot.title = element_text(hjust = 0.5)) ,title = "CGGA-LGGs",legend.labs=c("High","Low"))

splots[[4]]=ggsurvplot(sfit4,data = risk_score,
                       surv.median.line = "hv",
                       legend.title=" ",xlab ="Time(years)",
                       conf.int = TRUE,
                       pval = TRUE,palette = "jco",
                       ggtheme = theme_bw(),title = "TCGA LGG Overall Survival")

# Arrange multiple ggsurvplots and print the output
arrange_ggsurvplots(splots, print = TRUE,  ncol = 1, nrow = 3)
dev.off()

strata=list()
younger=risk_score[grep(pattern="younger",risk_score[,13]),] 

sfit5=survfit(Surv(OS.time,vital_status) ~ risk_level, data = younger)

strata[[1]]=ggsurvplot(sfit5,data = younger,
                       surv.median.line = "hv",
                       legend.title="Three-genes Signature",xlab ="Time(years)",
                       conf.int = TRUE,
                       pval = TRUE,palette = c("orchid2","dodgerblue2"),
                       ggtheme = theme_bw()+theme(plot.title = element_text(hjust = 0.5)) ,title = "CGGA LGGs-Younger",legend.labs=c("High risk","Low risk"))

old=risk_score[grep(pattern="old",risk_score[,13]),] 

sfit6=survfit(Surv(OS.time,vital_status) ~ risk_level, data = old)

strata[[2]]=ggsurvplot(sfit6,data = old,
                       surv.median.line = "hv",
                       legend.title="Three-genes Signature",xlab ="Time(years)",
                       conf.int = TRUE,
                       pval = TRUE,palette = c("orchid2","dodgerblue2"),
                       ggtheme = theme_bw()+theme(plot.title = element_text(hjust = 0.5)) ,title = "CGGA LGGs-Old",legend.labs=c("High risk","Low risk"))


male=risk_score[grep(pattern="Male",risk_score[,5]),] 

sfit7=survfit(Surv(OS.time,vital_status) ~ risk_level, data = male)

strata[[3]]=ggsurvplot(sfit7,data = male,
                       surv.median.line = "hv",
                       legend.title="Three-genes Signature",xlab ="Time(years)",
                       conf.int = TRUE,
                       pval = TRUE,palette = c("orchid2","dodgerblue2"),
                       ggtheme = theme_bw()+theme(plot.title = element_text(hjust = 0.5)) ,title = "CGGA LGGs-Male",legend.labs=c("High risk","Low risk"))


female=risk_score[grep(pattern="Female",risk_score[,5]),] 

sfit8=survfit(Surv(OS.time,vital_status) ~ risk_level, data = female)

strata[[4]]=ggsurvplot(sfit8,data = female,
                       surv.median.line = "hv",
                       legend.title="Three-genes Signature",xlab ="Time(years)",
                       conf.int = TRUE,
                       pval = TRUE,palette = c("orchid2","dodgerblue2"),
                       ggtheme = theme_bw()+theme(plot.title = element_text(hjust = 0.5)) ,title = "CGGA LGGs-Female",legend.labs=c("High risk","Low risk"))


G3=risk_score[grep(pattern="III",risk_score[,4]),] 

sfit9=survfit(Surv(OS.time,vital_status) ~ risk_level, data = G3)

strata[[5]]=ggsurvplot(sfit9,data = G3,
                       surv.median.line = "hv",
                       legend.title="Three-genes Signature",xlab ="Time(years)",
                       conf.int = TRUE,
                       pval = TRUE,palette = c("orchid2","dodgerblue2"),
                       ggtheme = theme_bw()+theme(plot.title = element_text(hjust = 0.5)) ,title = "CGGA LGGs-Grade III",legend.labs=c("High risk","Low risk"))

risk_score$grade=gsub("III","3",risk_score$grade)
G2=risk_score[grep(pattern="II",risk_score[,4]),]

sfit10=survfit(Surv(OS.time,vital_status) ~ risk_level, data = G2)

strata[[6]]=ggsurvplot(sfit10,data = G2,
                       surv.median.line = "hv",
                       legend.title="Three-genes Signature",xlab ="Time(years)",
                       conf.int = TRUE,
                       pval = TRUE,palette = c("orchid2","dodgerblue2"),
                       ggtheme = theme_bw()+theme(plot.title = element_text(hjust = 0.5)) ,title = "CGGA LGGs-Grade II",legend.labs=c("High risk","Low risk"))



IDH_Wild=risk_score[grep(pattern="Wildtype",risk_score[,11]),] 

sfit11=survfit(Surv(OS.time,vital_status) ~ risk_level, data = IDH_Wild)

strata[[7]]=ggsurvplot(sfit11,data = IDH_Wild,
                       surv.median.line = "hv",
                       legend.title="Three-genes Signature",xlab ="Time(years)",
                       conf.int = TRUE,
                       pval = TRUE,palette = c("orchid2","dodgerblue2"),
                       ggtheme = theme_bw()+theme(plot.title = element_text(hjust = 0.5)) ,title = "CGGA-LGGs IDH_wildtype group",legend.labs=c("High risk","Low risk"))

IDH_Mutant=risk_score[grep(pattern="Mutant",risk_score[,11]),] 

sfit12=survfit(Surv(OS.time,vital_status) ~ risk_level, data = IDH_Mutant)

strata[[8]]=ggsurvplot(sfit12,data = IDH_Wild,
                       surv.median.line = "hv",
                       legend.title="Three-genes Signature",xlab ="Time(years)",
                       conf.int = TRUE,
                       pval = TRUE,palette = c("orchid2","dodgerblue2"),
                       ggtheme = theme_bw()+theme(plot.title = element_text(hjust = 0.5)) ,title = "CGGA-LGGs IDH_Mutant group",legend.labs=c("High risk","Low risk"))


p19q_Non=risk_score[grep(pattern="Non-codel",risk_score[,12]),] 

sfit13=survfit(Surv(OS.time,vital_status) ~ risk_level, data = p19q_Non)

strata[[9]]=ggsurvplot(sfit13,data = p19q_Non,
                       surv.median.line = "hv",
                       legend.title="Three-genes Signature",xlab ="Time(years)",
                       conf.int = TRUE,
                       pval = TRUE,palette = c("orchid2","dodgerblue2"),
                       ggtheme = theme_bw()+theme(plot.title = element_text(hjust = 0.5)) ,title = "CGGA-LGGs 1p19qNon-codel group",legend.labs=c("High risk","Low risk"))

Codel=risk_score[grep(pattern="Codel",risk_score[,12]),] 

sfit14=survfit(Surv(OS.time,vital_status) ~ risk_level, data = Codel)

strata[[10]]=ggsurvplot(sfit13,data = p19q_Codel,
                        surv.median.line = "hv",
                        legend.title="Three-genes Signature",xlab ="Time(years)",
                        conf.int = TRUE,
                        pval = TRUE,palette = c("orchid2","dodgerblue2"),
                        ggtheme = theme_bw()+theme(plot.title = element_text(hjust = 0.5)) ,title = "CGGA-LGGs 1p19q-codel group",legend.labs=c("High risk","Low risk"))





arrange_ggsurvplots(strata, print = TRUE,  ncol = 3, nrow = 2)

#############################################################################################################
#######risk_distribution
rm(list=ls())
options(stringsAsFactors = F)
Rdata_dir='./Rdata/'
Figure_dir='./figures/'

load( file = 
        file.path(Rdata_dir,'CGGA-LGG-RNAseq-example2.Rdata')
)
load( file = 
        file.path(Rdata_dir,'CGGA-159LGG-RNAseq-multicox_input.Rdata')
)


new_dat=multicox_input

colnames(new_dat) 
s=Surv(OS.time,vital_status) ~ WEE1+CRTAC1+SEMA4G
model <- coxph(s, data = new_dat )
summary(model,data=new_dat)
options(scipen=1)
ggforest(model, data =new_dat, 
         main = "Hazard ratio", 
         cpositions = c(0.10, 0.22, 0.4), 
         fontsize = 1.0, 
         refLabel = "1", noDigits = 4)





fp <- predict(model,new_dat,type="risk");boxplot(fp)
fp <- predict(model,new_dat,type="expected") ;boxplot(fp)
plot(fp,phe$days)
fp <- predict(model,new_dat) ;boxplot(fp)
basehaz(model) 
library(Hmisc)
options(scipen=200)
with(new_dat,rcorr.cens(fp,Surv(OS.time,vital_status)  ))
# http://dni-institute.in/blogs/cox-regression-interpret-result-and-predict/

library(cowplot)
library(pheatmap)
# https://cran.r-project.org/web/packages/cowplot/vignettes/introduction.html
fp_dat=data.frame(s=1:length(fp),v=as.numeric(sort(fp )))
sur_dat=data.frame(s=1:length(fp),
                   t=meta[names(sort(fp )),'OS.time'] ,
                   e=meta[names(sort(fp )),'vital_status']  ) 
sur_dat=na.omit(sur_dat)
sur_dat$e=ifelse(sur_dat$e==0,'alive','death')
exp_dat=new_dat[names(sort(fp )),c(14,23,20)]             ####感兴趣基因所在列
plot.point=ggplot(fp_dat,aes(x=s,y=v))+geom_point();print(plot.point)
plot.sur=ggplot(sur_dat,aes(x=s,y=t))+geom_point(aes(col=e));print(plot.sur)

mycolors <- colorRampPalette(c("MediumBlue","white","red"))(256)
tmp=t(scale(exp_dat))
tmp[tmp > 1] = 1
tmp[tmp < -1] = -1
plot.h=pheatmap(tmp,col= mycolors,show_colnames = F,cluster_cols = F,cluster_rows = F,clustering_distance_rows = "correlation")
plot.h=pheatmap(tmp,col= mycolors,show_colnames = F,cluster_cols = F,clustering_distance_rows = "correlation")
plot.h$gtable



plot_grid(plot.point, plot.sur, plot.h$gtable,
          labels = c("A", "B","C"),
          align = 'v',ncol = 1)  


#####################
########boxplot
rm(list=ls())
options(stringsAsFactors = F)
library("tidyverse")
Rdata_dir='./Rdata/'
Figure_dir='./figures/'

load( file = 
        file.path(Rdata_dir,'TCGA-LGG-RNAseq-example.Rdata')
)
load( file = 
        file.path(Rdata_dir,'CGGA-159LGG-RNAseq-multicox_input.Rdata')
)

multicox_meta=merge(meta,multicox,by="sample")

dat=data.frame(WEE1=meta_expr$WEE1,CRTAC1=meta_expr$CRTAC1,SEMA4G=meta_expr$SEMA4G,
               Grade=meta_expr$grade)
save(dat,file = 'data_for_boxplot.Rdata')
head(dat)
drop_outliers=function(df_name,dep_col){
  a=df_name[,dep_col]
  iqr=IQR(a)
  q1=as.numeric(quantile(a,0.25))
  q3=as.numeric(quantile(a,0.75))
  bottem=q1-1.5*iqr
  top=q3+1.5*iqr
  df_return=df_name %>%
    filter(df_name[,dep_col]>=bottem & df_name[,dep_col]<=top)
  return(df_return)
}
mydat=dat
mydat=drop_outliers(dat,"WEE1") %>% drop_outliers("WEE1")%>% drop_outliers("WEE1")%>% drop_outliers("WEE1")%>% drop_outliers("WEE1")%>% drop_outliers("WEE1")%>% drop_outliers("WEE1")%>% drop_outliers("WEE1")%>% drop_outliers("WEE1")%>% drop_outliers("WEE1")%>% drop_outliers("WEE1")%>% drop_outliers("WEE1")%>% drop_outliers("WEE1")%>% drop_outliers("WEE1")%>% drop_outliers("WEE1")%>% drop_outliers("WEE1")%>% drop_outliers("WEE1")%>% drop_outliers("WEE1")%>% drop_outliers("WEE1")%>% drop_outliers("WEE1")%>% drop_outliers("WEE1")%>% drop_outliers("WEE1")%>% drop_outliers("WEE1")%>% drop_outliers("WEE1")%>% drop_outliers("WEE1")%>% drop_outliers("WEE1")%>% drop_outliers("WEE1")

mydat=drop_outliers(mydat,"CRTAC1") %>% drop_outliers("CRTAC1")%>% drop_outliers("CRTAC1")%>% drop_outliers("CRTAC1")%>% drop_outliers("CRTAC1")%>% drop_outliers("CRTAC1")%>% drop_outliers("CRTAC1")%>% drop_outliers("CRTAC1")%>% drop_outliers("CRTAC1")%>% drop_outliers("CRTAC1")%>% drop_outliers("CRTAC1")%>% drop_outliers("CRTAC1")%>% drop_outliers("CRTAC1")%>% drop_outliers("CRTAC1")%>% drop_outliers("CRTAC1")%>% drop_outliers("CRTAC1")%>% drop_outliers("CRTAC1")%>% drop_outliers("CRTAC1")%>% drop_outliers("CRTAC1")%>% drop_outliers("CRTAC1")%>% drop_outliers("CRTAC1")%>% drop_outliers("CRTAC1")%>% drop_outliers("CRTAC1")%>% drop_outliers("CRTAC1")%>% drop_outliers("CRTAC1")%>% drop_outliers("CRTAC1")


mydat=drop_outliers(mydat,"SEMA4G") %>% drop_outliers("SEMA4G")%>% drop_outliers("SEMA4G")%>% drop_outliers("SEMA4G")%>% drop_outliers("SEMA4G")%>% drop_outliers("SEMA4G")%>% drop_outliers("SEMA4G")%>% drop_outliers("SEMA4G")%>% drop_outliers("SEMA4G")%>% drop_outliers("SEMA4G")%>% drop_outliers("SEMA4G")%>% drop_outliers("SEMA4G")%>% drop_outliers("SEMA4G")%>% drop_outliers("SEMA4G")%>% drop_outliers("SEMA4G")%>% drop_outliers("SEMA4G")%>% drop_outliers("SEMA4G")%>% drop_outliers("SEMA4G")%>% drop_outliers("SEMA4G")%>% drop_outliers("SEMA4G")%>% drop_outliers("SEMA4G")



mydat=gather(mydat,key="Gene",value="Expression",WEE1:SEMA4G)
mydat=mydat[order(mydat[,1]),]


compare_means( Expression ~ Grade, data = mydat, 
               group.by = "Gene")
# Box plot facetted by "gene"
p=ggboxplot(mydat, x = "Grade", y = "Expression",
            color = "Grade", palette = "jco",
            add = "jitter",title = "CGGA LGGs",
            facet.by = "Gene", short.panel.labs = FALSE)+theme(plot.title=element_text(hjust=0.5))
# Use only p.format as label. Remove method name.
p + stat_compare_means()



ggboxplot(mydat, x = "gene", y = "Expression",
          color = "Grade", palette = "nejm",
          add = "jitter")
p + stat_compare_means(aes(Grade = Grade))


if(require('ggpubr')){
  library(ggpubr)
  # google search : ggpubr boxplot add p-value
  # http://www.sthda.com/english/articles/24-ggpubr-publication-ready-plots/76-add-p-values-and-significance-levels-to-ggplots/
  p <- ggboxplot(dat, x = "Grade", y = "WEE1",
                 color = "red", palette = "jco",
                 add = "jitter")
  #  Add p-value
  p + stat_compare_means()
}



###############################################################################################
#########ROC


rm(list=ls())
options(stringsAsFactors = F)

Rdata_dir='./Rdata/'
Figure_dir='./figures/'

load( file = 
        file.path(Rdata_dir,'CGGA-LGG-RNAseq-example2.Rdata')
)
load( file = 
        file.path(Rdata_dir,'CGGA-159LGG-RNAseq-multicox_input.Rdata')
)



risk_score=read.table(file = "risk_score.txt",sep = "\t",header = T,check.names = F) 
meta=risk_score


library(lars) 
library(glmnet) 
library(survival)
library(survminer)
library(timeROC)

########################### 1 year ROC comparison
library("survivalROC")
cutoff <- 365
ROC1= survivalROC(Stime=meta$OS.time, 
                  status=meta$vital_status,      
                  marker = meta$risk_score,     
                  predict.time =  cutoff, method="KM")

meta$risk_score2=0.597*meta_expr$EMP3+0.825*meta_expr$GSX2+0.942*meta_expr$EMILIN3

cutoff <- 365
ROC1b= survivalROC(Stime=meta$OS.time, 
                   status=meta$vital_status,      
                   marker = meta$risk_score2,     
                   predict.time =  cutoff, method="KM")

meta$risk_score3=0.012*meta_expr$IFI44+0.068*meta_expr$KIF2C+0.035*meta_expr$GNG12+0.02*meta_expr$EMP3

cutoff <- 365
ROC1c= survivalROC(Stime=meta$OS.time, 
                   status=meta$vital_status,      
                   marker = meta$risk_score3,     
                   predict.time =  cutoff, method="KM")

plot(ROC1$FP, ROC1$TP, type="l",col="red", xlim=c(0,1), ylim=c(0,1),   
     xlab=paste( "FP", "\n", "AUC = ",round(ROC1$AUC,3)), 
     ylab="TP",main="Time dependent ROC 1 year")
abline(0,1)
lines(ROC1b$FP, ROC1b$TP, type="l",col="blue",xlim=c(0,1), ylim=c(0,1))
lines(ROC1c$FP, ROC1c$TP, type="l",col="yellow",xlim=c(0,1), ylim=c(0,1))
legend(0.5,0.4,c(paste("AUC at our study =",round(ROC1$AUC,3)),paste("AUC at Chen X.P et al =",round(ROC1b$AUC,3)),paste("AUC at Chuang Zhang et al =",round(ROC1c$AUC,3))),x.intersp=1, y.intersp=0.8,
       lty= 1, lwd= 2, col=c("red","blue","yellow"),
       bty = "n",
       seg.len=1, cex=0.8) 
####################5 years ROC comparison

cutoff <- 1825
ROC5= survivalROC(Stime=meta$OS.time, 
                  status=meta$vital_status,      
                  marker = meta$risk_score,     
                  predict.time =  cutoff, method="KM")

meta$risk_score2=0.597*meta_expr$EMP3+0.825*meta_expr$GSX2+0.942*meta_expr$EMILIN3

cutoff <- 1825
ROC5b= survivalROC(Stime=meta$OS.time, 
                   status=meta$vital_status,      
                   marker = meta$risk_score2,     
                   predict.time =  cutoff, method="KM")

meta$risk_score3=0.012*meta_expr$IFI44+0.068*meta_expr$KIF2C+0.035*meta_expr$GNG12+0.02*meta_expr$EMP3

cutoff <- 1825
ROC5c= survivalROC(Stime=meta$OS.time, 
                   status=meta$vital_status,      
                   marker = meta$risk_score3,     
                   predict.time =  cutoff, method="KM")

plot(ROC5$FP, ROC5$TP, type="l",col="red", xlim=c(0,1), ylim=c(0,1),   
     xlab=paste( "FP", "\n", "AUC = ",round(ROC1$AUC,3)), 
     ylab="TP",main="Time dependent ROC 5 years")
abline(0,1)
lines(ROC5b$FP, ROC5b$TP, type="l",col="blue",xlim=c(0,1), ylim=c(0,1))
lines(ROC5c$FP, ROC5c$TP, type="l",col="yellow",xlim=c(0,1), ylim=c(0,1))
legend(0.5,0.4,c(paste("AUC at our study =",round(ROC5$AUC,3)),paste("AUC at Chen X.P et al =",round(ROC5b$AUC,3)),paste("AUC at Chuang Zhang et al =",round(ROC5c$AUC,3))),x.intersp=1, y.intersp=0.8,
       lty= 1, lwd= 2, col=c("red","blue","yellow"),
       bty = "n",
       seg.len=1, cex=0.8) 
#########################################################3 years ROC comparison
cutoff <- 1095
ROC3= survivalROC(Stime=meta$OS.time, 
                  status=meta$vital_status,      
                  marker = meta$risk_score,     
                  predict.time =  cutoff, method="KM")

meta$risk_score2=0.597*meta_expr$EMP3+0.825*meta_expr$GSX2+0.942*meta_expr$EMILIN3

cutoff <- 1095
ROC3b= survivalROC(Stime=meta$OS.time, 
                   status=meta$vital_status,      
                   marker = meta$risk_score2,     
                   predict.time =  cutoff, method="KM")

meta$risk_score3=0.012*meta_expr$IFI44+0.068*meta_expr$KIF2C+0.035*meta_expr$GNG12+0.02*meta_expr$EMP3

cutoff <- 1095
ROC3c= survivalROC(Stime=meta$OS.time, 
                   status=meta$vital_status,      
                   marker = meta$risk_score3,     
                   predict.time =  cutoff, method="KM")

plot(ROC3$FP, ROC3$TP, type="l",col="red", xlim=c(0,1), ylim=c(0,1),   
     xlab=paste( "FP", "\n", "AUC = ",round(ROC1$AUC,3)), 
     ylab="TP",main="Time dependent ROC 3 years")
abline(0,1)
lines(ROC3b$FP, ROC3b$TP, type="l",col="blue",xlim=c(0,1), ylim=c(0,1))
lines(ROC3c$FP, ROC3c$TP, type="l",col="yellow",xlim=c(0,1), ylim=c(0,1))
legend(0.5,0.4,c(paste("AUC at our study =",round(ROC3$AUC,3)),paste("AUC at Chen X.P et al =",round(ROC3b$AUC,3)),paste("AUC at Chuang Zhang et al =",round(ROC3c$AUC,3))),x.intersp=1, y.intersp=0.8,
       lty= 1, lwd= 2, col=c("red","blue","yellow"),
       bty = "n",
       seg.len=1, cex=0.8) 





#####data preparation for Univariable and multivariable analysis 

rm(list=ls())
options(stringsAsFactors = F)



Rdata_dir='./Rdata/'
Figure_dir='./figures/'

load( file = 
        file.path(Rdata_dir,'CGGA-LGG-RNAseq-example2.Rdata')
)
load( file = 
        file.path(Rdata_dir,'CGGA-159LGG-RNAseq-multicox_input.Rdata')
)



risk_score=read.table(file = "risk_score.txt",sep = "\t",header = T,check.names = F)
colnames(risk_score)
colnames(risk_score)=c("Gene_Name","PRS_type","type","grade","gender","age","OS.time","vital_status",         
                       "radition","chemotherapy","IDH_status","p19q_status",
                       "age_group","risk_score","risk_level")


select_gene=expr_matirx[,c("Gene_Name","WEE1","CRTAC1","SEMA4G")]



Coxanaylye=merge(risk_score,select_gene,by="Gene_Name")

###########Univariable and multivariable analysis
library(survival)

miRNA=Coxanaylye
coxR=data.frame()
coxf<-function(x){
  fmla1 <- as.formula(Surv(OS.time,vital_status)~miRNA[,x])
  mycox <- coxph(fmla1,data=miRNA)
}
for(a in colnames(miRNA[,c(4:5,9:13,15)])){
  mycox=coxf(a)
  coxResult = summary(mycox)
  coxR=rbind(coxR,cbind(miRNAname=a,HR=coxResult$coefficients[,"exp(coef)"],
                        P=coxResult$coefficients[,"Pr(>|z|)"]))
}
write.table(coxR,"coxResult.txt",sep="\t",row.names=F,quote=F)

grade=coxph(Surv(OS.time,vital_status)~grade,data = miRNA)
summary(grade)

gender=coxph(Surv(OS.time,vital_status)~gender,data = miRNA)
summary(gender)

miRNA$radition=gsub("1","YES",miRNA$radition) 
miRNA$radition=gsub("0","NO",miRNA$radition)

radition=coxph(Surv(OS.time,vital_status)~radition,data = miRNA)
summary(radition)

miRNA$chemotherapy=gsub("1","YES",miRNA$chemotherapy) 
miRNA$chemotherapy=gsub("0","NO",miRNA$chemotherapy)

chemotherapy=coxph(Surv(OS.time,vital_status)~chemotherapy,data = miRNA)
summary(chemotherapy)

table(miRNA$IDH_status)
IDH_status=coxph(Surv(OS.time,vital_status)~IDH_status,data = miRNA)
summary(IDH_status)
table(miRNA$p19q_status)
p19q_status=coxph(Surv(OS.time,vital_status)~p19q_status,data = miRNA)
summary(p19q_status)

age_group=coxph(Surv(OS.time,vital_status)~age_group,data = miRNA)
summary(age_group)

risk_level=coxph(Surv(OS.time,vital_status)~risk_level,data = miRNA)
summary(risk_level)

table(miRNA$type)

Anaplastic_astrocytomas=miRNA[grep(pattern="AA|rAA",miRNA[,3]),]    ####提取meta中第7列含有“YES或者NO”两种字符的行
table(Anaplastic_astrocytomas$radition)
Anaplastic_astrocytomas_cox=coxph(Surv(OS.time,vital_status)~radition,data = Anaplastic_astrocytomas)
summary(Anaplastic_astrocytomas_cox)

Anaplastic_oligoastrocytomas=miRNA[grep(pattern="AOA|rAOA",miRNA[,3]),]    ####提取meta中第7列含有“YES或者NO”两种字符的行
table(Anaplastic_oligoastrocytomas$radition)
Anaplastic_oligoastrocytomas_cox=coxph(Surv(OS.time,vital_status)~radition,data = Anaplastic_oligoastrocytomas)
summary(Anaplastic_oligoastrocytomas_cox)



miRNA$type=gsub("AOA","Anaplastic_oligoastrocytomas",miRNA$type)                 
miRNA$type=gsub("rAOA","RAnaplastic_oligoastrocytomas",miRNA$type)


Anaplastic_oligodendrogliomas=miRNA[grep("AO|rAO",miRNA[,3]),]    ####提取meta中第7列含有“YES或者NO”两种字符的行
table(Anaplastic_oligodendrogliomas$radition)
Anaplastic_oligodendrogliomas_cox=coxph(Surv(OS.time,vital_status)~radition,data = Anaplastic_oligodendrogliomas)
summary(Anaplastic_oligodendrogliomas_cox)

table(miRNA$type)

which(miRNA$type=="A")
which(miRNA$type=="rA")
Astrocytoma=miRNA[c(which(miRNA$type=="A"),which(miRNA$type=="rA")),]
Astrocytoma$radition=gsub("1","Yes",Astrocytoma$radition)                        

table(Astrocytoma$radition)
Astrocytoma_cox=coxph(Surv(OS.time,vital_status)~radition,data = Astrocytoma)
summary(Astrocytoma_cox)

which(miRNA$type=="O")
Oligodendroglioma=miRNA[c(which(miRNA$type=="O")),]
table(Oligodendroglioma$radition)
Oligodendroglioma=coxph(Surv(OS.time,vital_status)~radition,data = Oligodendroglioma)

#######multivariable analysis
colnames(miRNA)
miRNA$radition=gsub("1","Yes",miRNA$radition)  
miRNA$chemotherapy=gsub("1","Yes",miRNA$chemotherapy)

fmla1 <- as.formula(Surv(OS.time,vital_status)~p19q_status+grade+gender+radition+chemotherapy+
                      IDH_status+age_group+risk_level)
mycox <- coxph(fmla1,data=miRNA)


summary(mycox)

table(miRNA$risk_level)