###一、模型变量筛选
#单因素生存分析
#1.载入包
#1-1.批量单因素Cox包
library(survival)
library(plyr)
#2.清理工作环境
rm(list = ls()) 
#3.读入数据
setwd("D:/zhang/Documents/R/training/R-DATA-analysis/NCmajorrevise_code")
aa<- read.csv('Recurrece_metastasis_univariate.csv')
#查看数据性质
str(aa)
#4.构建生存分析的y
y<- Surv(time = aa$time,event = aa$status==1)
#5.Uni_cox_model：批量单因素Cox模型，提取HR，95%CI,P值
Uni_cox_model<- function(x){FML <- as.formula(paste("y~",x))
cox<- coxph(FML,data=aa)
cox1<-summary(cox)
HR <- round(cox1$coefficients[,2],2)
PValue <- round(cox1$coefficients[,5],3)
CI5 <-round(cox1$conf.int[,3],2)
CI95 <-round(cox1$conf.int[,4],2)
Uni_cox_model<- data.frame('Characteristics' = x,
                           'HR' = HR,
                           'CI5' = CI5,
                           'CI95' = CI95,
                           'P' = PValue)
return(Uni_cox_model)}  
#将数据中想进行批量单因素回归的变量写在下面
names(aa)
variable.names<- colnames(aa)[c(3:14)] 
#提取变量及亚变量名。{注，要跟上面一致}
name<-coxph(Surv(time,status==1)~ Gender + Age + Mitotic + Necrosis + Plemorphism + Cellularity + Tumor_size + Tumor_site + Specimen_type + WHO_classification +Radiotherapy + Chemotherapy,data=aa)
name<-summary(name)
names<-data.frame(round(name$coefficients[,2],2))
names
#输出批量单因素结果
Uni_cox<- lapply(variable.names,Uni_cox_model)
Uni_cox<- ldply(Uni_cox,data.frame)
Uni_cox
#调整表格
#增加表格HR(95% CI)
Uni_cox$HR.CI95<-paste0(Uni_cox$HR," (",Uni_cox$CI5,"-",Uni_cox$CI95,")")
#p值改一下格式，p=0的改为<0.001
Uni_cox$P[Uni_cox$P==0]<-"<0.001"
#将单因素回归的第一列变量名删除，替换为提取的变量名
Uni_cox<-Uni_cox[,-1]
rownames(Uni_cox)<-rownames(names)
Uni_cox<-tibble::rownames_to_column(Uni_cox, var ="Characteristics")
#查看表格
result<- data.frame(Uni_cox);result
write.csv(result,"UniCoxSur.csv")


#LASSO回归筛选
rm(list = ls()) 
setwd("D:/zhang/Documents/R/training/R-DATA-analysis/NCmajorrevise_code")
library("glmnet")
library("survival")
library("Matrix")
SFT <- read.csv("SYSUCC101_20230728_5.csv")
SFT <- na.omit(SFT)
x=as.matrix(SFT[,c(3:ncol(SFT))])
y=data.matrix(Surv(SFT$time,SFT$status))
fit = glmnet(x, y, family = "cox", alpha = 1, nlambda = 100)
plot(fit, xvar = "lambda", label = TRUE)
cvfit=cv.glmnet(x, y, family="cox", maxit = 1000)
plot(cvfit)
cvfit $ lambda.min
cvfit $ lambda.1se
coef2 <- coef(cvfit $ glmnet.fit, s = 0.04055804, exact = F)
coef1 <- coef(cvfit $ glmnet.fit, s = 0.08537078, exact = F)
coef1
coef2

#随机生存森林
rm(list=ls()) 
library(randomForestSRC)
library(survival)
library(ranger)
library(ggplot2)
library(dplyr)
library(ggfortify)
library(ggRandomForests)
SFT <- read.csv("SYSUCC101_20230728_5.csv")
rf.model <- rfsrc(Surv(time, status) ~ ., data = SFT,ntree = 10000,tree.err = TRUE, importance = TRUE)
print(rf.model)
plot(rf.model)
rf.model[["importance"]]
rfsrc_error <- gg_error(rf.model)
rfsrc_error <- na.omit(rfsrc_error)
plot(rfsrc_error)
#VIMP
plot(gg_vimp(rf.model))
#minimal_depth
varsel_pbc <- var.select(rf.model)
gg_md <- gg_minimal_depth(varsel_pbc)
plot(gg_md)
#VIMP&minimal_depth
plot(gg_minimal_vimp(gg_md))
+theme(legend.position = c(0.8,0.2))


###Baier score
rm(list = ls()) 
setwd("D:/zhang/Documents/R/training/R-DATA-analysis/NCmajorrevise_code")
###1.加载R包
library(survex)
library(survival)
aa<- read.csv('SYSUCC101_2_20230728.csv')
fml = as.formula(paste('Surv(time,status)~',paste0(colnames(aa)[3:7],collapse = "+")))
# Cox 回归
cox <- coxph(fml,
             data = aa,
             model = TRUE,
             x=TRUE)

# Cox回归模型解释器 
cox_exp <- explain(cox)
###Brier score
# 提取模型数据
y <- cox_exp$y
times <- cox_exp$time
# 计算预测值
surv <- cox_exp$predict_survival_function(cox,
                                          cox_exp$data,
                                          times)
# 计算integrated Brier score
integrated_brier_score(y, surv = surv, times = times)
#计算不同时间点的brier score
brier_score(y,surv = surv,times = times)

###二、构建模型
rm(list = ls())
setwd("D:/zhang/Documents/R/training/R-DATA-analysis/NCmajorrevise_code")
library(openxlsx)
library(survival)
library(lattice)
library(Formula)
library(ggplot2)
library(Hmisc)
library(rms)
#install.packages("SparseM")
library("SparseM")
#####普通列线图
lxrdata <- read.csv("SYSUCC_101_FY71_XH84_XH55_ROC.csv")
dd = datadist(lxrdata)
options(datadist = "dd")
lxrdata $ Ki67 <- factor(lxrdata $ Ki67,levels = c(0,1),labels = c("<454.7 Cells/(mm²)","≥454.7 Cells/(mm²)"))
lxrdata $ MTOR <- factor(lxrdata $ MTOR,levels = c(0,1),labels = c("Other","Damaging_mutation"))
lxrdata $ CD163 <- factor(lxrdata $ CD163,levels = c(0,1),labels = c("<929.3 Cells/(mm²)","≥929.3 Cells/(mm²)"))
lxrdata $ Mitotic <- factor(lxrdata $ Mitotic,levels = c(0,1),labels = c("<4/10HPF","≥4/10HPF"))
f <- cph(Surv(time,status==1)~  Ki67 + MTOR + Mitotic + CD163,data = lxrdata, x = TRUE,y = TRUE,surv = TRUE)
f
validate(f, method = "boot", B = 1000,dxy = T)
###计算C指数方法一
rcorrcens(Surv(time,status==1) ~ predict(f), data = lxrdata)
###计算C指数方法二
fit <- coxph(Surv(time,status==1)~Ki67 + MTOR + CD163 +Mitotic,data = lxrdata)
sum.surv <- summary(fit)
c_index <- sum.surv $ concordance
c_index
Survival <- Survival(f)
Survival1 <- function(x)Survival(60,x)
Survival2 <- function(x)Survival(96,x)
nom <- nomogram(f,fun = list(Survival1,Survival2),
                fun.at = c(.0001,.01,.05,seq(.1,.9,by = .1),.95,.99,.999),
                funlabel = c("5-year survival probability", "8-year survival probability"))
nom
plot(nom, xfrac = 2)

rcorrcens(Surv(time,status) ~ predict(f), data = lxrdata)



#####外部队列验证
be <- read.csv("SYSUCC_101_FY71_XH84_XH55_ROC.csv")
dd1 = datadist(be)
options(datadist = "dd1")
be $ Ki67  <- factor(be $ Ki67 ,levels = c(0,1),labels = c("<454.7 Cells/(mm²)","≥454.7 Cells/(mm²)"))
be $ Mitotic <- factor(be $ Mitotic,levels = c(0,1),labels = c("<4/10HPF","≥4/10HPF"))
be $ CD163 <- factor(be $ CD163,levels = c(0,1),labels = c("<929.3 Cells/(mm²)","≥929.3 Cells/(mm²)"))
be $ MTOR <- factor(be $ MTOR,levels = c(0,1),labels = c("Other","Damaging_mutation"))
pre1<-predict(f, newdata=be)#设定预测值
f1 <- cph(Surv(time, status==1)~pre1,
          x=T, y=T, surv=T, data=be, time.inc=60)#
validate(f1, method="boot", B=1000, dxy=T)
rcorrcens(Surv(time, status) ~pre1, data = be)
####基于regplpot包绘制交互式列线图模
#install.packages("regplot")
library(regplot)
regplot(f,
        plots = c("density", "no plot"),
        observation=be[74,],
        failtime = c(60, 96), prfail = TRUE)

###ROC
rm(list = ls())
setwd("D:/zhang/Documents/R/training/R-DATA-analysis/NCmajorrevise_code")
#instalstate_PFSl.packages("pROC")     # 下载 pROC 包
#install.packages("ggplot2")  # 下载 ggplot2 包
library(pROC)      # 加载pROC包
library(ggplot2)  # 调用ggplot2包以利用ggroc函数
###导入数据
aSAH <- read.csv("SYSUCC_101_FY71_XH84_XH55_ROC.csv")
###建立曲线
roc1 <- roc(aSAH$status,aSAH$Nomogram101);
roc1  # Build a ROC object and compute the AUC
roc2 <- roc(aSAH$status, aSAH$WHO_classification);
roc2  # Create a few more curves for the next examples
roc3 <- roc(aSAH$status, aSAH$mDemicco);
roc3
roc4 <- roc(aSAH$status, aSAH$G_score);
roc4
###统计分析
###分别计算roc1、roc2的AUROC和95%CI
auc(roc1);
ci(roc1);
auc(roc2);
ci(roc2);
auc(roc3);
ci(roc3);
auc(roc4);
ci(roc4);
###plot函数绘制多条曲线
plot(roc1,
     print.auc=TRUE, print.auc.x=0.5, print.auc.y=0.5,
     # 图像上输出AUC值,坐标为（x，y）
     auc.polygon=TRUE, auc.polygon.col="#fff7f7", # 设置ROC曲线下填充色
     max.auc.polygon=FALSE,  # 填充整个图像
     grid=c(0.1, 0.2), grid.col=c("black", "black"),  # 设置间距为0.1，0.2，线条颜色
     print.thres=TRUE, print.thres.cex=0.9, # 图像上输出最佳截断值，字体缩放倍数
     smooth=F, # 绘制不平滑曲线
     main="Comparison of two ROC curves", # 添加标题
     col="#FF2E63",  # 曲线颜色
     legacy.axes=TRUE)   # 使横轴从0到1，表示为1-特异度
###添加roc2曲线
plot.roc(roc2,
         add=T,  # 增加曲线
         col="black", # 曲线颜色为红色
         print.thres=TRUE, print.thres.cex=0.9,  # 图像上输出最佳截断值，字体缩放倍数
         print.auc=TRUE, print.auc.x=0.2,print.auc.y=0.2,
         # 图像上输出AUC值,坐标为（x，y）
         smooth = F)  # 绘制不平滑曲线
###添加roc3曲线
plot.roc(roc3,
         add=T,  # 增加曲线
         col="#003399", # 曲线颜色为红色
         print.thres=TRUE, print.thres.cex=0.9,  # 图像上输出最佳截断值，字体缩放倍数
         print.auc=TRUE, print.auc.x=0.3,print.auc.y=0.3,
         # 图像上输出AUC值,坐标为（x，y）
         smooth = F)  # 绘制不平滑曲线
###添加roc3曲线
plot.roc(roc4,
         add=T,  # 增加曲线
         col="#339933", # 曲线颜色为红色
         print.thres=TRUE, print.thres.cex=0.9,  # 图像上输出最佳截断值，字体缩放倍数
         print.auc=TRUE, print.auc.x=0.4,print.auc.y=0.4,
         # 图像上输出AUC值,坐标为（x，y）
         smooth = F)  # 绘制不平滑曲线

###Survival curve
rm(list=ls()) 
setwd("D:/zhang/Documents/R/training/R-DATA-analysis/NCmajorrevise_code")
library("Rcpp")
library("survival")
library("survminer")
library("dplyr")
library("ggplot2")
library("ggpubr")
SFT <- read.csv("SYSUCC_101_fy71_xh84_xh55_1.csv")
SFT <- na.omit(SFT)
x=as.matrix(SFT[,c(4:ncol(SFT))])
y=data.matrix(Surv(SFT$time,SFT$status))
fit <- survfit(Surv(time,status) ~ Tumor_site,data = SFT)
plot(fit)
p1 <- ggsurvplot(fit)
p1
ggsurvplot(fit,data = SFT,
           risk.table = TRUE,
           pval = TRUE,
           conf.int = F,
           size = 1.2,
           palette = c("#266bb5","#ee1c46","#f47d38","#97cd79","#f177b0"),
           risk.table.col = "strata",
           legend.labs = c("1","2","3","4","5"),
           surv.median.line = "hv",
           title = "PFS",
           ylab = "Progressive-free survival(percentage)",xlab = "Time(months)",
           legend.title = "SYSUCC_101_Nomogram101")


###boxline
#install.packages("colorspace")
library(ggplot2)
library(colorspace)
library(ggpubr)
setwd("D:/zhang/Documents/R/training/R-DATA-analysis/box-line")
data1 <- read.csv("violin.csv",header = T)
p <- ggplot(data1,aes(x = group, y = value)) +
  stat_boxplot(geom = "errorbar", size = 1, width = 0.5, linetype = "solid",
               col = c("#58b485","#97cd79","#bfded6","#cb83ad","#fee540","#f47d38",'#f177b0','#f27278',"#266bb5","#ee1c46")) +
  geom_boxplot(size = 1, fill = "white", linetype = "solid", col = c("#58b485","#97cd79","#bfded6","#cb83ad","#fee540","#f47d38",'#f177b0','#f27278',"#266bb5","#ee1c46"))+
  geom_jitter(width = 0.1,shape = 20,alpha = 1, size = 3) +
  aes(color = group) +
  scale_color_manual(values = c("#58b485","#97cd79","#bfded6","#cb83ad","#fee540","#f47d38",'#f177b0','#f27278',"#266bb5","#ee1c46")) 
#  Add p-value
p + stat_compare_means()
p + stat_compare_means( aes(label = ..p.signif..))

###sankey
setwd("D:/zhang/Documents/R/training/R-DATA-analysis/sankey")
library("ggalluvial")
library("ggplot2")
a <- read.csv("sankey1_CNS.csv")
p1 <- ggplot(a, aes(y = value, axis2 = source, axis1 = target))+
  geom_alluvium(aes(fill = source), width = 0, reverse = F, size = 4, discern = T)+
  geom_stratum(width = 1/8, reverse = F, discern = T)+
  geom_text(stat = "stratum", aes(label = after_stat(stratum)), reverse = F, size = 4, angle = 0, discern = T)+
  scale_x_continuous(breaks = 1:2, labels = c("Source", "Target"))+
  theme(legend.position = "none")
p1


###upset
rm(list=ls()) 
setwd("D:/zhang/Documents/R/training/R-DATA-analysis/upset")
#install.packages("URtools")
library(UpSetR) 
#library(Rtools)
data <- read.csv("upset.csv",header=TRUE)
head(data,6)
upset(data, 
      sets = c("WHO", "Nomogram", "Demicco", "Gscore"),#查看特定的几个集合 
      mb.ratio = c(0.55, 0.45),#控制上方条形图以及下方点图的比例 
      order.by = "freq", #如何排序，这里freq表示从大到小排序展示 
      keep.order = TRUE, #keep.order按照sets参数的顺序排序 
      number.angles = 30, #调整柱形图上数字角度 
      point.size = 2, line.size = 1, #点和线的大小 
      mainbar.y.label = "Genre Intersections", sets.x.label = "Movies Per Genre", #坐标轴名称 
      text.scale = c(1.3, 1.3, 1, 1, 1.5, 1)) #六个数字，分别控制c(intersection size title, intersection size tick labels, set size title, set size tick labels, set names, numbers above bars)