R_生存分析

基本知识

生存分析：研究各个因素与生存时间有无关系以及关联程度大小。可拓展到疾病复发时间，机器的故障时间等。
起始事件(start point)：反应研究对象开始生存过程的起始特征事件。
终点事件(死亡事件, end point)：出现研究者所关心的特定结局。如“病人因该疾病死亡”。
观察时间：从研究开始观察到研究观察结束的时间。
生存时间(survival time)：观察到的存活时间，用符号t表示。
完全数据：从观察起点到死亡事件所经历的时间，生存时间是完整的。
截尾数据（删失值）：观察时间不是由于终点事件而结束的，而是由于(1)失访(2)死于非研究因素(3)观察结束而对象仍存活以上三种原因结束的。常在截尾数据的右上角放一个“+”表示其实该对象可能活的更久。

变量：生存分析的变量有两个：生存时间t和结局变量(0-1)。其中结局变量1表示死亡事件，0表示截尾。
生存曲线：以生存时间为横轴、生存率为纵轴绘制一条生存曲线。
中位生存时间：又称为生存时间的中位数，表示刚好有50%的个体其存活期大于该时间。

生存分析研究的部分内容

1、描述生存过程：研究生存时间的分布特点，估计生存率及其标准误、绘制生存曲线。常用方法：乘积极限法(K-M法)和寿命表法，此处简单介绍K-M法。
注：生存曲线为单因素分析（两个或者多个水平），用中位生存时间表示生存时间的平均水平；

2、比较生存过程：获得生存率及其标准误的估计值后，可以进行两组或多组生存曲线的比较，常用方法有对数秩检验(log-rank检验)。如曲线交叉，可能存在混杂因素。
注：两个或者多个生存曲线的比较（单因素两个或者多个水平）。

3、影响生存时间的因素分析：Cox比例风险回归模型等。

R与生存分析

1. 数据：使用文章Comprehensive molecular comparison of BRCA1 hypermethylated and BRCA1 mutated triple negative breast cancers中的生存分析数据。

2. 目的：还原生存分析结果图，文章中Fig3g的第三张图：

   
   Fig3g 第三张
3. 实现代码：

### R packages
library(survival)
library(survminer)
library(openxlsx)
### Data
pdata82 = openxlsx::read.xlsx("~/41467_2020_17537_MOESM4_ESM.xlsx", 2)
colnames(pdata82) = pdata82[1,]
pdata82 = pdata82[-1,]
### extract the data used to plot (patients with OS/IDFS, treatment: "Adjuvant")
surdata1 = pdata62[,c(1:2,41:42)]
colnames(surdata1) = c("Sample", "class", "IDFS", "Time")
surdata1$class = ifelse(surdata1$class == "BRCA1null", 0, 1) #convert the 'chr' into 'num'
surdata1$Time = as.numeric(surdata1$Time)
surdata1$IDFS = as.numeric(surdata1$IDFS)
### univariate cox
res.cox1 = coxph(Surv(Time,IDFS) ~ class, data = surdata1 ) 
summary(res.cox1)
### Kaplan-Meier fit
km = survfit(Surv(Time,IDFS) ~ class, data = surdata1)
plot(km, main = "Simple plot")
survminer::ggsurvplot(km,conf.int =F, risk.table =F, 
                      xlab = "Time", 
                      ylab = "IDFS",
                      ncensor.plot = FALSE, 
                      pval.method = T,
                      surv.plot.height = 0.1,
                      legend = c(0.2, 0.1),
                      legend.title = "Class", 
                      palette = c("#000000", "#0000CD"),
                      censor.shape = "|",
                      censor.size = 3,
                      legend.labs = c( "BRCA1-null n = 19","BRCA1-hypermethylated n = 43"),
                      pval = "Log-rank test p = 0.29")

结果展示

4. 扩展：多因素分析，log-rank

### multivariate cox
res.cox = coxph(Surv(Time, survivalEvent) ~ 因素1 + 因素2 + 因素3 + 因素4,  data = data ) 
summary(res.cox)
### test the model
cox.zph(res.cox)  
survminer::ggcoxzph(cox.zph(res.cox))

### comparison analysis
### compute the p value using log-rank test
survdiff(Surv(Time, IDFS) ~ class, data = surdata1, rho = 0) # rho = 0 use long-rank test or Mantel-Haenszel test
(p = 0.29 )

---

• R: survival analysis