大家在进行生存预后分析时发现结果不显著，是不是当头一棒！两眼一黑！难不成这就代表我们的研究没意义吗？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
小小总结
- survminer包的surv_cutpoint函数
- 基于
maxstat包的maxstat.test函数计算出Maximally Selected Rank Statistics，通过最大化差异来确定最佳阈值。 - 可以同时计算多个变量的最佳截断值。
- 一次只能找到一个最佳截断值，无法同时找到多个截断值。
- 基于
- cutoff包的logrank函数
- 在自由度固定的情况下，计算log-rank卡方值，通过检验分组之间的差异来确定最佳阈值。
- 提供了不同的函数用于计算多种模型的截断点，例如cox、linear、logit等，每种模型都使用对应的统计检验来确定阈值。
- 一次只能找到一个最佳截断值。
- X-tile软件
- 基于log-rank卡方值来确定最佳截断值，将数据分为不同组，可以选择找到一个或两个最佳截断值，以将数据分成两组或三组。
- 一次只能计算一个变量的最佳阈值。
大家自行选择哟！依据个人经验，俺更推荐survminer包的surv_cutpoint函数，它在多数情况下表现都还不戳！个人观点，仅供参考！最终解释权归小蛮要所有！
文末碎碎念
那今天的分享就到这里啦！我们下期再见哟！
最后顺便给自己推荐一下嘿嘿嘿！
如果我的分享对你有用的话，欢迎关注点赞在看转发分享阿巴阿巴阿巴阿巴巴巴！这可是我的第一原动力！
蟹蟹你们的喜欢和支持！！！
啊对！如果小伙伴们有需求的话，也可以加入我们的交流群：一定要知道 | 永久免费的生信交流群终于来啦！
还有兴趣的话，也可以看看我掏心掏肺的干货满满 | 给生信小白的入门小建议 | 掏心掏肺版！绝对干货满满！
如果有小伙伴对付费分析有需求的话，可以看看这里：个性化科研服务 | 付费分析试营业正式启动啦！定制你的专属生信分析！可提供1v1答疑！
入群链接后续可能会不定期更新，主要是因为群满换码或是其他原因，如果小伙伴点开它之后发现，咦，怎么失效啦！不要慌！咱们辛苦一下动动小手去主页的要咨询那里，点击进交流群即可入群！
参考资料
- https://cloud.tencent.com/developer/article/1875291
- https://mp.weixin.qq.com/s/pOgbzHZNQC8z7KdrrrNd1A
- https://mp.weixin.qq.com/s/TB43jWO7CX8_o2eUXWl1Fw
