setwd('F:/') getwd() library(glmnet) library(survival) library(survminer) install.packages('rms') library(rms) library(readr) #ICC analysis install.packages('irr') library(irr) fpath <- choose.files("feature1.csv") feature_1 <- read_csv(fpath) fpath <- choose.files("feature2.csv") feature_2 <- read_csv(fpath) len <- 3562 icc_val<-vector(length=3562) thr <- 0.75 selected <- feature_1[feature_1$casename %in% feature_2$casename,] for (i in 2:3562){ ratings <- cbind(selected[,i],feature_2[,i]) icc <- icc(ratings, model = "twoway", type = "agreement", unit = "single", r0 = 0, conf.level = 0.95) icc_val[i] <- icc$value } Index <- which(icc_val > thr) dim(icc_val)=c(1,3562) write.csv(icc_val,file='F:/icc_val.csv') a<-dt[,1:8] b<-dt[,9:2836] b<-as.matrix(b) c<-scale(b,center=T,scale=T) dt3<-as.data.frame(cbind(a,c)) dt8<-read.csv('F:/vali-dt2.csv',header=T,skipNul=T) a<-dt8[,1:8] b<-dt8[,9:2836] b<-as.matrix(b) c<-scale(b,center=T,scale=T) dt9<-as.data.frame(cbind(a,c)) # univariate cox set.seed(99) outTab<-data.frame() outTab_UniCox<-data.frame() y<-Surv(time=train$OS,event=train$Status==2) for (i in colnames(train[,9:2836])) { cox<-coxph(y~train[,i],data = train) coxSummary<-summary(cox) outTab<-rbind(outTab,cbind(variables=i,HR=coxSummary$coefficients[,"exp(coef)"], z=coxSummary$coefficients[,"z"], pvalue=coxSummary$coefficients[,"Pr(>|z|)"])) if(coxSummary$coefficients[,"Pr(>|z|)"]>0.05) print(i) else outTab_UniCox=rbind(outTab_UniCox,cbind(i=i,coxSummary$coefficients[,"Pr(>|z|)"])) } dim(outTab_UniCox) dt_new <- data.frame(train$OS,train$Status) for (o in outTab_UniCox$i) { dt_new=cbind(dt_new,cbind(subset(train,select = o))) } #LASSO set.seed(99) x<-as.matrix(dt_new[,c(3:815)]) y<-data.matrix(Surv(dt_new$train.OS,dt_new$train.Status)) fitcv<-cv.glmnet(x,y,family="cox",alpha=1,nfolds = 10) print(fitcv) tiff(file="lambda.tiff",width = 18,height = 12,units = "cm", compression = "lzw",bg="white",res = 600) plot(fit,xvar="lambda",label=T)##xvar为x轴命名 dev.off() plot(fitcv) fitcv$lambda.min fitcv$lambda.1se which(coef(fitcv,s="lambda.min")!=0) which(coef(fitcv,s="lambda.1se")!=0) z<-dt_new[,1:2] b<-x[,c(1 , 6 , 9, 73 , 76, 107, 116, 159 ,200 ,225, 265, 340, 425, 468, 482, 494, 513, 546, 549, 559, 561 ,571,613, 618 ,624, 668, 669 ,693, 766, 783, 812)] lasso_dt<-cbind(z,b) ftrain<-coxph(Surv(train.OS,train.Status)~.,data = lasso_dt) sum.surv<-summary(ftrain) sum.surv train_c_index<-sum.surv$concordance train_c_index<-round(train_c_index,3) train_c_index ftest<-coxph(Surv(OS,Status)~predict(ftrain,newdata=test),data=test) sum.surv<-summary(ftest) sum.surv test_c_index<-sum.surv$concordance test_c_index<-round(test_c_index,3) test_c_index summary(ftest) fvali<-coxph(Surv(OS,Status)~predict(ftrain,newdata=dt9),data=dt9) summary(fvali) sum.surv<-summary(fvali) sum.surv vali_c_index<-sum.surv$concordance vali_c_index<-round(test_c_index,3) vali_c_index # high-low risk stratification chart riskScore<-predict(ftrain,type = "risk",newdata =lasso_dt ) dt2<-cbind(riskScore,lasso_dt) p<-median(dt2$riskScore) m<-dt1[,1] group<-ifelse(m>p,"high","low") datacut2<-cbind(group,dt2[,2:3]) fit2<- survfit(Surv(OS, Status) ~ group, data = datacut2) theme_set(theme_light()) ggsurvplot( fit2, title = "Survival Curve", ggtheme = theme(plot.title = element_text(hjust = 0.5), legend.text=element_text(size=14)), #调整字体大小 linetype = 1, surv.median.line = "hv", xlab = "Time", ylab = "Survival probability", legend = c(0.8,0.85), font.legend=c(17), legend.title = "riskscore", legend.labs = c("high", "low") , break.time.by = 1000, xlim = c(0,5000), break.y.by = .25, axes.offset = T, conf.int = TRUE, conf.int.alpha = .3, pval = TRUE, pval.size = 6, pval.coord = c(100,.2), palette = c("#E7B800","#2E9FDF"), # 设置调色板 risk.table = T, risk.table.fontsize=6, risk.table.col = c(1), fontsize = 4, font.x=c(14,'bold','darkred'), font.y=c(14,'bold','darkgreen'), font.xtickslab=c(12,'plain','blue'), font.title=c(20,'bold','black'), font.ytickslab=c(12,'plain','blue'), tables.theme=theme(axis.text.x = element_text(size=16), axis.text.y = element_text(size=16), axis.title.y = element_text(size=14), axis.title.x = element_text(size=14), plot.title = element_text(size=15)) ) dev.off() riskScore1<-predict(fvali,type = "risk",newdata =dt9) dt10<-cbind(riskScore1,dt9) m<-dt10[,1] group<-ifelse(m>1.143,"high","low") datacut2<-cbind(radscore,dt8[,8:9]) fit2<- survfit(Surv(OS, Status) ~ group, data = datacut2) ggsurvplot( fit2, linetype = 1, surv.median.line = "hv", palette = "aaas" , xlab = "Time", ylab = "Survival probability", title = "", legend = c(0.8,0.9), legend.title = "radscore", legend.labs = c("high", "low") , break.time.by = 1000, xlim = c(0,2500), break.y.by = .25, axes.offset = T, conf.int = TRUE, conf.int.alpha = .3, pval = TRUE, pval.size = 3, pval.coord = c(100,.3), risk.table = T, risk.table.col = c(1), fontsize = 4, ) dev.off() # multivariate analysis fit<-coxph(Surv(train.OS,train.Status==2)~age+IDH.Status+MGMT.Status+Grade, data =dt11) summary(fit) # subgroup risk stratification map windowsFonts(Times_New_Roman=windowsFont("Times New Roman")) tiff(file="F:/gender-female.tiff",width = 10,height = 8,units = "cm", compression = "lzw",bg="white",res = 600) dt1<-read.csv('F:/female.csv',header=T,skipNul=T) p<-median(dt1$riskScore) m<-dt1[,1] group<-ifelse(m>p,"high","low") datacut2<-cbind(group,dt1[,2:3]) fit2<- survfit(Surv(OS, Status) ~ group, data = datacut2) theme_set(theme_light()) ggsurvplot(fit2, data =datacut2, title = "Gender:female (Cut-off=1.143)", ggtheme = theme(plot.title = element_text(hjust = 0.5)), legend.labs = c("High Risk","Low Risk"), legend = 'top', legend.title='',# 图例位置 xlab = "Time(days)", break.time.by = 1000, xlim = c(0,3000), pval = TRUE, pval.size = 3, pval.coord = c(100,.10), pval.method = TRUE, pval.method.coord = c(100,.2), # 统计学方法文本的位置 palette = c("#E7B800","#2E9FDF"), # 设置调色板 font.x=c(14,'bold','darkred'), font.y=c(14,'bold','darkgreen'), font.xtickslab=c(12,'plain','blue'), font.title=c(16,'bold','black'), font.ytickslab=c(12,'plain','blue') ) dev.off() #Nomogram,calibration diagram,decision curve drawing dd<-datadist(dt5) options(datadist='dd') dt5<-read.csv('F:/dt5.csv',header=T,skipNul=T) str(dt5) for(i in names(dt5)[c(6:9)]){dt5[,i]<-as.factor(dt5[,i])} coxmodel<-cph(Surv(OS,Status) ~Grade+age+IDH.Status+MGMT.Status+riskScore ,data = dt5,x=T,y=T,surv=T) med<-Quantile(coxmodel) str(med) survival<-Survival(coxmodel) nom <- nomogram(coxmodel, fun=list(function(x) survival(365, x), function(x) survival(730, x), function(x) survival(1095, x)), lp=F, funlabel=c("1-year survival", "2-year survival", "3-year survival"), maxscale = 100, fun.at=c(0.05,seq(0.1,0.9,by=0.1),0.95) ) plot(nom，col.grid=c("pink","gray")) dev.off() dt_cal<-read.csv("dt15.csv") f1<-cph(formula = Surv(OS,Status)~age+IDH+MGMT+Grade+riskScore,data=dt_cal,x=T,y=T,surv = T,na.action=na.delete,time.inc = 365) cal1<-calibrate(f1, cmethod="KM", method="boot",u=365,m=10,B=1000) tiff(file="cal_1y.tiff",width = 14,height = 12,units = "cm", compression = "lzw",bg="white",res = 600) plot(cal1, lwd = 1,#error bar的粗细 lty = 6,#error bar的类型，可以是0-6 errbar.col = c("#228B22"),#error bar的颜色 xlim = c(0,1),ylim= c(0,1), xlab = "Nomogram-prediced OS (%)",ylab = "Observed OS (%)", col = c("#228B22"), cex.lab=1.0, cex.axis=1.0, cex.main=1.0, cex.sub=0.5) #字的大小 lines(cal1[,c('mean.predicted',"KM")], type = 'b', #连线的类型，可以是"p","b","o" lwd = 2, #连线的粗细 pch = 16, #点的形状，可以是0-20 col = c("#228B22")) #连线的颜色 mtext("") box(lwd = 1) #边框粗细 abline(0,1,lty = 3, #对角线为虚线 lwd = 2, #对角线的粗细 col = c("#224444")#对角线的颜色 ) dev.off() dt_cal<-read.csv("dt16.csv") f2<-cph(formula = Surv(OS,Status)~age+IDH+MGMT+Grade+riskScore,data=dt_cal,x=T,y=T,surv = T,na.action=na.delete,time.inc = 730) cal2<-calibrate(f2, cmethod="KM", method="boot",u=730,m=10,B=1000) tiff(file="cal_2y.tiff",width = 14,height = 12,units = "cm", compression = "lzw",bg="white",res = 600) plot(cal2, lwd = 1,#error bar的粗细 lty = 6,#error bar的类型，可以是0-6 errbar.col = c("#DAA520"),#error bar的颜色 xlim = c(0,1),ylim= c(0,1), xlab = "Nomogram-prediced OS (%)",ylab = "Observed OS (%)", col = c("#DAA520"), cex.lab=1.0, cex.axis=1, cex.main=1.0, cex.sub=0.5) #字的大小 lines(cal2[,c('mean.predicted',"KM")], type = 'b', #连线的类型，可以是"p","b","o" lwd = 2, #连线的粗细 pch = 16, #点的形状，可以是0-20 col = c("#DAA520")) #连线的颜色 mtext("") box(lwd = 1) #边框粗细 abline(0,1,lty = 3, #对角线为虚线 lwd = 2, #对角线的粗细 col = c("#224444")#对角线的颜色 ) dev.off() dt_cal<-read.csv("dt17.csv") f3<-cph(formula = Surv(OS,Status)~age+IDH+MGMT+Grade+riskScore,data=dt_cal,x=T,y=T,surv = T,na.action=na.delete,time.inc = 1095) cal3<-calibrate(f3, cmethod="KM", method="boot",u=1095,m=10,B=1000) tiff(file="cal_3y.tiff",width = 14,height = 12,units = "cm", compression = "lzw",bg="white",res = 600) plot(cal3, lwd = 1,#error bar的粗细 lty = 6,#error bar的类型，可以是0-6 errbar.col = c("#B22222"),#error bar的颜色 xlim = c(0,1),ylim= c(0,1), xlab = "Nomogram-prediced OS (%)",ylab = "Observed OS (%)", col = c("#B22222"), cex.lab=1.0, cex.axis=1, cex.main=1.0, cex.sub=0.5) #字的大小 lines(cal3[,c('mean.predicted',"KM")], type = 'b', #连线的类型，可以是"p","b","o" lwd = 2, #连线的粗细 pch = 16, #点的形状，可以是0-20 col = c("#B22222")) #连线的颜色 mtext("") box(lwd = 1) #边框粗细 abline(0,1,lty = 3, #对角线为虚线 lwd = 2, #对角线的粗细 col = c("#224444")#对角线的颜色 ) dev.off() library(devtools) library(ggDCA) library(caret) dt<-read.csv('F:/train-clin.csv',header=T,skipNul=T) radiomics<-coxph(Surv(OS,Status)~riskscore,data=dt) clinic_gene<-coxph(Surv(OS,Status)~age+Grade+IDH+MGMT,data=dt) radiomics_clinic_gene<-coxph(Surv(OS,Status)~age+ Grade+IDH+MGMT+riskscore,data=dt) dca_cph<-ggDCA::dca(radiomics,clinic_gene,radiomics_clinic_gene,times=365) library(ggsci) library(ggplot2) windowsFonts(Times_New_Roman=windowsFont("Times New Roman")) tiff(file="DCA.tiff",width = 18,height = 14,units = "cm", compression = "lzw",bg="white",res = 600) ggplot(dca_cph,linetype=1,lwd=0.5)+ scale_color_jama(name="", labels=c("radiomics","clinic_gene","radiomics_clinic_gene","All","None"))+ theme_bw(base_size = 14)+theme(legend.position = c(0.83,0.88), legend.background = element_blank())+ theme(legend.text=element_text(family="Times_New_Roman",face = "plain", size = 12, color = "black"))+ theme(axis.ticks.x=element_blank())+theme(axis.ticks.y=element_blank())+ theme(axis.text = element_text(family="Times_New_Roman",face = "plain", size = 11, color = "black"), axis.title.x = element_text(family="Times_New_Roman", face = "bold", size = 14, color = "black"), axis.title.y = element_text(family="Times_New_Roman",face = "bold", size = 14, color = "black"))+ scale_x_continuous(limits=c(0,1),expand=expansion(mult=c(0,0),add=c(0,0.02)))+scale_y_continuous(limits=c(-0.025,0.4), expand=expansion(mult=c(0,0),add=c(0,0))) dev.off() windowsFonts(Times_New_Roman=windowsFont("Times New Roman")) tiff(file="DCA.tiff",width = 18,height = 14,units = "cm", compression = "lzw",bg="white",res = 600) ggplot(dca_cph,linetype= 1,lwd=1 ) + scale_color_jama(name="", labels=c("radiomics model","clinical-molecular model","radiomics-clinical-molecular model","All","None"))+ theme_bw(base_size=14)+theme(legend.position = c(0.68,0.82), legend.background = element_blank())+ theme(legend.text=element_text(family="Times_New_Roman",face = "plain", size = 18, color = "black"))+ theme(axis.ticks.x=element_blank())+theme(axis.ticks.y=element_blank())+ theme(axis.text = element_text(family="Times_New_Roman",face = "plain", size = 15, color = "black"), axis.title.x = element_text(family="Times_New_Roman", face = "bold", size = 18, color = "black"), axis.title.y = element_text(family="Times_New_Roman",face = "bold", size = 18, color = "black"))+ scale_x_continuous(limits=c(0,1),expand=expansion(mult=c(0,0),add=c(0,0.02)))+scale_y_continuous(limits=c(-0.025,0.4), expand=expansion(mult=c(0,0),add=c(0,0))) dev.off() #AUC dt<-read.csv('F:/train-clic.csv',header=T,skipNul=T) library(timeROC) library(survival) library(ggplot2) time_roc_res <- timeROC(T = dt$OS, delta = dt$Status, marker = dt$riskScore1, cause = 2,weighting="marginal", times = c(1 * 365, 2 * 365, 3 * 365), ROC = TRUE,iid = TRUE) summary(time_roc_res) #Set working directory setwd('F:/') getwd() #Load the required R installation package library(glmnet) library(survival) library(survminer) install.packages('rms') library(rms) library(readr) # ICC analysis install.packages('irr') library(irr) fpath <- choose.files("feature1.csv") feature_1 <- read_csv(fpath) fpath <- choose.files("feature2.csv") feature_2 <- read_csv(fpath) len <- 3562 icc_val<-vector(length=3562) thr <- 0.75 selected <- feature_1[feature_1$casename %in% feature_2$casename,] for (i in 2:3562){ ratings <- cbind(selected[,i],feature_2[,i]) icc <- icc(ratings, model = "twoway", type = "agreement", unit = "single", r0 = 0, conf.level = 0.95) icc_val[i] <- icc$value } Index <- which(icc_val > thr) dim(icc_val)=c(1,3562) write.csv(icc_val,file='F:/icc_val.csv') #Load clinical and feature data after ICC analysis and normalize the data dt<-read.csv('F:/data.csv',header=T,skipNul=T) dt8<-read.csv('F:/External-validation.csv',header=T,skipNul=T) a<-dt[,1:8] b<-dt[,9:2836] b<-as.matrix(b) c<-scale(b,center=T,scale=T) dt3<-as.data.frame(cbind(a,c)) e<-dt8[,1:8] # Clinical variables f<-dt8[,9:2836] # Feature variables f<-as.matrix(f) g<-scale(f,center=T,scale=T) dt9<-as.data.frame(cbind(e,g)) #univariate cox analysis for selecting features set.seed(99) outTab<-data.frame() outTab_UniCox<-data.frame() y<-Surv(time=train$OS,event=train$Status==2) for (i in colnames(train[,9:2836])) { cox<-coxph(y~train[,i],data = train) coxSummary<-summary(cox) outTab<-rbind(outTab,cbind(variables=i,HR=coxSummary$coefficients[,"exp(coef)"], z=coxSummary$coefficients[,"z"], pvalue=coxSummary$coefficients[,"Pr(>|z|)"])) if(coxSummary$coefficients[,"Pr(>|z|)"]>0.05) print(i) else outTab_UniCox=rbind(outTab_UniCox,cbind(i=i,coxSummary$coefficients[,"Pr(>|z|)"])) } dim(outTab_UniCox) dt_new <- data.frame(train$OS,train$Status) for (o in outTab_UniCox$i) { dt_new=cbind(dt_new,cbind(subset(train,select = o))) } #Further select the optimal feature subset using LASSO set.seed(99) x<-as.matrix(dt_new[,c(3:815)]) y<-data.matrix(Surv(dt_new$train.OS,dt_new$train.Status)) fitcv<-cv.glmnet(x,y,family="cox",alpha=1,nfolds = 10) print(fitcv) tiff(file="lambda.tiff",width = 18,height = 12,units = "cm", compression = "lzw",bg="white",res = 600) plot(fit,xvar="lambda",label=T)##xvar为x轴命名 dev.off() plot(fitcv) fitcv$lambda.min fitcv$lambda.1se which(coef(fitcv,s="lambda.min")!=0) which(coef(fitcv,s="lambda.1se")!=0) z<-dt_new[,1:2] b<-x[,c(1 , 6 , 9, 73 , 76, 107, 116, 159 ,200 ,225, 265, 340, 425, 468, 482, 494, 513, 546, 549, 559, 561 ,571,613, 618 ,624, 668, 669 ,693, 766, 783, 812)] lasso_dt<-cbind(z,b) #Calculate the C index for each cohort ftrain<-coxph(Surv(train.OS,train.Status)~.,data = lasso_dt) sum.surv<-summary(ftrain) sum.surv train_c_index<-sum.surv$concordance train_c_index<-round(train_c_index,3) train_c_index ftest<-coxph(Surv(OS,Status)~predict(ftrain,newdata=test),data=test) sum.surv<-summary(ftest) sum.surv test_c_index<-sum.surv$concordance test_c_index<-round(test_c_index,3) test_c_index summary(ftest) fvali<-coxph(Surv(OS,Status)~predict(ftrain,newdata=dt9),data=dt9) summary(fvali) sum.surv<-summary(fvali) sum.surv vali_c_index<-sum.surv$concordance vali_c_index<-round(test_c_index,3) vali_c_index # high-low risk stratification of three cohorts riskScore<-predict(ftrain,type = "risk",newdata =lasso_dt ) dt2<-cbind(riskScore,lasso_dt) p<-median(dt2$riskScore) m<-dt2[,1] group<-ifelse(m>p,"high","low") datacut2<-cbind(group,dt2[,2:3]) fit2<- survfit(Surv(OS, Status) ~ group, data = datacut2) ggsurvplot( fit2, linetype = 1, surv.median.line = "hv", palette = "aaas" , xlab = "Time", ylab = "Survival probability", title = "", legend = c(0.8,0.9), legend.title = "radscore", legend.labs = c("high", "low") , break.time.by = 1000, xlim = c(0,2500), break.y.by = .25, axes.offset = T, conf.int = TRUE, conf.int.alpha = .3, pval = TRUE, pval.size = 3, pval.coord = c(100,.3), risk.table = T, risk.table.col = c(1), fontsize = 4, ) dev.off() riskScore<-predict(ftest,type = "risk",newdata =test ) dt7<-cbind(riskScore,test) n<-dt7[,1] group<-ifelse(n>1.143,"high","low") datacut1<-cbind(radscore,dt7[,8:9]) fit1<- survfit(Surv(OS, Status) ~ group, data = datacut1) ggsurvplot( fit1, linetype = 1, surv.median.line = "hv", palette = "aaas" , xlab = "Time", ylab = "Survival probability", title = "", legend = c(0.8,0.9), legend.title = "radscore", legend.labs = c("high", "low") , break.time.by = 1000, xlim = c(0,2500), break.y.by = .25, axes.offset = T, conf.int = TRUE, conf.int.alpha = .3, pval = TRUE, pval.size = 3, pval.coord = c(100,.3), risk.table = T, risk.table.col = c(1), fontsize = 4, ) dev.off() riskScore1<-predict(fvali,type = "risk",newdata =dt9) dt10<-cbind(riskScore1,dt9) m<-dt10[,1] group<-ifelse(m>1.143,"high","low") datacut2<-cbind(radscore,dt8[,8:9]) fit2<- survfit(Surv(OS, Status) ~ group, data = datacut2) ggsurvplot( fit2, linetype = 1, surv.median.line = "hv", palette = "aaas" , xlab = "Time", ylab = "Survival probability", title = "", legend = c(0.8,0.9), legend.title = "radscore", legend.labs = c("high", "low") , break.time.by = 1000, xlim = c(0,2500), break.y.by = .25, axes.offset = T, conf.int = TRUE, conf.int.alpha = .3, pval = TRUE, pval.size = 3, pval.coord = c(100,.3), risk.table = T, risk.table.col = c(1), fontsize = 4, ) dev.off() # Multivariate analysis for selecting clinical variables and constructing clinical model fit<-coxph(Surv(train.OS,train.Status==2)~age+IDH.Status+MGMT.Status+Grade, data =dt11) summary(fit) #Calculate the AUC values of the radiomics model for predicting 1 year,2 year and 3 year dt<-read.csv('F:/data.csv',header=T,skipNul=T) library(timeROC) library(survival) library(ggplot2) time_roc_res <- timeROC(T = dt$OS, delta = dt$Status, marker = dt$riskScore1, cause = 2,weighting="marginal", times = c(1 * 365, 2 * 365, 3 * 365), ROC = TRUE,iid = TRUE) summary(time_roc_res) ## subgroup risk stratification map windowsFonts(Times_New_Roman=windowsFont("Times New Roman")) tiff(file="F:/gender-female.tiff",width = 10,height = 8,units = "cm", compression = "lzw",bg="white",res = 600) dt1<-read.csv('F:/female.csv',header=T,skipNul=T) p<-median(dt1$riskScore) m<-dt1[,1] group<-ifelse(m>p,"high","low") datacut2<-cbind(group,dt1[,2:3]) fit2<- survfit(Surv(OS, Status) ~ group, data = datacut2) theme_set(theme_light()) ggsurvplot(fit2, data =datacut2, title = "Gender:female (Cut-off=1.143)", ggtheme = theme(plot.title = element_text(hjust = 0.5)), legend.labs = c("High Risk","Low Risk"), legend = 'top', legend.title='',# 图例位置 xlab = "Time(days)", break.time.by = 1000, xlim = c(0,3000), pval = TRUE, pval.size = 3, pval.coord = c(100,.10), pval.method = TRUE, pval.method.coord = c(100,.2), # 统计学方法文本的位置 palette = c("#E7B800","#2E9FDF"), # 设置调色板 font.x=c(14,'bold','darkred'), font.y=c(14,'bold','darkgreen'), font.xtickslab=c(12,'plain','blue'), font.title=c(16,'bold','black'), font.ytickslab=c(12,'plain','blue') ) dev.off() #Nomogram,calibration diagram,decision curve drawing dd<-datadist(dt5) options(datadist='dd') dt5<-read.csv('F:/dt5.csv',header=T,skipNul=T) str(dt5) for(i in names(dt5)[c(6:9)]){dt5[,i]<-as.factor(dt5[,i])} coxmodel<-cph(Surv(OS,Status) ~Grade+age+IDH.Status+MGMT.Status+riskScore ,data = dt5,x=T,y=T,surv=T) med<-Quantile(coxmodel) str(med) survival<-Survival(coxmodel) nom <- nomogram(coxmodel, fun=list(function(x) survival(365, x), function(x) survival(730, x), function(x) survival(1095, x)), lp=F, funlabel=c("1-year survival", "2-year survival", "3-year survival"), maxscale = 100, fun.at=c(0.1,seq(0.1,0.9,by=0.2),0.95) ) plot(nom,col.grid=c("pink","gray")) dev.off() dt_cal<-read.csv("dt15.csv") f1<-cph(formula = Surv(OS,Status)~age+IDH+MGMT+Grade+riskScore,data=dt_cal,x=T,y=T,surv = T,na.action=na.delete,time.inc = 365) cal1<-calibrate(f1, cmethod="KM", method="boot",u=365,m=10,B=1000) tiff(file="cal_1y.tiff",width = 14,height = 12,units = "cm", compression = "lzw",bg="white",res = 600) plot(cal1, lwd = 1, lty = 6, errbar.col = c("#228B22"), xlim = c(0,1),ylim= c(0,1), xlab = "Nomogram-prediced OS (%)",ylab = "Observed OS (%)", col = c("#228B22"), cex.lab=1.0, cex.axis=1.0, cex.main=1.0, cex.sub=0.5) lines(cal1[,c('mean.predicted',"KM")], type = 'b', lwd = 2, pch = 16, col = c("#228B22")) mtext("") box(lwd = 1) abline(0,1,lty = 3, lwd = 2, col = c("#224444") ) dev.off() dt_cal<-read.csv("dt16.csv") f2<-cph(formula = Surv(OS,Status)~age+IDH+MGMT+Grade+riskScore,data=dt_cal,x=T,y=T,surv = T,na.action=na.delete,time.inc = 730) cal2<-calibrate(f2, cmethod="KM", method="boot",u=730,m=10,B=1000) tiff(file="cal_2y.tiff",width = 14,height = 12,units = "cm", compression = "lzw",bg="white",res = 600) plot(cal2, lwd = 1, lty = 6, errbar.col = c("#DAA520"), xlim = c(0,1),ylim= c(0,1), xlab = "Nomogram-prediced OS (%)",ylab = "Observed OS (%)", col = c("#DAA520"), cex.lab=1.0, cex.axis=1, cex.main=1.0, cex.sub=0.5) lines(cal2[,c('mean.predicted',"KM")], type = 'b', lwd = 2, pch = 16, col = c("#DAA520")) mtext("") box(lwd = 1) abline(0,1,lty = 3, lwd = 2, col = c("#224444") ) dev.off() dt_cal<-read.csv("dt17.csv") f3<-cph(formula = Surv(OS,Status)~age+IDH+MGMT+Grade+riskScore,data=dt_cal,x=T,y=T,surv = T,na.action=na.delete,time.inc = 1095) cal3<-calibrate(f3, cmethod="KM", method="boot",u=1095,m=10,B=1000) tiff(file="cal_3y.tiff",width = 14,height = 12,units = "cm", compression = "lzw",bg="white",res = 600) plot(cal3, lwd = 1, lty = 6, errbar.col = c("#B22222"), xlim = c(0,1),ylim= c(0,1), xlab = "Nomogram-prediced OS (%)",ylab = "Observed OS (%)", col = c("#B22222"), cex.lab=1.0, cex.axis=1, cex.main=1.0, cex.sub=0.5) lines(cal3[,c('mean.predicted',"KM")], type = 'b', lwd = 2, pch = 16, col = c("#B22222")) mtext("") box(lwd = 1) abline(0,1,lty = 3, lwd = 2, col = c("#224444") ) dev.off() library(devtools) library(ggDCA) library(caret) dt<-read.csv('F:/train-clin.csv',header=T,skipNul=T) radiomics<-coxph(Surv(OS,Status)~riskscore,data=dt) clinic_gene<-coxph(Surv(OS,Status)~age+Grade+IDH+MGMT,data=dt) radiomics_clinic_gene<-coxph(Surv(OS,Status)~age+ Grade+IDH+MGMT+riskscore,data=dt) dca_cph<-ggDCA::dca(radiomics,clinic_gene,radiomics_clinic_gene,times=365) library(ggsci) library(ggplot2) windowsFonts(Times_New_Roman=windowsFont("Times New Roman")) tiff(file="DCA.tiff",width = 18,height = 14,units = "cm", compression = "lzw",bg="white",res = 600) ggplot(dca_cph,linetype=1,lwd=0.5)+ scale_color_jama(name="", labels=c("radiomics","clinic_gene","radiomics_clinic_gene","All","None"))+ theme_bw(base_size = 14)+theme(legend.position = c(0.83,0.88), legend.background = element_blank())+ theme(legend.text=element_text(family="Times_New_Roman",face = "plain", size = 12, color = "black"))+ theme(axis.ticks.x=element_blank())+theme(axis.ticks.y=element_blank())+ theme(axis.text = element_text(family="Times_New_Roman",face = "plain", size = 11, color = "black"), axis.title.x = element_text(family="Times_New_Roman", face = "bold", size = 14, color = "black"), axis.title.y = element_text(family="Times_New_Roman",face = "bold", size = 14, color = "black"))+ scale_x_continuous(limits=c(0,1),expand=expansion(mult=c(0,0),add=c(0,0.02)))+scale_y_continuous(limits=c(-0.025,0.4), expand=expansion(mult=c(0,0),add=c(0,0))) dev.off() windowsFonts(Times_New_Roman=windowsFont("Times New Roman")) tiff(file="DCA.tiff",width = 18,height = 14,units = "cm", compression = "lzw",bg="white",res = 600) ggplot(dca_cph,linetype= 1,lwd=1 ) + scale_color_jama(name="", labels=c("radiomics model","clinical-molecular model","radiomics-clinical-molecular model","All","None"))+ theme_bw(base_size=14)+theme(legend.position = c(0.68,0.82), legend.background = element_blank())+ theme(legend.text=element_text(family="Times_New_Roman",face = "plain", size = 18, color = "black"))+ theme(axis.ticks.x=element_blank())+theme(axis.ticks.y=element_blank())+ theme(axis.text = element_text(family="Times_New_Roman",face = "plain", size = 15, color = "black"), axis.title.x = element_text(family="Times_New_Roman", face = "bold", size = 18, color = "black"), axis.title.y = element_text(family="Times_New_Roman",face = "bold", size = 18, color = "black"))+ scale_x_continuous(limits=c(0,1),expand=expansion(mult=c(0,0),add=c(0,0.02)))+scale_y_continuous(limits=c(-0.025,0.4), expand=expansion(mult=c(0,0),add=c(0,0))) dev.off() #heatmap library(ggplot2) install.packages('ggthemes') library(ggthemes) install.packages('paletteer') library(paletteer) library(psych) library(tidyverse) dt1<-read.csv('F:/clinical.csv',header=T) dt2<-read.csv('F:/feature.csv',header=T) cor.result<-corr.test(dt1,dt2,method="spearman") cor.result$p cor.result$r cor.result$p%>% as.data.frame()%>% rownames_to_column()%>% pivot_longer(!rowname)%>% mutate(p_value=case_when( value>0.05~"A", value>0.01&value<=0.05~"B", value>0.001&value<=0.01~"D", value<=0.001~"E" ))%>% mutate(abs_cor=abs(value))-> wangge tiff(file="F:/heatmap.tiff",width = 18,height = 14,units = "cm", compression = "lzw",bg="white",res = 600) windowsFonts(A=windowsFont("Times New Roman")) p1<-ggplot()+geom_tile(data=wangge, aes(x=rowname,y=name,fill=p_value,alpha=p_value))+ scale_fill_manual(name='P_value', values=c("#F8F8F8","#D8D8D8", "#B0B0B0","#787878"), label=c(">0.05", "0.01~0.05", "0.001~0.01", "<0.01"))+ scale_alpha_manual(values=c(0,1,1,1))+ guides(alpha=F)+ theme_bw()+ theme(text = element_text(family = "A",face='bold'), axis.title=element_blank(), axis.ticks = element_blank(), legend.key=element_rect(colour="black"), axis.text.x=element_text(angle=45,hjust=1,vjust=1),)+ coord_equal()+ geom_point(data=wangge,aes(x=rowname,y=name, size=abs_cor, color=value))+ scale_color_paletteer_c(name="Spearman's r", limits=c(-0.25,1), palette="ggthemes::Classic Red-Blue", direction=-1)+ labs(size="Spearman's r") print(p1) dev.off() #Subgroup stratification modification library(survival) library(ggplot2) library(glmnet) library(survminer) install.packages('rms') library(rms) windowsFonts(Times_New_Roman=windowsFont("Times New Roman")) tiff(file="F:/idh-wild.tiff",width = 10,height = 8,units = "cm", compression = "lzw",bg="white",res = 600) dt1<-read.csv('F:/IDH-wild.csv',header=T,skipNul=T) p<-median(dt1$riskScore) m<-dt1[,1] group<-ifelse(m>p,"high","low") datacut2<-cbind(group,dt1[,2:3]) fit2<- survfit(Surv(OS, Status) ~ group, data = datacut2) theme_set(theme_light()) ggsurvplot(fit2, data =datacut2, title = "IDH:wild (Cut-off=1.143)", ggtheme = theme(plot.title = element_text(hjust = 0.5), legend.text=element_text(size=14)), #调整字体大小 legend.labs = c("High Risk","Low Risk"), legend = 'top', legend.title='',# 图例位置 xlab = "Time(days)", break.time.by = 1000, xlim = c(0,3000), pval = TRUE, pval.size = 6, pval.coord = c(2000,0.9), #调整P值位置 pval.method = FALSE, palette = c("#E7B800","#2E9FDF"), # 设置调色板 font.x=c(14,'bold','darkred'), font.y=c(14,'bold','darkgreen'), font.xtickslab=c(12,'plain','blue'), font.title=c(16,'bold','black'), font.ytickslab=c(12,'plain','blue') ) dev.off()