零基础入门转录组数据分析——预后模型之多因素cox模型

---

---

1. 预后模型和多因素cox模型基础知识

1.1 预后模型是什么？
预后模型（Prognostic Model）是一种基于统计学和机器学习方法的预测工具，它通过分析患者的临床特征、生物标志物、治疗方案等多种因素，来预测患者未来发生某种不良事件（如疾病复发、死亡等）的概率。

1.2 预后模型的特点是什么？

• 时间维度： 预后模型通常具有时间维度，即预测的是未来某一时刻或时间段内发生的事件，例如：在1000天的时候疾病患者的生存的概率。
• 多因素分析： 预后模型考虑多种影响因素，包括患者的临床特征、生物学标志物、治疗方案等，以提高预测的准确性和可靠性。

1.3 构建预后模型对于临床应用有哪些帮助？

• 指导临床决策： 预后模型可以帮助医生预测疾病的可能进展，从而为患者提供更加个性化的治疗建议。例如，在癌症治疗中，预后模型可以预测患者的复发风险和生存时间，帮助医生制定更加合理的治疗方案。

• 疾病风险评估： 预后模型能够识别哪些患者更有可能出现不良结局，使医生可以进行早期干预，降低患者的风险。

• 提供研究基础： 预后研究为后续的临床试验提供了有价值的数据，为疾病的进一步研究打下基础，例如：后续研究可以重点关注某些指标/特征。

1.4 风险评分是什么？
风险评分是一种量化评估患者疾病风险或治疗结果的方法。它通过综合考虑患者的多个临床特征和生物标志物，计算出一个数值（即风险评分），用于预测患者的预后（如生存率、复发率等）或治疗反应。

1.5 风险评分和预后模型的关系是什么？
预后模型往往通过计算风险评分来实现对患者预后的预测，风险评分是预后模型中的核心组成部分，它反映了模型对患者临床结局的量化评估。 简单理解为：风险评分就是预后模型的量化指标，后续也是通过风险评分来对预后模型进行验证

1.6 多因素cox模型是什么？
多因素Cox模型，即Cox比例风险模型（Cox Proportional Hazards Model），是一种可以同时分析多个特征/变量对患者生存状态影响的统计模型。

1.7 多因素cox分析和之前的单因素cox分析的区别？
区别就在于： 单因素cox分析是对每个特征/变量单独进行分析，来研究该特征/变量是否与患者生存有联系，局限性就在于只分析了一种特征/变量。然鹅不同特征/变量之间或多或少都是存在联系的，而多因素cox分析就是同时分析多个特征/变量对于生存时间的影响，来看不同特征/变量对于患者生存的影响是否独立。

举个栗子： 有8个基因通过了单因素cox筛选，被认为是预后基因 （可以显著影响患者生存或复发） ，此时将8个基因进行多因素cox分析，有3个基因多因素cox计算的p < 0.05，此时这3个基因被认为是独立预后基因——虽然多个基因之间存在相互干扰的情况，但是这3个基因仍然显著影响患者生存或复发 。

1.8 多因素cox分析筛选独立预后基因和多因素cox构建预后模型有什么区别？
目的不一样：前者是用来筛选基因的，后者是用输入的基因构建模型并计算风险评分

注：多因素cox计算风险评分这里有两种思路（在本章中会重点介绍第一种方式，在下一章中会介绍第二种方式）：
（1）第一种就是直接用输入的全部基因去构建多因素cox模型并计算风险评分。
（2）第二种是在第一种多因素cox模型的基础上加入逐步回归，剔除一些相对不重要的基因，剩下的基因构建的模型被认为是最优模型，拿这个模型去计算风险评分

综上所述： 预后模型就是一种预测工具，它可以通过多种特征/变量来预测未来某段时间患者的生存/复发概率，在预后模型中有一个量化指标叫做风险评分。在确认要输入的特征/变量之后，就可以去构建多因素cox预后模型并计算风险评分。

注意：
（1）与预后相关的分析输入的都是疾病样本，和对照样本无关！！
（2）输入的特征/变量必须是与预后相关的特征/变量，如何筛选预后特征，可见之前的章节：单因素cox筛选预后相关特征

---

---

2. 多因素cox预后模型（Rstudio）——代码实操

本项目以TCGA——肺腺癌为例展开分析
物种：人类（Homo sapiens）
R版本：4.2.2
R包：tidyverse，survival

废话不多说，代码如下：

2. 1 数据处理

设置工作空间：

rm(list = ls()) # 删除工作空间中所有的对象
setwd('/XX/XX/XX') # 设置工作路径
if(!dir.exists('./15_risk_model_multivariate_cox')){
  dir.create('./15_risk_model_multivariate_cox')
} 
setwd('./15_risk_model_multivariate_cox/') 

加载包：

library(tidyverse)
library(survival)

导入要分析的表达矩阵train_data ，并对train_data 的列名进行处理（这是因为在读入的时候系统会默认把样本id中的“-”替换成“.”，所以要给替换回去）

train_data <- read.csv("./data_fpkm.csv", row.names = 1, check.names = F)  # 行名为全部基因名，每列为样本名
colnames(train_data) <- gsub('.', '-', colnames(train_data), fixed = T)

train_data 如下图所示，行为基因名（symbol），列为样本名

导入分组信息表group

group <- read.csv("./data_group.csv", row.names = 1) # 为每个样本的分组信息（tumor和control）
colnames(group) <- c('sample', 'group')

group 如下图所示，第一列sample为样本名，第二列为样本对应的分组 （分组为二分类变量：disease和control）

导入样本对应的生存信息表survival_data

survival_data <- read.csv('./data_survival.csv', row.names = 1, check.names = F) # 样本生存信息（生存时间，生存状态）

survival_data 如下图所示，第一列sample为样本名，第二列为样本对应的生存状态（0表示存活，1表示死亡） ，第三列为样本对应的生存时间（单位是天）

接下来从group，train_data 和 survival_data中筛选出疾病样本

注：涉及到复发/生存的分析点基本都是用的疾病样本，不要附加对照样本！！！！

group <- group[group$group == 'disease', ] ## 筛选疾病样本
train_data <- train_data[, group$sample] ## 从全部表达矩阵中同样筛选出疾病样本
survival_data <- survival_data[survival_data$sample %in% group$sample, ] ## 从生存信息表中提取出疾病样本的生存信息

导入要用于构建模型的基因hub_gene (10个基因，假设这10个基因均通过了前面的单因素cox筛选是预后基因)

hub_gene <- data.frame(symbol = gene <- c('VPS13D', 'MFF', 'ACSL1', 'VDAC1', 'PRELID1', 'BAK1',
                                          'CYCS', 'BCL2L10', 'MPV17L', 'PHB'))
colnames(hub_gene) <- "symbol"

hub_gene 如下图所示，只有一列：10个基因的基因名

从train_data中取出这10个基因对应的表达矩阵，并且与之前准备的生存信息表从survival_data进行合并

dat <- train_data[rownames(train_data) %in% hub_gene$symbol, ] %>%
  t() %>% 
  as.data.frame() # 整理后行为样本名，列为基因名
dat$sample <- rownames(dat)
dat <- merge(survival_data, dat, var = 'sample')
dat <- column_to_rownames(dat, var = "sample") %>% as.data.frame()

dat 如下图所示，行为样本名，第一列OS为生存状态，OS.time为生存时间，后面的列均为基因的表达量。

2. 2 构建多因素cox模型（用输入的全部10个基因）

通过as.formula函数将要进行分析的生存状态（OS）和生存时间（OS.time）以及基因转换为R的公式对象

## 准备多因素cox的公式
multivariable_cox_formula <- as.formula(paste0('Surv(OS.time, OS)~',
                              paste(hub_gene$symbol,
                                    sep = '',
                                    collapse = '+')))

multivariable_cox_formula如下图所示

根据准备的公式拟合Cox比例风险回归模型

## 多因素cox分析
multivariabl_cox <- coxph(multivariable_cox_formula,
                          data = dat)

在查看多因素cox结果multivariabl_cox前需要进行ph假定检验

## ph假设检验
multivariabl_cox_ph <- cox.zph(multivariabl_cox)
multivariabl_cox_ph_table <- multivariabl_cox_ph$table[-nrow(multivariabl_cox_ph$table),] %>% as.data.frame()

PH假定是Cox比例风险回归模型的一个核心假设，它要求不同个体之间的风险比率（hazard ratio）在整个随访期间内保持恒定，即不受时间的影响。换句话说，不同个体的风险在整个时间范围内以相同的比率增加或减少

举个栗子：假设正在研究两种不同治疗方法（治疗A和治疗B）对癌症患者生存时间的影响。招募了100名癌症患者，其中50名接受治疗A，50名接受治疗B，目标是确定哪种治疗方法能让患者活得更长。此时要控制不论是在研究开始时、中途还是结束时，治疗A组相对于治疗B组的死亡风险都应该保持不变，否则无法说明是治疗效果差诱导的死亡还是由于其他因素干扰导致的死亡，这时模型的估计结果可能会受到偏差，导致无法准确判断哪种治疗方法更有效。

multivariabl_cox_ph_table 如下图所示

• chisq： 这一列显示的是卡方（Chi-squared）统计量的值。在进行PH假定检验时，这个统计量用于量化协变量效应是否随时间变化。具体来说，它是通过比较不同时间点的风险比率来计算的。如果chisq值很小（并且对应的p值也很大），那么可以认为PH假定成立，即风险比率在整个随访期间内保持恒定
• df： 这一列显示的是自由度（Degrees of Freedom）。在PH假定检验的上下文中，这个值通常等于1，因为在检验的是单个协变量是否满足PH假定。
• p： 这一列显示的是与chisq统计量相关联的p值。如果p值大于0.05，则没有足够的证据拒绝PH假定的零假设，即认为PH假定成立——风险比率在整个随访期间内保持恒定，即不受时间的影响。

从multivariabl_cox_ph_table 中可以看到第三列对应的p值大于0.05，说明所有基因均满足ph假定检验，可以往下查看多因素cox的结果，否则要剔除p < 0.05的基因重新构建模型。

接下来就可以提取多因素cox分析结果的详细信息了

## 提取多因素cox分析结果的详细信息
multivariabl_cox_res_all <- summary(multivariabl_cox)

multivariable_cox_res <- data.frame(gene = rownames(data.frame(multivariabl_cox_res_all$coefficients)),
                                    coef = data.frame(multivariabl_cox_res_all$coefficients)[, 1],
                                    p_value = signif(as.matrix(multivariabl_cox_res_all$coefficients)[,5],2),
                                    HR = signif(as.matrix(multivariabl_cox_res_all$coefficients)[,2],2), 
                                    HR.95L = signif(multivariabl_cox_res_all$conf.int[,3],2), 
                                    HR.95H = signif(multivariabl_cox_res_all$conf.int[,4],2))

write.csv(multivariable_cox_res, './multivariable_cox_res.csv')

multivariable_cox_res 如下图所示:

• gene——表示基因名称
• coef——为每个基因对应的系数（这个系数就是后面参与风险评分计算的）
• p_value——表示多因素cox计算的p值，小于0.05被认为该基因是独立预后基因
• HR——表示风险比值比（如果HR大于1表示该基因是个危险因素，如果HR小于1表示该基因是个安全因素
• HR.95L和HR.95H——表示HR波动的95上下置信区间，可以理解成波动范围
  

2. 3 计算风险评分

接下来就是计算该模型对应的风险评分，来量化评估该模型对患者未来健康状况或疾病进展的风险。

使用predict函数来根据之前构建的多因素cox模型来预测风险评分

riskScore <- predict(multivariabl_cox, newdata = dat, type = "lp") %>% as.data.frame() # 行为每个样本名，第一列为其对应的风险评分
colnames(riskScore) <- 'riskscore'

riskScore 如下图所示，行名为样本名，第一列riskscore就是对应样本的风险评分——就是给每个样本分配一个具体的分数，以反映其未来发生不良事件（如疾病复发、死亡等）的风险大小。

将风险评分和生存信息，基因表达量等整理到一张表格中，后续验证会用到

riskScore$sample <- rownames(riskScore)

cox_data <- dplyr::select(dat, OS, OS.time, hub_gene$symbol)
cox_data$sample <- rownames(cox_data)

risk <- merge(cox_data, riskScore, by = 'sample')
risk$risk <- ifelse(risk$riskscore > median(risk$riskscore), 1, 0)

write.csv(risk, './risk.csv')

risk 如下图所示，第一列sample为样本名，第二列OS就是样本对应的生存状态，OS.time就是样本对应的生存时间，中间的列就是那10个预后基因对应的表达量，倒数第二列riskscore就是前一步计算出来的风险评分，最后一列risk就是对应的风险分组，1表示高风险组，0表示低风险组，区分度是风险评分的中位值。

计算出来的风险评分，就是模型对患者未来健康状况或疾病进展风险的量化评估，后续将针对计算出来的风险评分去研究如何对预后模型进行验证，在本章中不做介绍，会在全部预后模型介绍完之后统一进行说明。

---

---

结语：

以上就是构建预后模型之多因素cox模型的所有过程，如果有什么需要补充或不懂的地方，大家可以私聊我或者在下方评论。

如果觉得本教程对你有所帮助，希望广大学习者能够花点自己的小钱支持一下（点赞旁的打赏按钮）作者创作（可以的话一杯蜜雪奶茶即可），感谢大家的支持~~~~~~ ^_^ !!!

祝大家能够开心学习，轻松学习，在学习的路上少一些坎坷~~~

---

• 目录部分跳转链接：零基础入门生信数据分析——导读