Step 4：Hub gene expression in the training set and validation set
setwd("-")
load("SD_exp_combined-human.Rdata")
exp2=exp2[,-42]
exp2=as.data.frame(t(exp2))
dat1=cbind(group_list,exp2)
hub_gene<-as.data.frame(dat1[,c('CDKN1A','HSPA5','NR4A1','PRKAB1')])
dat2=cbind(group_list,hub_gene) 
library(ggstatsplot)
ggbetweenstats(data=dat2,x=group_list,y=CDKN1A)
ggbetweenstats(data=dat2,x=group_list,y=HSPA5)
ggbetweenstats(data=dat2,x=group_list,y=NR4A1)
ggbetweenstats(data=dat2,x=group_list,y=PRKAB1)


library(pcutils)
library(ggstatsplot)
library(PMCMRplus) 
library(ggplot2) 
library(cowplot) 

set.seed(123)

data <- read.table("表达矩阵.txt",sep = "\t",check.names = F,stringsAsFactors = F,header = T)

{
  plist = list()
  for (i in 1:9) {
    plist[[i]] = group_box(tab = dat2[c("PRKAB1")], group = "group_list", metadata = dat2, 
                           mode = i,p_value1 = "t.test") +#t.test，wilcox.test，anova，kruskal.test等方法进行比较
      ggtitle(paste0("mode", i)) +
      theme_classic() + 
      theme(legend.position = "none")
  }
  
  plot_grid(plotlist = plist, ncol = 3)
}
plot(plist[[4]])

library(ggplot2)
library(ggpubr)
pdf("test.pdf",w=8,h=6)
p<-ggplot(data=dat3,aes(x=sample,y=value))+geom_boxplot(aes(fill=group))+#绘制分组箱线图
  labs(x=NULL,y="Gene expression",title="Expression of Hub genes in GSE9443")+ #设置y轴标题和图片标题
  scale_fill_manual(values=c("turquoise1","violetred1"))+#设置分组颜色
  stat_compare_means (aes(group=group), #使用R包ggpubr进行统计学检验    
                      label = "p.signif", #设置星标    
                      size = 8, #显著性星标的大小    
                      label.y = c(6,5.5,6.5))+ # 显著性标题放在图中对应y轴的哪个位置  
  theme_classic() +  
  theme(plot.title = element_text(hjust = 0.5, vjust = 0.5), #设置图片标题居中展示    
        text = element_text(size = 15),  #设置图中全部字体的大小    
        axis.text.x = element_text(angle = 75, hjust = 1), #调整横坐标字体倾斜  适用于名儿比较长或者x轴数据多 横着排不开的时候 也会设置字体倾斜    
        legend.position="top",legend.title = element_blank())#图例放在图的上边 不展示图例标题
print(p)
dev.off()

#########Nomogram曲线(数据类型为行为样本，列为基因加status的数据框，均为数字类型)
load("GSE37667_dat_group.Rdata")
dat1=as.data.frame(exp2)
sample_status=group_list
marker=c('PRKAB1','NR4A1','HSPA5','CDKN1A')
dat1=dat1[marker,]
dat1=dat1[,-42]
important_gene_exp=as.data.frame(t(dat1))
status=rep(0,ncol(dat1))
loc1=which(group_list=="SD")
status[loc1]=1
important_gene_exp$status=status


library(rms)
nom_data=important_gene_exp
dd=datadist(nom_data)
options(datadist="dd") 
mod.glm=lrm(status~.,data=nom_data)
library(regplot)
regplot(mod.glm,
        #observation=nom_data[21,],
        points=TRUE,
        odds=F,
        leftlabel=T,
        prfail=TRUE,
        showP=T,
        droplines=T,
        #colors="red",
        rank="range",
        interval="confidence",
        title="Nomogram")  



library(pROC)
mod_roc=roc(nom_data$status~predict(mod.glm),ci=TRUE)
plot(mod_roc,  
     print.auc=TRUE, 
     print.auc.x=0.5, print.auc.y=0.5, 
     auc.polygon=TRUE, 
     auc.polygon.col="skyblue",  
     grid= FALSE, 
     legacy.axes=TRUE,main="GSE9442&GSE33302")  


sig_exp=as.data.frame(predict(mod.glm))
Group=group_list
sig_exp$Group=Group
colnames(sig_exp)=c("Riskscore","Group")
sig_exp$Group=factor(sig_exp$Group,levels = c("Control","SD"))
library(ggprism)
library(ggpubr)
library(ggbiplot)

ggplot(data = sig_exp, aes(x = Group, y =Riskscore, fill = Group)) +
  geom_violin(position = position_dodge(0.5))+
  geom_boxplot(width=0.1,
               position = position_dodge(0.5))+
  theme_prism()+
  scale_fill_manual(values = c("steelblue","#ff7f00"))+ 
  stat_compare_means()

library(rmda)  
fit1 <- decision_curve(status~PRKAB1+NR4A1+HSPA5+CDKN1A,data=nom_data,
                       bootstraps = 50    
)

plot_decision_curve(fit1, curve.names = "Nomogram Model",
                    cost.benefit.axis = T, 
                    confidence.intervals = T
)


simple<- decision_curve(status~CDKN1A+HSPA5+NR4A1+PRKAB1,data=nom_data,
                        
                        family = binomial(link ='logit'),
                        thresholds= seq(0,1, by = 0.01),
                        confidence.intervals = 0.95,
)#研究类型：case-control or cohort 


plot_decision_curve(simple,curve.names=c('Nomogram Model'), 
                    cost.benefit.axis =T,
                    col= c('red'),
                    confidence.intervals=FALSE, 
                    standardize = FALSE)






mod.glm=lrm(status~.,data=nom_data,x=TRUE,y=TRUE)
cal1 <- calibrate(mod.glm, method = "boot", B=1000)
plot(cal1,
     xlim = c(0,1),
     xlab = "Predicted Probability",
     ylab = "Observed Probability",
     legend = FALSE,
     subtitles = FALSE)
abline(0,1,col = "black",lty = 2,lwd = 2)
lines(cal1[,c("predy","calibrated.orig")], type = "l",lwd = 2,col="#ff7f00",pch =16)

lines(cal1[,c("predy","calibrated.corrected")], type = "l",lwd = 2,col="steelblue",pch =16)
legend(0.55,0.35,
       c("Apparent","Ideal","Bias-corrected"),
       lty = c(2,1,1),
       lwd = c(2,1,1),
       col = c("black","#ff7f00","steelblue"),
       bty = "n") # "o"



library(pROC)
load("SD_exp_combined-human.Rdata")
exp2=exp2[,-42]
exp2=as.data.frame(t(exp2))
dat1=cbind(Group,exp2)
hub_gene<-as.data.frame(dat1[,c('CDKN1A','HSPA5','NR4A1','PRKAB1')])
dat2=cbind(Group,hub_gene) 
p1=roc(Group ~ CDKN1A, dat2)
p2=roc(Group ~ HSPA5, dat2)
p3=roc(Group ~ NR4A1, dat2)
p4=roc(Group ~ PRKAB1, dat2)


#绘制ROC曲线，col控制roc曲线的颜色
p1$auc
p2$auc
p3$auc
p4$auc
p5$auc
p6$auc


pdf(file="CDKN1A.pdf",width=5,height=5)
plot(p1,
     print.auc=TRUE,
     auc.polygon=TRUE,
     grid=c(0.1,0.2),
     grid.col=c("green", "red"),
     max.auc.polygon=TRUE, 
     auc.polygon.col="skyblue",
     print.thres=TRUE)
dev.off()

pdf(file="NR4A1.pdf",width=5,height=5)
plot(p2,
     print.auc=TRUE,
     auc.polygon=TRUE,
     grid=c(0.1,0.2),
     grid.col=c("green", "red"),
     max.auc.polygon=TRUE, 
     auc.polygon.col="skyblue",
     print.thres=TRUE)
dev.off()

pdf(file="HSPA5.pdf",width=5,height=5)
plot(p3,
     print.auc=TRUE,
     auc.polygon=TRUE,
     grid=c(0.1,0.2),
     grid.col=c("green", "red"),
     max.auc.polygon=TRUE, 
     auc.polygon.col="skyblue",
     print.thres=TRUE)
dev.off()



library(rms)
load("SD_exp_combined-human.Rdata")
sample_status=group_list
dat1=exp2
dat1=as.data.frame(dat1)


marker=c('NR4A1','HSPA5','CDKN1A')
dat1=dat1[marker,]
important_gene_exp=as.data.frame(t(dat1))
status=rep(0,ncol(dat1))
loc1=which(sample_status=="SD")
status[loc1]=1
important_gene_exp$status=status
nom_data=important_gene_exp
nom_data=nom_data[-42,]
dd=datadist(nom_data)
options(datadist="dd") 

nom_data$CDKN1A <- as.numeric(as.character(nom_data$CDKN1A))
nom_data$HSPA5 <- as.numeric(as.character(nom_data$HSPA5))
nom_data$NR4A1 <- as.numeric(as.character(nom_data$NR4A1))

mod.glm=lrm(status~.,data=nom_data)


library(regplot)
mod.glm=lrm(status~.,data=nom_data,x=TRUE,y=TRUE



library(pROC)
mod_roc=roc(nom_data$status~predict(mod.glm),ci=TRUE)
plot(mod_roc,   
     print.auc=TRUE, 
     print.auc.x=0.5, print.auc.y=0.5, 
     auc.polygon=TRUE, 
     auc.polygon.col="skyblue",  
     grid= FALSE, 
     legacy.axes=TRUE,main="GSE33302&GSE9442")  



library(ggprism)
library(ggpubr)
library(rmda)
sig_exp=as.data.frame(predict(mod.glm))
Group=sample_status
sig_exp$Group=Group
colnames(sig_exp)=c("Riskscore","Group")
sig_exp$Group=factor(sig_exp$Group,levels = c("control","SD"))
ggplot(data = sig_exp, aes(x = Group, y =Riskscore, fill = Group)) +
  geom_violin(position = position_dodge(0.5))+
  geom_boxplot(width=0.1,
               position = position_dodge(0.5))+
  theme_prism()+
  scale_fill_manual(values = c("steelblue","#ff7f00"))+ 
  stat_compare_means()


fit1 <- decision_curve(status~CDKN1A+HSPA5+NR4A1+PRKAB1,data=nom_data,
                       bootstraps = 50 # 
)

plot_decision_curve(fit1, curve.names = "Nomogram Model",
                    cost.benefit.axis = T, 
                    confidence.intervals = T
)

simple<- decision_curve(status~CDKN1A+HSPA5+NR4A1,data=nom_data,
                        
                        family = binomial(link ='logit'),
                        thresholds= seq(0,1, by = 0.01),
                        confidence.intervals = 0.95,
)


plot_decision_curve(simple,curve.names=c('Nomogram Model'), 
                    cost.benefit.axis =T,
                    col= c('red'),
                    confidence.intervals=FALSE, 
                    standardize = FALSE)







mod.glm=lrm(status~.,data=nom_data,x=TRUE,y=TRUE)
par(mar = c(5, 5, 4, 2) + 0.1) 

mod.glm <- lrm(status ~ ., data = nom_data, x = TRUE, y = TRUE)
cal1 <- calibrate(mod.glm, method = "boot", B = 1000)

plot(0, 0, type = "n",
     xlim = c(0, 1), ylim = c(0, 1),
     xlab = "Predicted Probability",
     ylab = "Observed Probability",
     main = "Calibration Curve")


lines(cal1[, "predy"], cal1[, "calibrated.orig"],
      type = "l", lwd = 2, col = "#ff7f00")
lines(cal1[, "predy"], cal1[, "calibrated.corrected"],
      type = "l", lwd = 2, col = "steelblue")


abline(0, 1, col = "black", lty = 2, lwd = 2)

legend("bottomright",
       c("Ideal", "Apparent", "Bias-corrected"),
       lty = c(2, 1, 1),
       lwd = c(2, 2, 2),
       col = c("black", "#ff7f00", "steelblue"),
       bty = "n")

dat2 <- nom_data[,-4]

dat2=cbind(group_list,dat2) 
rt=dat2
dat2$CDKN1A <- as.numeric(as.character(dat2$CDKN1A))
dat2$HSPA5 <- as.numeric(as.character(dat2$HSPA5))
dat2$NR4A1 <- as.numeric(as.character(dat2$NR4A1))


y=colnames(rt)[1]
rt[, y] <- factor(rt[, y], levels = c("Control", "SD"))
rt[, 2] <- as.numeric(as.character(rt[, 2]))

bioCol=c("red","blue","green","yellow")
if(ncol(rt)>3){
  bioCol=rainbow(ncol(rt))}
pdf("ROC(M).pdf",width=5,height=5)
roc1=roc(rt[,y], as.vector(rt[,2]))
aucText=c( paste0(colnames(rt)[2],", AUC=",sprintf("%0.3f",auc(roc1))) )
plot(roc1, col=bioCol[1])

for(i in 3:ncol(rt)){
  roc1=roc(rt[,y], as.vector(rt[,i]))
  lines(roc1, col=bioCol[i-1])
  aucText=c(aucText, paste0(colnames(rt)[i],", AUC=",sprintf("%0.3f",auc(roc1))) )
}

legend("bottomright", aucText,lwd=2,bty="n",col=bioCol[1:(ncol(rt)-1)])

dev.off()