TCGA生存分析理论与R语言实践指南

TCGA生存分析理论与R语言实践指南

引用

CSDN

1.

https://blog.csdn.net/swangee/article/details/141503303

TCGA（The Cancer Genome Atlas）生存分析是生物信息学领域中常用的数据分析方法之一，主要用于研究肿瘤患者的生存时间和预后情况。本文将从生存分析的基本概念出发，结合R语言实践，详细介绍TCGA生存分析的理论和方法。

一、生存分析的基本概念

1. 生存时间：从规定的观察起点至事件终点出现所经历的时间长度。

2. 生存数据：① 完全数据：准确的生存时间；② 删失数据/截尾数据：由于某种原因无法观察到事件终点出现的数据。删失数据的原因包括：1）研究结束时终点事件未发生；2）失访；3）研究对象因其他原因死亡或无法耐受而终止观察。

1号和5号 观察期间内患者死亡，为完全数据；

2号在观察结束后死亡、3号失访、4号始终未观察到终点事件，这3个数据为删失数据。

二、 生存分析的常用指标

1. 中位生存期（Median survival time）：50%个体存活的时间。

2. 生存函数（Survival function）：某时刻观察对象仍然存活的概率，与时间相关，又称积累生存率。

3. 风险函数（Hazard function） ：某时刻观察对象的瞬时死亡率。

三、Kaplan-Meier 生存曲线

Kaplan-Meier法，即计算的是上述概念所提到的积累生存率，是一种单因素生存分析方法，可用于研究1个因素对于生存时间的影响。

• 曲线上的“+”表示删失数据；
• 曲线越高、下降越平缓表示生存率越高或预后越好；
• 曲线越低、下降越陡峭表示生存率越低或预后不佳；
• 蓝色圈表示中位生存期，肿瘤直径>5cm组中位生存期约23月，≤2cm组约45月；

四、生存率的比较

1. 单因素分析方法检验

非参数检验方法使用较为广泛，主要包括log-rank检验和Breslow检验，二者都属于卡方检验。

• log-rank检验给组间死亡的远期差别更大，即log-rank方法对远期差异敏感。
• Breslow检验给组间死亡的近期差别更大，即Breslow方法对近期差异敏感。

当生存曲线有交叉时，提示具有混杂因素的影响，这两种检验方法不适用。此时可以考虑R语言TSHRC包的Two-Stage方法进行检验，可以同时计算Kaplan-Meier和log-rank test的p值。

2. 多因素Cox回归模型检验

（1）cox回归模型的基本表达式：

h(t)即风险函数，指在t时刻观察对象的瞬时死亡率；

X1m是自变量，β1m是对应X的回归系数；

因此风险函数展示了多个自变量X对因变量h(t)的影响。

回归系数 β 与 HR (Hazard Ratio) 有对应关系，经过对数转换，自变量X取特殊值：X=0表示未暴露组，X=1表示暴露组，风险函数可以改写为

（2）风险比（HR, Hazard Ratio）

【举例】分析乳腺癌术后复发时间与是否化疗有关。X=0表示患者未接受化疗，X=1表示患者接受化疗，经过分析X的回归系数β= -0.38。故HR如下：

提示接受化疗的乳癌术后患者复发的风险为不接受化疗患者的0.68倍。

五、生存曲线绘制及假设检验

前文所述，对生存分析的基本概念和方法有了基本的认识，现在可以结合上一篇文章整理的临床数据来进行生存曲线的绘制。

生存曲线绘制的代码如下：

# 加载包，准备绘制KM生存曲线
library(survival)
library(survminer)
library(dplyr)
library(ggplot2)
# 之前代码写的matrix容易与R内置的matrix混淆，因此这里改成了df来存储数据
class(df)
colnames(df)
# 对患者进行分组
a <- select(df,c("month","Race","Gender","Status","Stage","Age","OS")) # 提取分析需要的列
survdata = Surv(time = month, event = Status =='Dead', data = a) # 存储生存数据
a$Age <- factor(cut(a$Age,breaks=c(-Inf,median(a$Age),Inf),
                  labels = c("Young","Old"))) # 按中位年龄对患者进行分层
sfit <- survfit(Surv(time = month, event = Status == 'Dead') ~ Age, data = a)
sfit
# 进行假设检验
# 两组间的检验可以使用survdiff函数
data.survdiff <- survdiff(survdata ~ Age,data = a)
data.survdiff
# 计算P值
p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1) 
p.val
# 分组多可以使用pairwise_survdiff函数进行两两检验
# data.survdiff2 <- pairwise_survdiff(Surv(month,Status)~Race,data = a)
# data.survdiff2
# data.survdiff2$p.value  

运行结果：

绘制简单的生存曲线：

plot1 <- ggsurvplot(sfit, data = a)
plot1  
![](https://wy-static.wenxiaobai.com/chat-rag-image/15099800127829284535)

优化生存曲线，添加p值：

plot2 <- ggsurvplot(sfit,
           data = a,
![](https://wy-static.wenxiaobai.com/chat-rag-image/6556698607226153620)
           pval = T)
plot2  

可以看到按中位年龄分组，两组的生存率没有明显的统计学差异，p=0.065

进一步优化生存曲线，添加95%置信区间、风险表等，让图更加美观！

plot3 <- ggsurvplot(sfit, # 创建的拟合对象
![](https://wy-static.wenxiaobai.com/chat-rag-image/1023332887588531530)
           data = a,  # 指定变量数据来源
           conf.int = TRUE, # 显示置信区间
           pval = TRUE, # 添加P值
           surv.median.line = "hv",  # 添加中位生存时间线
           risk.table = TRUE, # 添加风险表
           risk.table.height=.25,
           xlab = "Time in months",
           ylab = "Overall Survival Probability",
           legend = c(0.8,0.75), # 指定图例位置
           legend.title = "Age", # 设置图例标题“Age”
           legend.labs = c("Young","Old"), # 指定图例分组标签
           break.x.by = 20,  # 设置x轴刻度间距
           palette = "hue")  # 默认为hue，可以自定义调色板 "grey","npg","aaas","lancet","jco", "ucscgb","uchicago","simpsons"和"rickandmorty"
plot3  

还可以添加总生存曲线：

plot4 <- ggsurvplot(sfit, # 创建的拟合对象
           data = a,  # 指定变量数据来源
![](https://wy-static.wenxiaobai.com/chat-rag-image/3421019279601565757)
           conf.int = TRUE, # 显示置信区间
           pval = TRUE, # 添加P值
           surv.median.line = "hv",  # 添加中位生存时间线
           risk.table = TRUE, # 添加风险表
           risk.table.height=.25,
           xlab = "Time in months",
           ylab = "Overall Survival Probability",
           legend = c(0.8,0.75), # 指定图例位置
           legend.title = "Age", # 设置图例标题为“Age”
           legend.labs = c("All","Young","Old"), # 指定图例分组标签
           break.x.by = 20,  # 设置x轴刻度间距
           palette = "hue",  # 默认为hue，可以自定义调色板 "grey","npg","aaas","lancet","jco", "ucscgb","uchicago","simpsons"和"rickandmorty"
           add.all = TRUE) # 添加总患者生存曲线
plot4  

六、Cox回归分析及KM曲线优化

接下来我们利用cox模型计算HR

### cox回归分析计算HR
![](https://wy-static.wenxiaobai.com/chat-rag-image/18030702050041279111)
b <- select(df,c("month","Race","Gender","Status","Stage","Age","OS"))
res.cox <- coxph(Surv(time = month, event = Status == 'Dead') ~ Age, data = b)
coxSummary <- summary(res.cox)
coxSummary  

cox分析结果：HR=1.015（1.002-1.028），似然比检验和Wald检验p均小于0.05。

继续优化KM曲线，将HR添加到图形中：

# 首先绘制生存曲线
plot5 <- ggsurvplot(sfit1, # 创建的拟合对象
           data = a,  # 指定变量数据来源
           conf.int = TRUE, # 显示置信区间
           # pval = TRUE, # 添加P值
           surv.median.line = "hv",  # 添加中位生存时间线
           risk.table = TRUE, # 添加风险表
           risk.table.height=.25,
           xlab = "Time in months",
           ylab = "Overall Survival Probability",
           legend = c(0.8,0.75), # 指定图例位置
           legend.title = "Age", # 设置图例标题“Age”
           legend.labs = c("Young","Old"), # 指定图例分组标签
           break.x.by = 20,  # 设置x轴刻度间距
           palette = "hue")  # 默认为hue，可以自定义调色板 "grey","npg","aaas","lancet","jco", "ucscgb","uchicago","simpsons"和"rickandmorty"
# 修改图例，添加HR
library(ggplot2)
theme <- theme_bw() + 
  theme(plot.title = element_text(hjust = 0.5,size=20),
        axis.text.x = element_text(hjust = 0.5,size=16), 
        axis.text.y = element_text(size = 16),
        axis.title.x = element_text(size = 18), 
        axis.title.y = element_text(size = 18), 
        axis.line = element_line(colour = "black",size = 1), 
        legend.text = element_text(colour = "black",size = 16),
        legend.title = element_blank(),
        legend.position = c(0.65, 0.9))
plot5$plot <- plot5$plot + theme + 
  #在图上添加HR等信息：
  annotate("text", x = max(b$month)/60, y = 0.1,
           label= bquote("HR = 1.015 (1.002-1.028)\t, " * italic(p) * " = 0.00149"),size = 5,hjust = 0) # HR和置信区间来自cox分析的输出结果，p值为log-rank检验的输出结果
plot5$table <- plot5$table + theme + theme(plot.title = element_blank(),
                                     legend.position = "none")
plot5  

后面将继续分享TCGA数据分析的更多内容，未完待续...