R基于TCGA数据画生存曲线

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

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

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

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

前面小编给大家介绍过

☞ 生存曲线怎么看，怎么画

☞ 零代码生存曲线—ENCORI篇

☞ 零代码生存曲线—GEPIA篇

☞ 零代码发文神器—UALCAN

今天我们来看看怎么利用R来绘制生存曲线。

1. 重点概念理解

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

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

2. 数据准备

首先从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)

3. 开始画图

得到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")

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

R基于TCGA数据画生存曲线