专栏首页生信补给站R生存分析|关心的变量KM曲线不显著,还有救吗?

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 statisticTP53     748.3  2.928906WHSC1   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 WHSC1GSM50986 69.24     0  low   lowGSM50988 66.43     0 high   lowGSM50989 66.50     0 high  highGSM50990 42.67     1 high  highGSM50991 65.00     0  low   lowGSM50992 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

本文分享自微信公众号 - 生信补给站(Bioinfo_R_Python),作者:生信补给站

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2021-08-25

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • TCGA的28篇教程- 对TCGA数据库的任意癌症中任意基因做生存分析

    生存分析一般来说是针对RNA表达数据,可以说mRNA-seq的转录组数据,也可以说miRNA-seq数据,或者基因表达芯片的表达量值。

    生信技能树
  • 纯生信文章补几张免疫组化真的很重要!

    今天和大家分享的是2020年发表在Journal of cancer(IF:3.565)上的一篇文章,“Genome-wide Analysis of the ...

    科研菌
  • 生存分析凭什么不需要矫正P值

    虽然生存分析如此重要而且如此常见,但是仍然有一些未解之谜,不同数据库来源,病人的不同时期的记录信息,以及不同的阈值分组,拿到的结果居然是可以不一样的!虽然大家都...

    生信技能树
  • 纯生信免疫微环境末班车

    今天和大家分享的是2020年2月发表在Aging(IF:4.831)上的一篇文章,“Profiles of immune cell infiltration a...

    科研菌
  • 精心整理(含图PLUS版)|R语言生信分析,可视化

    为了能更方便的查看,检索,对文章进行了精心的整理(PLUS)。建议收藏,各取所需,当前没用也许以后就用到了呢!

    西游东行
  • 二十年前做科研你只需要检测一些基因在一些癌症细胞系表达量情况即可

    检测的基因是 CUL-5: 以及 other CUL family members (CUL-1, -2, -3, -4A, and -4B) 。实际上哺乳动...

    生信技能树
  • Oh my god!不做实验也能发3分SCI!

    大家好,本期给大家推荐的文献是Differentially Expressed lncRNAs in Gastric Cancer Patients: A Po...

    百味科研芝士
  • 生存分析|知道这些又没有坏处

    生存分析:研究各个因素与生存时间有无关系以及关联程度大小。可拓展到疾病复发时间,机器的故障时间等。

    西游东行
  • R|生存分析(1)

    生存分析:研究各个因素与生存时间有无关系以及关联程度大小。可拓展到疾病复发时间,机器的故障时间等。 起始事件:反应研究对象开始生存过程的起始特征事件。 终点事件...

    西游东行
  • 用R语言进行KM生存分析

    R是数据分析常用的软件之一,通过各种功能强大的R包,可以简单方便的实现各种分析。在R语言中,能够进行生存分析的R包很多,survival和survminer是其...

    生信修炼手册
  • 肿瘤panel测序研究不应该公开基因列表吗

    数据分析我们一般希望是从fastq的测序数据文件开始,但是因为并不是常规肿瘤外显子,所以使用agilent的v6不管用,很多流程都需要其panel对应的个性化的...

    生信技能树
  • 不容错过的6分+预后模型套路

    题目:Development of a four-gene prognostic model for clear cell renal cell carcino...

    百味科研芝士
  • 复现生信论文36.临床意义

    生信论文36是单基因分析的生信论文,单纯生信数据库的数据分析,没有湿实验验证,但是可以发表在接近5分的期刊上,很多分析做得很棒,值得借鉴。我们对文章数据进行复现...

    芒果先生聊生信
  • 生存分析详细解读

    探究变量之间的关系是数据挖掘中的一个基本分析内容,对于常规的离散型或者连续型变量,有很多的方法可以用于挖掘其中的关系,比如线性回归,逻辑回归等等。然而有一类数据...

    生信修炼手册
  • 分析不够思路来凑!好的idea简单分析也可以4+分

    大家好,今天茶叶蛋和大家分享的是一篇4+分的学习笔记;文末点击阅读原文可获得原文笔记。

    科研菌
  • 单基因简单生信分析加上3张免疫组化发4+

    大家好,今天和大家分享的是一月份发表在Frontiers in Oncology (IF:4.848)杂志上的一篇文章,“Prognostic Heteroge...

    科研菌
  • R语言使用限制平均生存时间RMST比较两条生存曲线分析肝硬化患者

    在比较性的纵向临床研究中,主要终点往往是发生特定临床事件的时间,如死亡、心衰住院、肿瘤进展等。_风险_比例估计值几乎被常规用于量化治疗差异。然而,当基础模型假设...

    拓端
  • 为何要劳民伤财做同样的数据

    本来呢,还在奇怪,TCGA数据库里面的乳腺癌患者的放化疗信息应该是没有那么全吧。等我看完摘要才明白,原来是研究者自己招募的病人队列,来自于Iceland bet...

    生信技能树
  • 25张图带你玩转表达量差异分析思路

    研究者首先做了一个bulk转录组,走了标准的差异分析,拿到了上下调基因以及注释它们的功能。然后把这些基因在自己的单细胞转录组各个亚群具体看其是否有表达差异,发现...

    生信技能树

扫码关注云+社区

领取腾讯云代金券