不同数据来源的生存分析比较

Merry Christmas

最近浙江大学的学徒咨询了我TCGA数据里面的生存分析的不确定性问题，略微指引了一下他，就让他写了这个教程！

对比2015.11.1的TCGA数据，最新的TCGA数据，GOBO数据三种数据来源的CCR1，CCL23两种基因在乳腺癌病人中的生存分析。

说在前面：

事情的起因是我（浙江大学学徒）在文献中看到这么一副生存分析曲线图（文献题目：CCL9 Induced by TGFb Signaling in Myeloid Cells Enhances Tumor Cell Survival in the Premetastatic Organ）

于是想重复一下，这篇文献的数据来源是GOBO，一个乳腺癌的专属数据库，所以我一开始选择了调用TCGA的数据，但是很可惜这个结果的癌症种类特异性是比较强的，试了几种癌症都没有这么显著的结果，要么就是相反的结果。

不过在曾老师的指引之下我顺便探索了一下不同数据来源的生存分析结果会有什么不同。因为有文献作为参考，所以决定还是以乳腺癌为例子，随便拿两个基因，比如CCR1，CCL23为例来给大家讲解一下我的探索过程！

2015.11.1 TCGA

1.数据获取（RTCGA）

RTCGA是一个可以调用TCGA数据并为画生存分析曲线做方便的数据准备的包，不同于常见的生存分析曲线的地方在于，这个包可以把两个基因的表达信息整合到一起。

除了本文要用到的clinical数据和rnaseq数据外，这个包还支持一系列TCGA数据的调用，但值得注意的是，只能调用2015年11月1日版本的TCGA数据，这是一个比较大的缺点（见下图）。

本文将以乳腺癌和CCL23，CCR1这两种基因的表达信息为例，展示一种癌症、两种基因的生存分析曲线画法。

参考来自原作者的教程：https://github.com/RTCGA/RTCGA/issues/97

2.包的安装

首先需要两个数据包：RTCGA.clinical和RTCGA.rnaseq.

3.数据预处理

rm(list = ls())
options(stringsAsFactors = F)
library(RTCGA.clinical)
library(RTCGA.rnaseq)
library(tidyverse)

# 提取生存情况信息
survivalTCGA(BRCA.clinical) -> BRCA.surv
# 提取两种基因的表达信息
expressionsTCGA(
  BRCA.rnaseq,
  extract.cols = c("CCL23|6368", "CCR1|1230")#这里需要symbol和entrez ID
) -> BRCA.rnaseq
# 看看数据长啥样
head(BRCA.surv); head(BRCA.rnaseq)
# 两个数据包的病人barcode长得不一样，需要统一
BRCA.rnaseq$bcr_patient_barcode=
  substr(BRCA.rnaseq$bcr_patient_barcode, 1, 12)
# 随后就可以融合两个数据
BRCA.surv %>%
  left_join(BRCA.rnaseq,
            by = "bcr_patient_barcode") ->
  BRCA.surv_rnaseq
BRCA.surv_rnaseq=na.omit(BRCA.surv_rnaseq)
head(BRCA.surv_rnaseq)

library(survminer)
# 改个名，方便一些
colnames(BRCA.surv_rnaseq)[5]="CCL23"
colnames(BRCA.surv_rnaseq)[6]="CCR1"
# 找到合适的“高低表达量”的定义，高于cut，就是高表达
BRCA.surv_rnaseq.cut <- surv_cutpoint(
  BRCA.surv_rnaseq,
  time = "times",
  event = "patient.vital_status",
  variables = c("CCL23", "CCR1")
)
summary(BRCA.surv_rnaseq.cut)
plot(BRCA.surv_rnaseq.cut, "CCL23", palette = "npg")
# 将数字的表达信息根据粗体转化为high or low的二元信息
BRCA.surv_rnaseq.cat <- surv_categorize(BRCA.surv_rnaseq.cut)
head(BRCA.surv_rnaseq.cat)

4.画图

library(survival)
# 常规的画法
fit <- survfit(Surv(times, patient.vital_status) ~ CCL23 + CCR1,data = BRCA.surv_rnaseq.cat)
ggsurvplot(
  fit,
  risk.table = TRUE,
  pval = TRUE,
  conf.int = FALSE,
  xlim = c(0,4000),
  break.time.by = 1000,
  ggtheme = theme_RTCGA(),
  risk.table.y.text.col = T,
  risk.table.y.text = FALSE
)

可以看到和文献结果基本一致。不过我这里采取的分组和文献中不完全相同，文献中是把两种基因的表达量整合到一起，而我选择了把所有可能的情况都列入分组。

值得注意的是：两个基因的表达量如何整合，其实是一个值得商榷的问题

最新 TCGA

用UCSC xena 浏览器来下载。

1.数据预处理

rm(list = ls())
options(stringsAsFactors = F)
# 下面的两个数据文件均是手动下载的，select_exp.txt是取了想要的两种基因的数据，因为原数据包含所有基因的表达信息，读进R里非常慢
exp=read.table("select_exp.txt",sep = '\t',header = T)
tmp=t(exp)
exp=data.frame(patient=rownames(tmp)[-1],CCR1=tmp[-1,1],
               CCL23=tmp[-1,2])
# 统一病人编号
exp$patient=gsub('[.]','-',exp$patient)
sul=read.table("TCGA-BRCA.survival.tsv",sep = '\t',header = T)
sul=data.frame(patient=sul$sample,OS=sul$OS,OS.time=sul$OS.time)
# 融合两个数据
tmp=merge(exp,sul,by="patient")
tmp$CCR1=as.numeric(tmp$CCR1)
tmp$CCL23=as.numeric(tmp$CCL23)
head(tmp)

# 依旧是用surv_cutpoint来定义high or low
library(survminer)
cut <- surv_cutpoint(
  tmp,
  time = "OS.time",
  event = "OS",
  variables = c("CCL23", "CCR1")
)
summary(cut)
plot(cut, "CCL23", palette = "npg")

dat <- surv_categorize(cut)
head(dat)

2.画图

library(survival)
fit <- survfit(Surv(OS.time, OS) ~ CCL23 + CCR1,
               data = dat)
ggsurvplot(
  fit,
  risk.table = TRUE,
  pval = TRUE,
  conf.int = FALSE,
  xlim = c(0,4000),
  break.time.by = 1000,
  ggtheme = theme_RTCGA(),
  risk.table.y.text.col = T,
  risk.table.y.text = FALSE
)

用R包cgdsr获取TCGA数据

rm(list = ls())
library(cgdsr)
library(DT)

mycgds = CGDS("http://www.cbioportal.org/")
#test(mycgds)
setVerbose(mycgds, TRUE)
all_TCGA_studies = getCancerStudies(mycgds)
DT::datatable(all_TCGA_studies)
all_gen_pro = getGeneticProfiles(mycgds,'brca_tcga')
all_tables <- getCaseLists(mycgds,'brca_tcga')

# 获得两种基因的表达水平
my_dataset <- 'brca_tcga_rna_seq_v2_mrna'
my_table <- "brca_tcga_all"
exp <- getProfileData(mycgds, c("CCR1","CCL23"), my_dataset, my_table)
exp$patient=rownames(exp)

# 获得临床信息
cli_dat = getClinicalData(mycgds,my_table)
myclinicaldata = cli_dat
DT::datatable(myclinicaldata,
              extensions = 'FixedColumns',
              options = list(                    #dom = 't',
                scrollX = TRUE,
                fixedColumns = TRUE
              ))
cli_dat=myclinicaldata[,c("OS_MONTHS","OS_STATUS")]
cli_dat$patient=rownames(myclinicaldata)

# 融合两种数据并微调
tmp=merge(exp,cli_dat,by="patient")
tmp=na.omit(tmp)
colnames(tmp)[4]="OS.time"
colnames(tmp)[5]="OS"
tmp$OS=ifelse(tmp$OS=="DECEASED",1,0)
tmp$OS.time=tmp$OS.time*30

# 后面的步骤就是一样了

两个数据来源都是和老版本TCGA数据库的结果有些许的差别，但大致的趋势是一致的。

GOBO

最后再用文献的数据来源试试。

我只找到了调用GOBO database的一个网页工具，是直接出图的，当选择左右乳腺癌亚型时，情况是这样的：

可以看到结果并不显著，随后我又看了每个亚型分开的图，其中只有一张比较符合文献，但是也没那么显著：

所以文章可能是对数据进行了更多方面的筛选。

总结

三种数据来源的结果大体趋势一致，但是显著性和一些细节上有差别。

你可能会需要：广州专场（全年无休）GEO数据挖掘课，带你飞（1.11-1.12）和 生信入门课全国巡讲2019收官--长沙站