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

Posted on 2020年3月31日 by ulwvfje

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

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

关于 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，多一点数据认知，让他们的科研上一个台阶：

Comments are closed.