如果想查看某些因素,如年龄,性别,分期,肿瘤数目,大小,实验室指标 或者 通过生信手(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进行分类得到的结果一致,此处不展示了。