零基础入门转录组数据分析——单因素cox筛选预后相关特征

零基础入门转录组数据分析——单因素cox筛选预后相关特征

---

---

1. 单因素cox基础知识

1.1 单因素cox是什么？
单因素Cox分析（Cox univariate analysis）是一种生存分析方法，用于评估某个单独的因素对患者生存（OS）或复发风险（Recurrence）的影响。

1.2 做单因素cox的目的是什么？
通过单因素Cox分析，可以确定某个因素是否与患者生存或复发显著相关。在转录组分析中，由于分析的都是基因的表达量，因此最后能够筛选出与患者生存显著相关的基因。

1.3 单因素cox的优点？

• 专注于单一因素： 单因素Cox分析能够精确地评估单个因素对生存时间或复发风险的影响，有助于深入理解该因素在疾病进程中的作用。
• 提供风险比： 通过估计风险比（hazard ratio, HR），单因素Cox分析能够量化因素对生存或复发风险的影响，为临床决策提供依据。
• 适用于多种数据类型： 无论是分类变量（如性别、治疗方法）还是连续变量（如年龄、基因表达水平），单因素Cox分析都能进行处理。

1.4 单因素cox的缺点？

• 无法考虑多个因素： 单因素Cox分析的主要局限性在于它只能评估单个因素的影响，而无法同时考虑多个因素的交互作用，这可能导致在分析复杂疾病时忽略掉一些重要的信息。
• 对样本量要求较高： 对于小样本量或事件数较少的数据集，单因素Cox分析可能产生不稳定的估计结果，因此，在进行分析前需要确保样本量足够大（关于这个数目选择没有一个硬性要求，但只能说是越多越好，如果非要说一个最低样本数目，个人建议不要低于30个疾病样本）。

1.5 什么是预后？
预后（Prognosis）在医学领域是一个重要的概念，它指的是基于疾病的发展规律、治疗情况以及患者的个体特征等因素，对疾病可能的发展过程、转归和结局进行的预测和评估。

1.6 预后的评估包括哪些方面？

预后评估通常涉及多个方面，包括但不限于：

• 病情评估（对疾病的严重程度、分期、病理类型等进行评估，了解疾病的当前状态）
• 治疗反应预测（基于患者的疾病特征和既往治疗经验，预测患者对特定治疗的反应和疗效）
• 生存期估计（根据疾病的特点和患者的个体情况，估计患者的生存时间和生存质量）
• 并发症和复发风险评估（评估患者发生并发症和疾病复发的风险，以便及时采取预防措施）

注：涉及到转录组分析中常用的就是生存数据和复发风险

1.7 为什么章节名称是“单因素cox筛选预后相关特征”？
前述单因素cox会评估某个单独的因素对患者生存（OS）或复发风险（Recurrence）的影响，而刚好预后评估的范畴包括生存期估计以及复发风险评估，所以通过单因素cox筛选的特征/变量通常被叫做是预后特征/变量。

举个栗子： 有8个基因，想看一下哪些基因与患者的生存/复发显著相关，就可以通过单因素cox进行筛选，通常p值选择0.05，最后结果显示有3个基因对应的单因素cox计算结果p值小于0.05，那么就认为这3个基因与该基因患者生存/复发有关联，在转录组中就可以叫做预后基因。

综上所述： 单因素cox就是筛选与患者生存相关的特征，并且可以提供每个特征对应的风险比，筛选后的特征或变量就被称为是预后特征/变量

注意：做单因素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('./13_univariate_cox')){
  dir.create('./13_univariate_cox')
} 
setwd('./13_univariate_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和normal）
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 (6个基因)

hub_gene <- data.frame(symbol = gene <- c('ADAMTS2', 'ADAMTS4', 'AGRN', 'COL5A1', 'CTSB', 'FMOD'))
colnames(hub_gene) <- "symbol"

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

从train_data中取出这6个基因对应的表达矩阵，并且与之前准备的生存信息表从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分析

注：这里以一个基因为例展示单因素cox分析，循环进行单因素cox分析的脚本到时候会放到配套资源中

选择一个基因AGRN进行单因素cox

• coxph： 这是R语言中survival包的一个函数，用于拟合Cox比例风险回归模型
• Surv(time = OS.time, event = OS)： 这是survival包中的一个函数，用于创建一个生存对象。time参数表示生存时间，这里是OS.time；event参数表示事件是否发生，这里是OS。如果OS是1，表示事件发生（例如，患者死亡）；如果OS是0，表示事件未发生（例如，患者存活）
• ~ AGRN： 这表示模型的公式部分，AGRN表示从数据框dat中选择名为AGRN的列作为模型中的变量
• data = dat： 指定数据集dat为模型的数据来源

## 拟合Cox比例风险回归模型
unicox <- coxph(Surv(time = OS.time, event = OS) ~ AGRN, data = dat)

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

cox_zph <- cox.zph(unicox) ## ph假定检验
cox_table <- data.frame(cox_zph$table) ## 将ph假定检验结果转成数据框形式

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

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

cox_table 如下图所示

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

从cox_table 中可以看到第三列对应的p值大于0.05，说明满足ph假定检验，可以往下查看单因素cox的结果，否则要剔除该基因。

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

## 提取单因素cox分析结果的详细信息
unisum <- summary(unicox)

最后一步就是进一步提取常规单因素cox分析会用到的一些信息，并将其储存起来。

## 提取单因素cox常规分析要用到的信息
univariate_cox_result <- data.frame(gene = 'AGRN',
                   HR = unisum$coefficients[, 2], # coef是公式中的回归系数b（有时候也叫beta值）, exp(coef)是cox模型中的风险比(HR), z代表wald统计量, 是coef除以其标准误se(coef)。
                   pvalue = unisum$coefficients[, 5],
                   L95CI = unisum$conf.int[, 3], # lower .95 upper .95则是exp(coef)的95%置信区间，可信区间越窄，可信度越高，你的实验越精确，越是真理。
                   H95CI = unisum$conf.int[, 4]
)

univariate_cox_result 如下图所示：

• gene——表示基因名称
• HR——表示风险比值比（如果HR大于1表示该基因是个危险因素，如果HR小于1表示该基因是个安全因素
• pvalue——表示单因素cox计算的p值，小于0.05被认为该基因与患者生存有联系
• L95CI和H95CI——表示HR波动的95上下置信区间，可以理解成波动范围

最后将单因素cox结果保存出去

## 导出结果
write.csv(univariate_cox_result, './uncox_res.csv')

注：以上过程只分析了一个基因的，单因素cox分析就是挨个对每个基因进行分析，因此往后就是重复的工作，可以用循环实现，再此不做展示。

---

---

结语：

以上就是单因素cox筛选预后相关特征的所有过程，如果有什么需要补充或不懂的地方，大家可以私聊我或者在下方评论。

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

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

---

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