TCGA数据库生存分析

正文

简介

生存分析(Survival analysis)是指根据试验或调查得到的数据对生物或人的生存时间进行分析和推断，研究生存时间和结局与众多影响因素间关系及其程度大小的方法，也称生存率分析或存活率分析。

生存分析适合于处理时间－事件数据,生存时间(survival time)是指从某起点事件开始到被观测对象出现终点事件所经历的时间，如从疾病的“确诊”到“死亡”。

类型

生存时间有两种类型：完全数据(complete data)指被观测对象从观察起点到出现终点事件所经历的时间;截尾数据(consored data)或删失数据，指在出现终点事件前，被观测对象的观测过程终止了。由于被观测对象所提供的信息是不完全的，只知道他们的生存事件超过了截尾时间。截尾主要由于失访、退出和终止产生。生存分析方法大体上可分为三类：非参数法、半参数方法和参数法，用Kaplan-Meier曲线(也称乘积极限法Product limit method)和寿命表法(Life table method)估计生存率和中位生存时间等是非参数的方法，半参数方法指Cox比例风险模型，参数方法指指数模型、Weibull模型、Gompertz模型等分析方法。

非参数法

寿命表时描述一段时间内生存状况、终点事件和生存概率的表格，需计算累积生存概率即每一步生存概率的乘积，可完成对病例随访资料在任意指定时点的生存状况评价。survival包中包括了所有生存分析所必须的函数，生存分析主要是把数据放入Surv object，通过Surv()函数做进一步分析。Surv object是将时间和生存状况的信息合并在一个简单的对象内，Surv(time, time2, event,type=c(‘right’, ‘left’, ‘interval’, ‘counting’, ‘interval2’, ‘mstate’),origin=0)，time为生存时间，time2为区间删失的结束时间，event为生存状况，生存状况变量必须是数值或者逻辑型的。如果时数值型，则有两个选项，0表示删失，1表示终点事件，或者1表示删失，2表示终点事件。如果时逻辑型的，则FALSE表示删失，True表示终点事件。type为删失的类型有右删失、左删失、区间删失、第一类区间删失、第二类区间删失。

实例：

使用一个包含免疫浸润以及复发随访记录的数据，探究该免疫指标具有无复发生存意义

1. setwd('E:\\multi')

2. 
   

3. imm_info <> read.csv('blca_immu.csv',header = T,row.names = 1)

4. 
   

5. 
   

6. library(survival)

7. library(survminer)

8. library(dplyr)

9. require('survival')

10. 
    

11. imm_info <> imm_info %>%

12.  dplyr::select(Lymphocyte_Infiltration_Signature.Score,PFI,PFI.Time)

13. names(imm_info)[1] <> 'infiltra'

整理数据，将代表复发时间，复发状态，以及想预测免疫指标提取出来。

1. imm_info <> na.omit( imm_info)

2. imm_info$infiltra <> ifelse(imm_info$infiltra   > median(imm_info[,'infiltra']),

3.                         'high','low')

4. 
   

5. fit <> survfit(Surv(PFI.Time, PFI) ~ infiltra, data = imm_info)

将免疫指标二分类，以中位值为cutoff，将数据分为高表达和低表达

1. ggsurv <> ggsurvplot(fit, data = imm_info,

2.                     pval = T,

3.                     xlim = c(0,2000),  

4.                     break.time.by = 500,  

5.                     xlab = 'Time in days',

6.                     palette = c('#E41A1C', '#377EB8'))

7. 
   

8. ggsurv <> ggpar( ggsurv,

9.                 font.y  = c(16, 'bold'),

10.                 font.x  = c(16, 'bold'),

11.                 legend = 'top',

12.                 font.legend = c(16, 'bold'))

13. 
    

14. ggsurv

作图

Kaplan-Meier曲线有助于可视化两个分类组之间的生存差异，当你设置参数pval = TRUE时，可以获得的对数秩检验值有助于探讨不同组之间的生存率是否存在差异。 但这并不能很好地评估连续性定量变量的对生存的影响。比如你的某一个node属性取值范围是0-33，这将导致生存曲线图上出现33条生存曲线。如果遇到分组过多或者想要评估多个变量如何协同以影响生存。 例如，比如当希望同时检查种族和社会经济状况对生存的影响时就可能需要换种生存分析方法。

Cox回归

Cox PH回归可以评估分类变量和连续变量的影响，并且可以一次模拟多个变量的影响。 coxph（）函数使用与lm（），glm（）等相同的语法。使用Surv（）创建的响应变量位于公式的左侧，用〜指定。Cox回归大致可以分成单因素和多因素两类。

让我们使用常见的肺癌数据并对性别进行单因素Cox回归分析。

1. > fit <> coxph(Surv(time, status)~sex, data=lung)

2. > fit

3. Call:

4. coxph(formula = Surv(time, status) ~ sex, data = lung)

5. 
   

6.      coef exp(coef) se(coef)     z      p

7. sex -0.531     0.588    0.167 -3.18 0.0015

8. 
   

9. Likelihood ratio test=10.6  on 1 df, p=0.00111

10. n= 228, number of events= 165

exp（coef）列是风险比 - 该变量对危险率的乘法效应（对于该变量的每个单位增加）。因此，对于像性别这样的分类变量，从男性到女性，死亡风险降低约40％。如果翻转coef的正负，比如exp（0.531），您可以将其解释为从女性到男性，死亡危险性增加1.7倍，或者男性死亡率约为每单位时间内女性1.7倍（女性死亡率为每单位时间男性的0.588倍）。

简单起见可以用下列来解释：

• HR = 1：无效

• HR> 1：危险增加

• HR <1：减少危害（保护性） 下一步让我们创建一个模型来分析数据集中的所有变量！="">

然后我们计算多因素cox，我们可以根据P值和exp(coef)评估变量对生存的影响。

1. > fit <> coxph(Surv(time, status)~sex+age+ph.ecog+ph.karno+pat.karno+meal.cal+wt.loss,

2. +              data=lung)

3. > fit

4. Call:

5. coxph(formula = Surv(time, status) ~ sex + age + ph.ecog + ph.karno +

6.    pat.karno + meal.cal + wt.loss, data = lung)

7. 
   

8.               coef exp(coef)  se(coef)     z      p

9. sex       -5.51e-01  5.76e-01  2.01e-01 -2.74 0.0061

10. age        1.06e-02  1.01e+00  1.16e-02  0.92 0.3591

11. ph.ecog    7.34e-01  2.08e+00  2.23e-01  3.29 0.0010

12. ph.karno   2.25e-02  1.02e+00  1.12e-02  2.00 0.0457

13. pat.karno -1.24e-02  9.88e-01  8.05e-03 -1.54 0.1232

14. meal.cal   3.33e-05  1.00e+00  2.60e-04  0.13 0.8979

15. wt.loss   -1.43e-02  9.86e-01  7.77e-03 -1.84 0.0652

16. 
    

17. Likelihood ratio test=28.3  on 7 df, p=0.000192

18. n= 168, number of events= 121

19.   (60 observations deleted due to missingness)

对比

我们尝试在肺数据集上基于年龄属性创建一个分类变量，其中以0,62（平均值）和Inf（无限大）为截断值。 基于截断值我们可以添加labels =选项来标记我们创建的分组，例如，'yong'和'old'。 最后，我们可以将结果分配给肺数据集中的新对象。

1. mean(lung$age)

2. hist(lung$age)

3. ggplot(lung, aes(age)) + geom_histogram(bins=20)

1. cut(lung$age, breaks=c(0, 62, Inf))

2. cut(lung$age, breaks=c(0, 62, Inf), labels=c('young', 'old'))

3. 
   

4. # the base r way:

5. lung$agecat <> cut(lung$age, breaks=c(0, 62, Inf), labels=c('young', 'old'))

6. 
   

7. # or the dplyr way:

8. lung <> lung %>%

9.  mutate(agecat=cut(age, breaks=c(0, 62, Inf), labels=c('young', 'old')))

10. 
    

11. head(lung)

12. > head(lung)

13.  inst time status age sex ph.ecog ph.karno pat.karno meal.cal

14. 1    3  306      2  74   1       1       90       100     1175

15. 2    3  455      2  68   1       0       90        90     1225

16. 3    3 1010      1  56   1       0       90        90       NA

17. 4    5  210      2  57   1       1       90        60     1150

18. 5    1  883      2  60   1       0      100        90       NA

19. 6   12 1022      1  74   1       1       50        80      513

20.  wt.loss agecat

21. 1      NA    old

22. 2      15    old

23. 3      15  young

24. 4      11  young

25. 5       0  young

26. 6       0    old

现在，通过截断值我们将age属性截断成“老”和“年轻”两类患者的曲线，同时探讨两者的生存曲线是否存在一些差异，老年患者的生存几率略差。 但是如下图所示p值 = 0.39时，这说明62岁以下和62岁以下人群的生存率差异不显着。

1. ggsurvplot(survfit(Surv(time, status)~agecat, data=lung), pval=TRUE)

但是，如果我们选择一个不同的切点，比如70岁，这大致是年龄分布的上四分位数的截止值。 结果现在发现P值有意义了。

1. # the base r way:

2. lung$agecat <> cut(lung$age, breaks=c(0, 70, Inf), labels=c('young', 'old'))

3. # or the dplyr way:

4. lung <> lung %>%

5.  mutate(agecat=cut(age, breaks=c(0, 70, Inf), labels=c('young', 'old')))

6. # plot!

7. ggsurvplot(survfit(Surv(time, status)~agecat, data=lung), pval=TRUE)

请记住，Cox回归是分析连续变量在其分布范围内，其中Kaplan-Meier图上的对数秩检验值可以根据您对连续变量的截断值分组而改变。 这两种生存分析方法以不同的方式回答了一个类似的问题：回归模型是在问 年龄对生存的影响是什么？，而生存表法回答的问题是， 组与组之间存在生存差异吗？比如在那些不到70岁的人群和70岁以上的人群？

TCGA

癌症基因组图谱（TCGA）是国家癌症研究所（NCI）和国家人类基因组研究所（NHGRI）之间的合作，收集了33种癌症类型的大量临床和基因组数据。 整个TCGA数据集的基因表达超过2PB，数据类型包括CNV分析，SNP基因分型，DNA甲基化，miRNA分析，外显子组测序和其他类型的数据。 可以在cancergenome.nih.gov上了解有关TCGA的更多信息。 数据现在位于Genomic Data Commons Portal。 有很多方法可以访问TCGA数据而无需实际下载和解析来自GDC的数据。 我们将在下面介绍更多这些内容。 但首先，让我们看一个R包，它提供方便，直接的TCGA数据访问。

1. # Try http:// if https:// doesn't work.

2. source('https://bioconductor.org/biocLite.R')

3. # Install the main RTCGA package

4. biocLite('RTCGA')

5. # Install the clinical and mRNA gene expression data packages

6. biocLite('RTCGA.clinical')

7. biocLite('RTCGA.mRNA')

让我们加载RTCGA包，并使用infoTCGA（）函数获取有关每种癌症类型可用数据类型的一些信息。

1. library(RTCGA)

2. infoTCGA()

RTCGA临床数据的生存分析接下来，让我们加载RTCGA.clinical软件包，并获得一些有关可用内容的帮助。

1. library(RTCGA.clinical)

2. ?clinical

这告诉我们所有可用于每种癌症类型的临床数据集。 如果我们只关注乳腺癌，请查看数据有多大！ 仅此数据中有1098行乘3703列。 我们来看一些变量名。 这里要小心使用View（），数据太大，使用view会可能卡。

我们将使用RTCGA包中的survivalTCGA（）函数从临床数据中提取生存信息。 它通过查看生命状态（死亡或活着）并创建一个时间变量来实现这一点，该变量既可以是死亡天数，也可以是随访天数。 查看有关生存TCGA的帮助以获取更多信息。 您可以为其提供一个临床数据集列表，以及要提取的变量的字符向量。 让我们来看看乳腺癌，卵巢癌和多形性胶质母细胞瘤。 我们只是提取癌症类型（admin.disease_code）。

1. # Create the clinical data

2. clin <> survivalTCGA(BRCA.clinical, OV.clinical, GBM.clinical,

3.                     extract.cols='admin.disease_code')

4. # Show the first few lines

5. head(clin)

6. > head(clin)

7.  times bcr_patient_barcode patient.vital_status admin.disease_code

8. 1  3767        TCGA-3C-AAAU                    0               brca

9. 2  3801        TCGA-3C-AALI                    0               brca

10. 3  1228        TCGA-3C-AALJ                    0               brca

11. 4  1217        TCGA-3C-AALK                    0               brca

12. 5   158        TCGA-4H-AAAK                    0               brca

13. 6  1477        TCGA-5L-AAT0                    0               brca

1. > # How many samples of each type?

2. > table(clin$admin.disease_code)

3. 
   

4. brca  gbm   ov

5. 1098  595  576

1. # Tabulate by outcome

2. xtabs(~admin.disease_code+patient.vital_status, data=clin) %>% addmargins()

3. 
   

4.                  patient.vital_status

5. admin.disease_code    0    1  Sum

6.              brca  994  104 1098

7.              gbm   149  446  595

8.              ov    279  297  576

9.              Sum  1422  847 2269

现在让我们针对疾病代码运行Cox PH模型。 默认情况下，它会将乳腺癌作为基线，因为按字母顺序排在第一位。 但是如果你想使用因子（），你可以重新排序。

1. > coxph(Surv(times, patient.vital_status)~admin.disease_code, data=clin)

2. Call:

3. coxph(formula = Surv(times, patient.vital_status) ~ admin.disease_code,

4.    data = clin)

5. 
   

6.                        coef exp(coef) se(coef)    z      p

7. admin.disease_codegbm  2.887    17.948    0.113 25.6 2e-16

8. admin.disease_codeov   1.547     4.697    0.115 13.4 2e-16

9. 
   

10. Likelihood ratio test=904  on 2 df, p=0

11. n= 2269, number of events= 847

这告诉我们，与基线brca组相比，GBM患者的危害增加了约18倍，卵巢癌患者的生存率降低了约5倍。 让我们创建一个生存曲线，用Kaplan-Meier图显示它，并显示前5年存活率的表格。

1. sfit <> survfit(Surv(times, patient.vital_status)~admin.disease_code, data=clin)

2. summary(sfit, times=seq(0,365*5,365))

3. 
   

4. Call: survfit(formula = Surv(times, patient.vital_status) ~ admin.disease_code,

5.    data = clin)

6. 
   

7.                admin.disease_code=brca

8. time n.risk n.event survival std.err lower 95% CI upper 95% CI

9.    0   1096       0    1.000 0.00000        1.000        1.000

10.  365    588      13    0.981 0.00516        0.971        0.992

11.  730    413      11    0.958 0.00851        0.942        0.975

12. 1095    304      20    0.905 0.01413        0.878        0.933

13. 1460    207       9    0.873 0.01719        0.840        0.908

14. 1825    136      14    0.799 0.02474        0.752        0.849

15. 
    

16.                admin.disease_code=gbm

17. time n.risk n.event survival std.err lower 95% CI upper 95% CI

18.    0    595       2   0.9966 0.00237       0.9920       1.0000

19.  365    224     257   0.5110 0.02229       0.4692       0.5567

20.  730     75     127   0.1998 0.01955       0.1649       0.2420

21. 1095     39      31   0.1135 0.01617       0.0858       0.1500

22. 1460     27       9   0.0854 0.01463       0.0610       0.1195

23. 1825     12       9   0.0534 0.01259       0.0336       0.0847

24. 
    

25.                admin.disease_code=ov

26. time n.risk n.event survival std.err lower 95% CI upper 95% CI

27.    0    576       0    1.000  0.0000        1.000        1.000

28.  365    411      59    0.888  0.0139        0.861        0.915

29.  730    314      55    0.761  0.0198        0.724        0.801

30. 1095    210      59    0.602  0.0243        0.556        0.651

31. 1460    133      49    0.451  0.0261        0.402        0.505

32. 1825     78      39    0.310  0.0260        0.263        0.365

1. ggsurvplot(sfit, conf.int=TRUE, pval=TRUE)

从上面图可以看出，不同肿瘤的生存曲线明显存在差异。最好的乳腺癌，其次卵巢癌，最差的是胶质母细胞瘤。