对“不同数据来源的生存分析比较”的补充说明
前面我的学徒的一个推文：不同数据来源的生存分析比较 , 代码细节和原理展现做的非常棒，但是因为学徒的TCGA数据库知识不熟悉，所以被捉到了一个bug，先更正一下：
有留言说：“TCGA里病人01-09是肿瘤，>10是正常的，他没有根据病人的barcode去掉正常组织。“在此向 ta 的提醒表示感谢。
关于 TCGA barcode
简单的描述可以看下面这张图：
如果想更详细地了解，请参考：https://gdc.cancer.gov/resources-tcga-users/tcga-code-tables
回到我的数据
和上次一样，先读取数据并预处理
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
)
最后的结果如下：
上次的结果如下：
比较之下差别还是很大的，以后要多多注意了。