有的时候我们会发现,单因素cox回归某个基因对生存有意义,当以中位值绘制生存曲线的时候,却变得没意义了。
这是因为中位值并不是这个基因最合适的截断值,因此,我们需要选择最合适的截断值。
选择截断值的软件目前比较火比较方便的便是X-tile软件。
X-tile 软件
X-tile软件长这样:
如果把x-lite软件得到的截断值写进论文,需要引用人家的原始论文,论文在软件官网可以查到。官网地址在软件开始页面左上角。
点击 Analyze 即可开始分析。
可以看到,分析所需要的数据分为三种:Censor,Survival Time,Marker。
接下来准备输入数据:
数据需要三列,列名没有要求。Censor一列为生存状态,生存为0,死亡为1。Survival Time一列为生存时间,单位为日、周、月、年均可。Marker为我们要研究的变量,这里以基因表达量为例。
数据的格式为txt(文本文件,制表符分隔)。相关代码如下:
write.table(a, file = "1.txt", row.names = F, sep = "\t")
准备好数据之后,选择文件开始分析。
将左侧的各项Load到右侧。
然后依次点击如下图:
得到如下五张图:
最上面一张图是各种计算值。
中间左侧一张图有一个大的三角形和一个矩形,三角形横轴表示低表达人数逐渐增加,纵坐标表示高表达人数逐渐增加。点击三角形任意位置可将总人数分成差异最大的三组,点击矩形的任意位置可将总人数分成差异最大的两组。分组的截断值展示在中间的条形图中。其生存曲线在最右侧的图中。
点击三角形的结果如下:
点击矩形的结果如下:
把cut-off值记下来,然后用R语言做生存曲线。
对于生存数据寻找截断值,上述软件可快速找到最优截断值。图中的截断值为0.8。
survminer::surv_cutpoint确定截断值
需要的生存数据同X-tile,三列数据:time为年、周、月、日均可,event为生存状态,生存为0,死亡为1,variables为变量,如基因表达量。
tmp <- surv_cutpoint(meta, time = "time", event = "event", variables = "gene") # cut_off <- as.numeric(tmp[["cutpoint"]][1, 1])
上述与X-lite相同的数据,survminer::surv_cutpoint得到的截断值为0.837。
tinyarray::point_cut确定截断值
生信技能树写的R包,用的也是survminer::surv_cutpoint函数,用了for循环批量得到基因的截断值。
这个R包中其他函数也有 cut.point = T 的参数,可以直接用。survivalROC::survivalROC确定截断值,并画出ROC曲线
需要的生存数据同X-tile,三列数据:time可谓年、周、月、日均可,event为生存状态,生存为0,死亡为1,variables为变量。
cutoff <- 60 # 随访截止时间,3年,5年或10年等,具体看time的格式是日、月、还是年。 Mayo4.2 = survivalROC(Stime = meta$month, # 生存时间,可以是日、月、年 status = meta$status, # 生存状态,生存为0,死亡为1 marker = meta$LINC01602, # 研究变量,如基因表达量 predict.time = cutoff, # 预测截止时间 method = "KM") cut.op = Mayo4.2$cut.values[which.max(Mayo4.2$TP - Mayo4.2$FP)] #计算最佳截断值 plot(Mayo4.2$FP, Mayo4.2$TP, type="l", xlim = c(0,1), ylim = c(0,1), xlab = paste( "FP", "\n", "AUC = ", round(Mayo4.2$AUC, 3)), ylab = "TP", main = "Mayoscore 4, Method = KM \n Year = 1") # 画ROC曲线,如果需要的话 abline(0,1) # 给ROC曲线加辅助线 legend("bottomright",c("Xgene expression"),col="red",lty=c(1,1))
上述与X-lite相同的数据,survivalROC::survivalROC得到的截断值为0.6762。
后续可以按照前述方法整理数据进行生存分析等。同时我们也使用R软件画出了survivalROC曲线,计算了曲线下面积(AUC):
survivalROC曲线,AUC=0.663。