对“不同数据来源的生存分析比较”的补充说明

前面我的学徒的一个推文：不同数据来源的生存分析比较 , 代码细节和原理展现做的非常棒，但是因为学徒的TCGA数据库知识不熟悉，所以被捉到了一个bug，先更正一下：

有留言说：“TCGA里病人01-09是肿瘤，>10是正常的，他没有根据病人的barcode去掉正常组织。“在此向 ta 的提醒表示感谢。

关于 TCGA barcode

简单的描述可以看下面这张图：

如果想更详细地了解，请参考：https://gdc.cancer.gov/resources-tcga-users/tcga-code-tables

下面以从 UCSC Xena 上下载的数据为例重新做一次生存分析（其他来源的数据也是一样的做法）

回到我的数据

和上次一样，先读取数据并预处理

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)
# 融合两个数据
for_surv=merge(exp,sul,by="patient")
for_surv$CCR1=as.numeric(for_surv$CCR1)
for_surv$CCL23=as.numeric(for_surv$CCL23)
head(for_surv)

生存分析中用到的数据长下面这个样子

> head(for_surv)
patient     CCR1      CCL23 OS OS.time
1 TCGA-3C-AAAU-01A 1.150008 0.09778235  0    4047
2 TCGA-3C-AALI-01A 1.940841 0.36061270  0    4005
3 TCGA-3C-AALJ-01A 2.319911 0.08252519  0    1474
4 TCGA-3C-AALK-01A 1.671542 0.79132710  0    1448
5 TCGA-4H-AAAK-01A 2.000762 0.18760070  0     348
6 TCGA-5L-AAT0-01A 1.630782 0.82857700  0    1477

第一列就是病人 barcode 了，根据 barcode 的含义把 sample 信息取出来看一下

sample_code = substr(for_surv$patient,14,15)

> table(sample_code)
sample_code
  01   06   11 
1075    7  112 

也就是说这 112 个正常组织的样本要去除

for_surv_tumor = for_surv[as.numeric(sample_code)>=0 
                          & as.numeric(sample_code)<=9,]
for_surv_normal = for_surv[as.numeric(sample_code)>=10 
                          & as.numeric(sample_code)<=19,]
for_surv_control = for_surv[as.numeric(sample_code)>=20 
                           & as.numeric(sample_code)<=29,]

选择 tumor 的数据继续走生存分析流程

library(survminer)
cut <- surv_cutpoint(
  for_surv_tumor,
  time = "OS.time",
  event = "OS",
  variables = c("CCL23", "CCR1")
)

> summary(cut)
       cutpoint statistic
CCL23 0.2994369  2.791561
CCR1  3.2532045  1.159042

dat <- surv_categorize(cut)
library(survival)
library(RTCGA)
myfit <- survfit(Surv(OS.time, OS) ~ CCL23 + CCR1,
               data = dat)
ggsurvplot(
  myfit,
  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
)

最后的结果如下：

上次的结果如下：

比较之下差别还是很大的，以后要多多注意了。

文末友情宣传

强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：

• 全国巡讲全球听（买一得五） ，你的生物信息学入门课
• 生信技能树的2019年终总结  ，你的生物信息学成长宝藏