生存预后不显著？最佳阈值来帮你！| 附完整代码 + 注释

大家在进行生存预后分析时发现结果不显著，是不是当头一棒！两眼一黑！难不成这就代表我们的研究没意义吗？NONONO！别慌！说不定还有救！快来看看最佳阈值能不能捞你一把！

对生存分析感兴趣的小伙伴可以查看：看完不会来揍我 | 生存分析详解 | 从基础概念到生存曲线绘制 | 代码注释 + 结果解读

生存分析不显著怎么办

通常情况下，为了确定一个二分变量（例如基因表达高/低）的最佳阈值，我们可能会使用中位数作为阈值，将样本分为两组，然后对生存曲线进行比较。但有时候使用中位数作为阈值可能并不足以找到显著的差异，这时可以考虑使用最佳阈值方法来寻找更加适合的阈值。来！咱们今天就一起去xio习一下！

老规矩！本文所用到的数据和代码，我已经上传到了GitHub，有需要的话，大家可以在公众号后台回复最佳阈值即可获得存放数据的链接，很多需要注意的问题也会在代码注释中进行详细说明。不过我在分享过程中也会把每一步的输入数据和输出结果进行展示，大家可以作为参考并调整自己的数据格式，然后直接用自己的数据跑，也是没有任何问题的！

# 生存预后不显著？最佳阈值来帮你！# 加载包，没安装的记得安装一下哟！
library(survival)
library(survminer)# 加载数据
best_threshold_data <- read.csv("./best_threshold_data.csv")
head(best_threshold_data)
#            sample risk_score OS OS.time
# 1 TCGA-06-0878-01   2.538694  0     218
# 2 TCGA-26-5135-01   3.736050  1     270
# 3 TCGA-06-5859-01   3.701219  0     139
# 4 TCGA-06-2563-01   3.318001  0     932
# 5 TCGA-41-2571-01   5.102783  1      26
# 6 TCGA-28-5207-01   6.899652  1     343# 数据包含了四列：样本名称（sample）、风险评分（risk_score）、生存状态（OS）和生存时间（OS.time）。# 使用中位数作为阈值来将样本分为高风险组和低风险组
# 如果风险评分大于等于中位数，则标记为"High Risk"，否则标记为"Low Risk"
best_threshold_data$group <- ifelse(best_threshold_data$risk_score >= median(best_threshold_data$risk_score), "High Risk", "Low Risk")# 创建生存对象（Surv对象），包含生存时间和生存状态信息
surv_obj <- Surv(time = best_threshold_data$OS.time, event = best_threshold_data$OS)# 使用survfit()函数拟合生存曲线
surv_fit <- survfit(surv_obj ~ best_threshold_data$group)# 绘制生存曲线
ggsurvplot(surv_fit, data = best_threshold_data, surv.median.line = "hv",pval = TRUE,  # 显示p值xlab = "Time (days)", ylab = "Survival Probability",  # x轴和y轴标签legend.title = "",  # 图例标题为空break.x.by = 1000,  # x轴刻度间隔color = "strata",  # 根据分组着色palette = c("#bc5148", "#3090a1"))  # 自定义颜色

哎！最恶毒的诅咒莫过于“祝你P > 0.05”！！！这不显著可咋整，咱这分析不就没意义了嘛！

且慢！我们的救兵来啦！下面请欣赏它的表演！

最佳阈值选择

目前比较常见的方法有：

• 使用survminer包中的surv_cutpoint函数实现
• 使用cutoff包中的logrank函数实现
• 基于X-Tile软件实现

下面，咱们就挨个介绍！

survminer包的surv_cutpoint函数

# survminer包的surv_cutpoint函数
# 使用survminer包的surv_cutpoint函数来寻找最佳阈值
best_threshold_surv <- surv_cutpoint(best_threshold_data,time = "OS.time",  # 生存时间列名event = "OS",      # 生存事件列名variables = "risk_score",  # 需要寻找阈值的变量列名minprop = 0.3,     # 最小比例，防止找到的阈值过于极端progressbar = TRUE)  # 显示进度条# 查看找到的最佳阈值的摘要统计信息
summary(best_threshold_surv)
#            cutpoint statistic
# risk_score 4.420631  1.748393# 我们这里只计算了一个变量的最佳阈值，还可以计算多个变量的最佳阈值，只需将`variables`参数设为`c("变量1, "变量2", "变量3")`。# 根据找到的最佳阈值对数据进行分组
best_threshold_data <- surv_categorize(best_threshold_surv)# 创建生存对象（Surv对象），包含生存时间和生存状态信息
surv_obj <- Surv(time = best_threshold_data$OS.time, event = best_threshold_data$OS)# 使用survfit()函数拟合生存曲线
surv_fit <- survfit(surv_obj ~ best_threshold_data$risk_score)# 绘制生存曲线
ggsurvplot(surv_fit, data = best_threshold_data, surv.median.line = "hv",pval = TRUE,                    # 显示p值xlab = "Time (days)", ylab = "Survival Probability",  # x轴和y轴标签legend.title = "",             # 图例标题为空break.x.by = 1000,             # x轴刻度间隔color = "strata",              # 根据分组着色palette = c("#bc5148", "#3090a1"))  # 自定义颜色

是不是比上面显著多啦！

cutoff包的logrank函数

# cutoff包的logrank函数# 重新加载数据
best_threshold_data <- read.csv("./best_threshold_data.csv")# 加载cutoff包
library(cutoff)# 使用cutoff包的logrank函数来寻找最佳阈值
best_threshold_surv_2 <- logrank(data = best_threshold_data,time = "OS.time",  # 生存时间列名y = "OS",          # 生存事件列名x = "risk_score",  # 需要寻找阈值的变量列名cut.numb = 1,      # 截点个数n.per = 0.2,       # 分组后每组样本量占总样本量的最小比例y.per = 0.1,       # 分组后每组中较少结果的最小比例p.cut = 0.1,       # p值截断round = 5)         # 保留几位小数# 打印找到的最佳阈值及相关信息
best_threshold_surv_2[order(best_threshold_surv_2$pvalue, decreasing = F), ]
#      cut1      n           n.per     y           y.per  pvalue p.adjust
# 1 4.420631 103/51 0.66883/0.33117 76/47 0.73786/0.92157 0.07617 7.08378 # 根据找到的最佳阈值对数据进行分组
best_threshold_data$Group = ifelse(best_threshold_data$risk_score >= best_threshold_surv_2[order(best_threshold_surv_2$pvalue, decreasing = F), ][1, 1],"High Risk","Low Risk")# 创建生存对象（Surv对象），包含生存时间和生存状态信息
surv_obj <- Surv(time = best_threshold_data$OS.time, event = best_threshold_data$OS)# 使用survfit()函数拟合生存曲线
surv_fit <- survfit(surv_obj ~ best_threshold_data$Group)# 绘制生存曲线
ggsurvplot(surv_fit, data = best_threshold_data, surv.median.line = "hv",pval = TRUE,                    # 显示p值xlab = "Time (days)", ylab = "Survival Probability",  # x轴和y轴标签legend.title = "",             # 图例标题为空break.x.by = 1000,             # x轴刻度间隔color = "strata",              # 根据分组着色palette = c("#bc5148", "#3090a1"))  # 自定义颜色

虽然不如刚刚那个，但也比原来好多啦是不！

X-Tile

这个我就不详细介绍啦！有兴趣的小伙伴们可以去点点点试试！

官网：https://medicine.yale.edu/lab/rimm/research/software/

教程：https://cloud.tencent.com/developer/news/283239

小小总结

1. survminer包的surv_cutpoint函数
   • 基于maxstat包的maxstat.test函数计算出Maximally Selected Rank Statistics，通过最大化差异来确定最佳阈值。
   • 可以同时计算多个变量的最佳截断值。
   • 一次只能找到一个最佳截断值，无法同时找到多个截断值。
2. cutoff包的logrank函数
   • 在自由度固定的情况下，计算log-rank卡方值，通过检验分组之间的差异来确定最佳阈值。
   • 提供了不同的函数用于计算多种模型的截断点，例如cox、linear、logit等，每种模型都使用对应的统计检验来确定阈值。
   • 一次只能找到一个最佳截断值。
3. X-tile软件
   • 基于log-rank卡方值来确定最佳截断值，将数据分为不同组，可以选择找到一个或两个最佳截断值，以将数据分成两组或三组。
   • 一次只能计算一个变量的最佳阈值。

大家自行选择哟！依据个人经验，俺更推荐survminer包的surv_cutpoint函数，它在多数情况下表现都还不戳！个人观点，仅供参考！最终解释权归小蛮要所有！

文末碎碎念

那今天的分享就到这里啦！我们下期再见哟！

最后顺便给自己推荐一下嘿嘿嘿！

如果我的分享对你有用的话，欢迎关注点赞在看转发分享阿巴阿巴阿巴阿巴巴巴！这可是我的第一原动力！

蟹蟹你们的喜欢和支持！！！

啊对！如果小伙伴们有需求的话，也可以加入我们的交流群：一定要知道 | 永久免费的生信交流群终于来啦！

还有兴趣的话，也可以看看我掏心掏肺的干货满满 | 给生信小白的入门小建议 | 掏心掏肺版！绝对干货满满！

如果有小伙伴对付费分析有需求的话，可以看看这里：个性化科研服务 | 付费分析试营业正式启动啦！定制你的专属生信分析！可提供1v1答疑！

入群链接后续可能会不定期更新，主要是因为群满换码或是其他原因，如果小伙伴点开它之后发现，咦，怎么失效啦！不要慌！咱们辛苦一下动动小手去主页的要咨询那里，点击进交流群即可入群！

参考资料

1. https://cloud.tencent.com/developer/article/1875291
2. https://mp.weixin.qq.com/s/pOgbzHZNQC8z7KdrrrNd1A
3. https://mp.weixin.qq.com/s/TB43jWO7CX8_o2eUXWl1Fw