专栏首页用户7627119的专栏R基于TCGA数据画生存曲线

R基于TCGA数据画生存曲线

生存分析是生物信息医学领域中最常见的一类分析方法。其应用主要包括几个方面:

一是研究某癌症类型中患者的生存情况;

二是研究biomarker在癌症中的预后效能;

三是研究不同分组之间患者的生存是否存在差异。

常见如对于同一种癌症类型使用放疗的患者与使用化疗的患者之间的生存是否存在显著差异,从而判断使用哪种治疗方法更有利于患者的生存。

01

重点概念理解

生存分析不同于其它多因素分析的主要区别点就是生存分析考虑了每个观测出现某一结局的时间长短。因此要顺利的画出生存曲线,首先需要理解两个概念。其一是生存时间,其二是终点事件。弄清楚了这两个概念。我们就开始为画图做准备工作啦。

  • 生存时间:从规定的观察起点开始到某一特定的终点事件发生的这段时间。
  • 终点事件:研究者所关心的特定结局。

02

数据准备

首先从TCGA下载临床数据。从TCGA下载数据有很多方法和教程这里就不多加赘述啦。教程虽然多,但是拿到数据如何处理为生存分析时需要的数据格式呢?上面我们说过生存资料的两个变量:结局事件和生存时间,要想画出生存曲线,至少需要包含这两列数据。下面以肾透明细胞癌KIRC数据为例进行代码实战。TCGA-KIRC.GDC_phenotype.tsv文件从https://gdc.xenahubs.net/download/TCGA-KIRC.GDC_phenotype.tsv.gz获取

kirc.phenotype <- read.csv("TCGA-KIRC.GDC_phenotype.tsv", header = T, sep = "\t")
# 提取肿瘤的表型信息
tumor.sample <- sapply(as.character(kirc.phenotype$submitter_id.samples), function(x){
  number <- as.numeric(unlist(strsplit(unlist(strsplit(as.character(x), split="-"))[4], split = "[A-Z]"))[1])
  if(number<=9){
    id <- as.character(x)
    return(id)
  }
})
tumor.sample.id <- unlist(tumor.sample)
tumor.kirc.phenotype <- kirc.phenotype[match(tumor.sample.id, kirc.phenotype$submitter_id.samples), ]
rownames(tumor.kirc.phenotype) <- tumor.kirc.phenotype$submitter_id.samples
# 获取唯一的肿瘤样本以及表型矩阵
uniq.sample.id <- unique(tumor.kirc.phenotype$submitter_id)
uniq.tumor.kirc.phenotype <- tumor.kirc.phenotype[match(uniq.sample.id, tumor.kirc.phenotype$submitter_id), ]
rownames(uniq.tumor.kirc.phenotype) <- uniq.tumor.kirc.phenotype$submitter_id
# 获取OS time
os.time <- rep(0, nrow(uniq.tumor.kirc.phenotype))
dead.index <- which(uniq.tumor.kirc.phenotype$days_to_death > 0)
alive.index <- which(uniq.tumor.kirc.phenotype$days_to_last_follow_up.diagnoses > 0)
os.time[dead.index] <- uniq.tumor.kirc.phenotype$days_to_death[dead.index]
os.time[alive.index] <- uniq.tumor.kirc.phenotype$days_to_last_follow_up.diagnoses[alive.index]
# 获取终点事件,并用0,1表示
os.index <- which(os.time > 0)
os.time <- os.time[os.index]
status <- rep(0, nrow(uniq.tumor.kirc.phenotype))
dead.index <- which(uniq.tumor.kirc.phenotype$vital_status.demographic == "Dead")
alive.index <- which(uniq.tumor.kirc.phenotype$vital_status.demographic == "Alive")
status[dead.index] <- 0
status[alive.index] <- 1
status <- status[os.index]
uniq.tumor.kirc.phenotype <- uniq.tumor.kirc.phenotype[os.index, ]
# 提取感兴趣的表型信息
interesting.tumor.kirc.data <- data.frame(gender = uniq.tumor.kirc.phenotype$gender.demographic,
                                          age = uniq.tumor.kirc.phenotype$age_at_initial_pathologic_diagnosis,
                                          status = status,
                                          OS = os.time)
rownames(interesting.tumor.kirc.data) <- rownames(uniq.tumor.kirc.phenotype)

(向左滑动查看更多)

03

开始画图

得到OS数据后,接下来就进入到画图阶段啦。让我们跟着下面的代码,一步一步的画出KM曲线。

# step1 加载R包
library(survival)
library(survminer)
# step2 使用Surv()函数创建生存数据对象(生存时间、终点事件)
# step3 再用survfit()函数对生存数据对象拟合生存函数,创建KM(Kaplan-Meier)生存曲线
plot.interesting.tumor.kirc.data <- survfit(Surv(interesting.tumor.kirc.data$OS,event=interesting.tumor.kirc.data$status)~gender,
                                            data=interesting.tumor.kirc.data)
# step4 使用ggsurvplot函数进行画图
ggsurvplot(plot.interesting.tumor.kirc.data ,
           pval = T,
           legend.title = "sex",
           legend="right",
           #font.main = c(16, "bold", "darkblue"),
           #censor.shape=26,
           palette = c("blue","red"),
           censor=F,
           xlab="overall survival",
           title="survival analysis")

(向左滑动查看更多)

下面我们基于M分期来画生存曲线。如果对肿瘤TNM分期还不了解的小伙伴可以参考肿瘤TNM分期

#pathologic_M的生存曲线,三个分期
interesting.tumor.kirc.data <- data.frame(pathologic_M = uniq.tumor.kirc.phenotype$pathologic_M,
                                          age = uniq.tumor.kirc.phenotype$age_at_initial_pathologic_diagnosis,
                                          status = status,
                                          OS = os.time)
rownames(interesting.tumor.kirc.data) <- rownames(uniq.tumor.kirc.phenotype)
index=pathologic_M %in% c("M0","M1","MX")
interesting.tumor.kirc.data=interesting.tumor.kirc.data[index,]
plot.interesting.tumor.kirc.data <- survfit(Surv(interesting.tumor.kirc.data$OS,event=interesting.tumor.kirc.data$status)~pathologic_M,
                                            data=interesting.tumor.kirc.data)
# step4 使用ggsurvplot函数进行画图
ggsurvplot(plot.interesting.tumor.kirc.data ,
           pval = T,
           conf.int = TRUE,
           legend.title = "pathologic_M",
           legend="right",
           surv.median.line = "hv",
           palette = c("blue","red","green"),
           censor=F,
           risk.table=T,
           xlab="overall survival",
           title="survival analysis")

以上就是快速上手画出生存曲线图的全部内容了,你会了吗?感兴趣的小伙伴们也可以动手试一试。

本文分享自微信公众号 - 生信交流平台(gh_d04ce007f7b8),作者:生信交流平台

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

原始发表时间:2020-06-18

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

我来说两句

0 条评论
登录 后参与评论

相关文章

  • R语言中%||%是什么意思?

    不知道大家有没有在R代码中见到过这样的运算符号%||%,看上去有点像或,却又不是。

    生信交流平台
  • 肿瘤研究常用数据库,总有一款适合你

    一个基于TCGA数据库,不需要注册登录即可进行可视化分析的网页,不需要代码,里面内容十分丰富,提供了最详细的TCGA在线分析展示。

    生信交流平台
  • 癌症数据库专题

    想研究肿瘤数据库,可不是只有TCGA、GEO能用。小编今天帮大家总结了一些没那么广为人知,但好用且仍在更新的癌症基因数据库。

    生信交流平台
  • express 与 express-generator

    其实express只是一个框架,那么npm install -g express 也仅仅是安装了这个框架,其作用是自己构建express项目的时候,库已经可以引...

    bear_fish
  • 自动化测试如何解析excel文件?

    自动化测试中我们存放数据无非是使用文件或者数据库,那么文件可以是csv,xlsx,xml,甚至是txt文件,通常excel文件往往是我们的首选,无论是编写测试...

    小老鼠
  • 自动化测试如何解析excel文件?

      自动化测试中我们存放数据无非是使用文件或者数据库,那么文件可以是csv,xlsx,xml,甚至是txt文件,通常excel文件往往是我们的首选,无论是编写测...

    py3study
  • hiphop原理分析1

    Hiphop是Facebook开发一款PHP二进制化的一个工具,最开始是由php转为C++,但是后来发现编译为c++的话,许多的时间会花费在编译代码上面,调试不...

    小小科
  • 【JVM】浅谈双亲委派和破坏双亲委派

    笔者曾经阅读过周志明的《深入理解Java虚拟机》这本书,阅读完后自以为对jvm有了一定的了解,然而当真正碰到问题的时候,才发现自己读的有多粗糙,也体会到只有实践...

    joemsu
  • 数值微分|有限差分法的误差分析

    以(1)为例,分子可能会为0。但是我们不能使h太大,因为这样截断错误将变得过大。为了解决这个矛盾,我们可以采取以下措施:

    fem178
  • React 16 加载性能优化指南(下)

    | 导语 本篇干货是接本周三React 16 加载性能优化指南(上)推文。 关于 React 应用加载的优化,其实网上类似的文章已经有太多太多了,随便一搜就是...

    腾讯NEXT学位

扫码关注云+社区

领取腾讯云代金券