R生存分析|关心的变量KM曲线不显著，还有救吗？

如果想查看某些因素，如年龄，性别，分期，肿瘤数目，大小，实验室指标 或者 通过生信手（tao）段（lu）构建的模型和评分是否对预后有影响时候，经常会把连续变量变为分类变量，然后绘制KM曲线或者列线图等。

这时候会有一些常用的方法：

（1）实验室指标：根据正常范围进行分类

（2）临床指标：根据临床意义进行分类

（3）生信模型评分：根据中位数，平均值等进行分类

（4）生信模型评分：根据统计上的最优cutoff来分类

本次主要介绍基于统计上的最优cutoff分类的方法，并与常见的中位数进行简单的比较。

一 载入数据，R包

为了复现方便，使用内置myeloma数据集

#载入所需的R包
library("survival")
library("survminer")
#查看myeloma数据集
data(myeloma)
head(myeloma)

二 KM-中位数分类

以TP53基因为例，按照常用的中位数分为表达量高，低两组

#按照中位数进行分类
myeloma <- myeloma %>%
 mutate(TP53_cat = ifelse(TP53 > median(TP53) ,"High" , "Low"))

head(myeloma)

构建模型，并绘制KM曲线

#构建模型
fit <- survfit(Surv(time, event) ~ TP53_cat, data = myeloma)

#绘制生存曲线并显示P值
ggsurvplot(fit,
          data = myeloma,
          palette=c("red", "blue"), #自定义颜色
          legend.labs=c("TP53_High","TP53_Low"), #自定义标签
          risk.table = TRUE, 
          break.x.by = 6, #横坐标刻度间隔
          pval = T) #是否显示P值

如图显示P值不显著，这时候可以试一下最优cutoff。

更多调整可参考R|生存分析 - KM曲线 ，必须拥有姓名和颜值

三 KM-最优cutoff分类

3.1 计算最优cutoff

使用surv_cutpoint函数找到最优cutoff

res.cut <- surv_cutpoint(myeloma, 
                        time = "time",
                        event = "event",
                        variables = c("TP53", "WHSC1")) #可以添加多列

summary(res.cut)#查看最佳cutoff

cutpoint statistic

TP53     748.3  2.928906

WHSC1   3205.6  3.361330

可以看到TP53 和 WHSC1 基因统计得到的最优cutoff。

3.2 根据最优cutoff分类

A：根据得到的最优cutoff 自行分类

myeloma <- myeloma %>%
 mutate(TP53_cutoff = ifelse(TP53 > 748.3 ,"High" , "Low"))

head(myeloma)

构建模型，并绘制KM曲线

#构建模型
fit <- survfit(Surv(time, event) ~ TP53_cutoff, data = myeloma)
#绘制生存曲线
ggsurvplot(fit,
          data = myeloma,
          palette=c("red", "blue"), #自定义颜色
          legend.labs=c("TP53_High","TP53_Low"), #自定义标签
          risk.table = TRUE,  
          break.x.by = 6,##横坐标间隔
          pval = T) #是否展示P值

可以看到P值 < 0.05了，变化还是很明显的。

B：根据surv_categorize函数获取重新构建的矩阵

此处推荐这种方法，能比较简单的获取重新构建的矩阵

##重新构建的矩阵
res.cat <- surv_categorize(res.cut)

head(res.cat)

time event TP53 WHSC1

GSM50986 69.24     0  low   low

GSM50988 66.43     0 high   low

GSM50989 66.50     0 high  high

GSM50990 42.67     1 high  high

GSM50991 65.00     0  low   low

GSM50992 65.20     0 high   low

构建模型，并绘制KM曲线

fit <- survfit(Surv(time, event) ~TP53, data = res.cat) 
#绘制生存曲线
ggsurvplot(fit,
          data = res.cat,
          palette=c("red", "blue"),
          legend.labs=c("TP53_High","TP53_Low"), #标签
          risk.table = TRUE,
          break.x.by = 6, ##横坐标间隔
          pval = T)

结果和自行根据最优cutoff，使用ifelse进行分类得到的结果一致，此处不展示了。

参考资料：https://www.rdocumentation.org/packages/survminer/versions/0.4.9

https://www.rdocumentation.org/packages/survminer/versions/0.4.9/topics/surv_cutpoint

◆ ◆ ◆  ◆ ◆

精心整理（含图PLUS版）|R语言生信分析，可视化（R统计，ggplot2绘图，生信图形可视化汇总）