专栏首页生信技能树生存分析就是一个任人打扮的小姑凉

生存分析就是一个任人打扮的小姑凉

最近接到粉丝的TCGA分析需求,想看看指定基因在指定癌症是否具有临床意义(也就是生存分析是否有统计学显著效果咯!)其实很早以前我在生信技能树就号召粉丝讨论过这个问题:集思广益-生存分析可以随心所欲根据表达量分组吗 这里我做实力演绎一下。

我这里选择最方便的 网页工具:https://xenabrowser.net/heatmap/ 选择合适的数据集及样本信息还有基因来演示一下,随便选择一个基因一个癌症吧,如下:

这个时候,我草率的制作了生存分析图如下:

(@ο@) 哇~非常显著,差点准备交差了,然后下意识的看了看病人数量,TCGA数据库的BRCA病人没道理居然快1200个了,肯定是有什么地方错误了,重新看了看,的确是因为没有顾虑到里面有正常组织测序的那些病人,怎么说呢,相当于把有正常组织测序的那一百多个病人,在我这个生存分析里面计算了两次,他们的生存时间信息,生存状态都重复计算了,所以实际上这个生存分析是错误的。

过滤一下,仅仅是保留tumor的表达量信息和病人临床信息,再次制作生存分析曲线,如下所示:

可以看到,之前明明是显著的结果消失了,而且不管是使用哪种表达量划分方式,都达不到统计学显著阈值。

是不是就没有办法了呢?

当然不是,还可以使用R包,一个非常棒的外国小哥博客写的很清楚:http://r-addict.com/2016/11/21/Optimal-Cutpoint-maxstat.html

还有专门的文章,这里就不细心讲解啦。

使用survminer包的surv_cutpoint函数找寻最近生存分析阈值

外国小哥博客写的很清楚:http://r-addict.com/2016/11/21/Optimal-Cutpoint-maxstat.html 我们现在就测试一下这个流程。

首先下载我们前面的数据文件:'PLEKHA5-BRCA.tsv' 内容如下:

总共6列,在前面的 网页工具:https://xenabrowser.net/heatmap/ 选择对应的信息下载即可:

然后是R代码读入上面的文件,主要是列名需要保证正确无误!!!

rm(list=ls())
options(stringsAsFactors = F)
# install.packages("survminer")
library(survminer)
a=read.table('PLEKHA5-BRCA.tsv',header = T,sep = '\t')
head(a)
dat=a[a$sample_type.samples=='Primary Tumor',4:6]
head(dat)
dat$vital_status.demographic=ifelse(dat$vital_status.demographic=='Alive',1,0)
surv_rnaseq.cut <- surv_cutpoint(
  dat,
  time = "OS.time",
  event = "vital_status.demographic",
  variables = c("ENSG00000052126.13")
)
summary(surv_rnaseq.cut)
plot(surv_rnaseq.cut, "ENSG00000052126.13", palette = "npg")

surv_rnaseq.cat <- surv_categorize(surv_rnaseq.cut) 

library(survival)
fit <- survfit(Surv(OS.time, vital_status.demographic) ~ ENSG00000052126.13,
               data = surv_rnaseq.cat)
ggsurvplot(
  fit,                     # survfit object with calculated statistics.
  risk.table = TRUE,       # show risk table.
  pval = TRUE,             # show p-value of log-rank test.
  conf.int = TRUE,         # show confidence intervals for 
  # point estimaes of survival curves.
  xlim = c(0,3000),        # present narrower X axis, but not affect
  # survival estimates.
  break.time.by = 1000,    # break X axis in time intervals by 500. 
  risk.table.y.text.col = T, # colour risk table text annotations.
  risk.table.y.text = FALSE # show bars instead of names in text annotations
  # in legend of risk table
)

重要的的列名是:

  time = "OS.time",
  event = "vital_status.demographic",
  variables = c("ENSG00000052126.13")

如果是你自己的数据集,需要稍微修改哦。

见证奇迹的时刻:

是不是统计学显著啦!!!

函数帮我们选择的分组;

本文分享自微信公众号 - 生信技能树(biotrainee),作者:生信技能树

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

原始发表时间:2019-10-23

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

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 转录本定量本来就不是一件容易的事情

    比如文章:2015 Nov 7. doi: 10.1080/01621459.2015.1040880 就尝试图解如下:

    生信技能树
  • 我测试了一下Jbrowse的安装及初步试用

    前些天我们公众号元老,熊,投稿了关于Jbrowse的史上最全介绍,如下: 可能是最全的JBrowse基因浏览器介绍(请点击阅读) 最为生物信息学痴的我当然不能...

    生信技能树
  • 全网第一个单细胞转录组数据分析实战视频教程

    回首年前开创的单细胞天地公众号,再看看单细胞转录组知识星球的精华资源,一年时间就这样过去了,感慨万千!

    生信技能树
  • 登录微软账号的Windows电脑如何远程?

    一般情况下,我们都使用的是Windows电脑的本地账户。但是随着Windows 10的推广,现在微软也开始主推微软账号登录Windows电脑了。

    BigYoung小站
  • Java中Synchronized的用法

    版权声明:署名,允许他人基于本文进行创作,且必须基于与原先许可协议相同的许可协议分发本文 (Creative Commons)

    Fisherman渔夫
  • 深度剖析Swagger原理swagger简介

    swagger确实是个好东西,可以跟据业务代码自动生成相关的api接口文档,尤其用于restful风格中的项目,开发人员几乎可以不用专门去维护rest api,...

    Java架构
  • 这家创业企业能帮助你投资印度艺术品,房地产和对冲基金

    ? 当人们都投资于快速发展且盈利丰厚的领域时,就会产生泡沫,如果情况发生变化,泡沫很容易就会崩溃。像印度这样的发展中经济中,情况正是如此,很多资金正蜂拥而来。...

    点滴科技资讯
  • 快速掌握R语言中类SQL数据库操作技巧

    在数据分析中,往往会遇到各种复杂的数据处理操作:分组、排序、过滤、转置、填充、移动、合并、分裂、去重、找重、填充等操作。这时候R语言就是一个很好的选择:R可以高...

    1480
  • MVC官方教程索引

    最近一直在学习MVC(MVC出来这么久了才开始学习,惭愧!不过我一向认为MS的东西不到RC版或至少第三个版本,基本上学了也是白学,按微软的风格,这个补丁那个bu...

    菩提树下的杨过
  • 基于滑动场景解析RecyclerView的回收复用机制原理

    最近在研究 RecyclerView 的回收复用机制,顺便记录一下。我们知道,RecyclerView 在 layout 子 View 时,都通过回收复用机制来...

    请叫我大苏

扫码关注云+社区

领取腾讯云代金券