一.对一个连续性变量做生存分析,我们首先需要把变量分成两组,有两组才能有比较,此处需要用到ROC曲线,因为ROC曲线能找到最有意义的一个点。
(要注意的是,因为生存资料的结果涉及生存状态和生存时间两个变量,所有不能用普通的ROC曲线,必须用时间依赖性ROC曲线,也就是survivalROC,这样才能把两个因素都分析进去)
library(survivalROC)#加载包
data<- read.csv('cutoff.csv',
sep = ',',
)#读取数据cutoff文件
nobs=NROW(data) #定义数据集的行数
data$OS<-data$OS*30#把OS这一列都乘以30变为天数
cutoff=1460 #设定为4年生存时间(可根据需要修改)
#生存日期转换为days
#delete"NA" 去掉缺失数据
#data=data[which(data$Status!="NA"),]
head(data)#此时文件为上图
#Fit survival ROC model with method of "KM"
SROC=survivalROC(Stime=data$OS,#生存时间
status=data$censored,#生存状态
marker=data$percentage,#需要分析对生存影响的数据
predict.time=cutoff,#预期时间,一般是5年
method="KM")#生存模型KM
cut.op=SROC$cut.values[which.max(SROC$TP-SROC$FP)]
cut.op#此为约登指数最大的点,截断值
#一般把其作为二分类的分割,把指标分成两组,此处可以实现连续型变量二分类分析对生存的影响
#plot survival ROC#画出时间依赖性ROC曲线
plot(SROC$FP,SROC$TP,type="l",xlim=c(0,1),ylim=c(0,1),
xlab=paste("FP","\n","AUC=",round(SROC$AUC,3)),
ylab="TP",main="4-year survival ROC",col="red")
abline(0,1)#加上对角线
legend("bottomright",c("percentage"),col="red",lty=c(1,1))#在图上加上分析的哪个变量数据
####找出了cut.op=SROC$cut.values[which.max(SROC$TP-SROC$FP)]
二.找到截断值分为两组之后,做生存分析,画图
##找出截断值分为两组后做生存分析,比较两组的生存是否有差异
library(survminer)
sur.cat <- surv_categorize(cut.op)
head(sur.cat)
fit1 <- survfit(Surv(OS, censored) ~percentage, data = sur.cat)
# 以sur.cat为截断值做生存分析
ggsurvplot(fit1, data = sur.cat, conf.int = F,pval = T,legend.title="percentage", legend.labs=c(">22","<=22"))
#画出生存分析图形
最后
感谢jimmy的生信技能树团队!
感谢导师岑洪老师!
感谢健明、孙小洁,慧美等生信技能树团队的老师一路以来的指导和鼓励!