零基础入门转录组数据分析——预后模型之随机生存森林模型

呆猪儿

于 2024-08-29 15:42:54 发布

阅读量392

点赞数 5

分类专栏：生信之转录组——下游分析文章标签：数据分析人工智能算法机器学习数据挖掘 r语言大数据

本文链接：https://blog.csdn.net/weixin_49878699/article/details/141064420

版权

生信之转录组——下游分析专栏收录该内容

20 篇文章 104 订阅

订阅专栏

零基础入门转录组数据分析——预后模型之随机生存森林模型

1. 预后模型和随机生存森林模型基础知识

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

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

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

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

指导临床决策： 预后模型可以帮助医生预测疾病的可能进展，从而为患者提供更加个性化的治疗建议。例如，在癌症治疗中，预后模型可以预测患者的复发风险和生存时间，帮助医生制定更加合理的治疗方案。
疾病风险评估： 预后模型能够识别哪些患者更有可能出现不良结局，使医生可以进行早期干预，降低患者的风险。
提供研究基础： 预后研究为后续的临床试验提供了有价值的数据，为疾病的进一步研究打下基础，例如：后续研究可以重点关注某些指标/特征。

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

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

1.6 随机生存森林模型是什么？
随机生存森林模型（Random Survival Forest，简称RSF）是一种结合了随机森林和生存分析思想的机器学习模型。具体来说，随机生存森林模型通过训练大量生存树（Survival Trees），并以表决的形式，从个体树中加权选举出最终的预测结果。这种模型在能够处理生存数据，预测个体的生存时间或生存概率，并评估不同因素对生存时间的影响。

1.7 随机森林的基础
随机森林是一种集成学习方法，它通过构建多棵决策树并将它们的预测结果进行汇总，来提高整体预测的准确性和稳定性。每棵决策树都是基于随机选取的样本和特征子集进行训练的，这样做的好处是可以减少模型对特定数据的过拟合，同时增加模型的泛化能力。

1.8 生存分析的概念
生存分析是研究生存时间和生存状态之间关系的统计方法。在生物医学领域，生存时间可能指的是患者的存活时间，而生存状态则是指患者是否存活或是否发生了某个特定事件（如疾病复发、死亡等）。生存分析的目标是预测个体的生存时间或生存概率，并评估不同因素对生存时间的影响。

1.9 随机生存森林模型的优势是什么？

高准确性： 能够处理高维数据和大规模数据集，并具有较高的预测准确性。
鲁棒性： 对于数据中的噪声、异常值和缺失值具有较好的稳定性，能够有效地处理复杂的现实问题。
可解释性： 通过特征的重要性排序，可以提供特征选择和变量重要性分析的结果，给每个变量的选择提供理论依据

综上所述： 随机生存森林模型就是采用随机森林的算法原理来评估不同因素对生存时间的影响。

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

2. 随机生存森林模型（Rstudio）——代码实操

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

废话不多说，代码如下：

2. 1 数据处理

设置工作空间：

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

加载包：

library(tidyverse)
library(randomForestSRC)

导入要分析的表达矩阵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 (6个基因，假设这6个基因均通过了前面的单因素cox筛选是预后基因)

hub_gene <- data.frame(symbol = gene <- c('VPS13D', 'MFF', 'ACSL1', 'VDAC1', 'PRELID1', 'BAK1'))
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 构建随机生存森林模型（用输入的全部6个基因）

设定种子数后直接通过rfsrc函数构建随机生存森林模型

set.seed(123) —— 设定随机种子数，这样做可以确保你的模型结果具有可重复性
Surv(OS.time, OS) ~ . —— 这里，Surv(OS.time, OS)是一个生存对象，用于指定生存时间（OS.time）和事件指示符（OS，通常1表示死亡，0表示存活）。~ .表示模型将使用数据框dat中的所有其他列作为预测变量（特征）。
data = dat —— 表明用到的是上述处理好的dat数据框
ntree = 100 —— 这个参数指定了随机森林中要构建的树的数量。增加树的数量可以提高模型的预测准确性，但也会增加计算时间和模型的复杂度，要找到一个平衡点。
mtry = 4 —— 在随机森林中，mtry指定了每次树分裂时随机选择的候选变量的数量。较小的mtry值可以增加模型的多样性，但可能会牺牲一些预测准确性；较大的mtry值则可能使模型过于复杂，导致过拟合，同样是要找到一个平衡点。

set.seed(123)  
RSF_model <- rfsrc(Surv(OS.time, OS) ~ ., 
                     data = dat,
                     ntree = 100, 
                     mtry = 4)

接下来从构建的模型中提取每个基因的重要性，以此来观测哪些基因更为重要

# 提取特征重要性
Importance <- vimp(RSF_model)$importance

# 转换为数据框
Importance_df <- data.frame(variable = names(importance), importance = importance)

Importance_df如下图所示，第一列为基因名，第二列为每个基因对应的重要性。
在这里插入图片描述

2. 3 计算风险评分

接下来就用这6个基因构成的随机生存森林模型计算风险评分，这个风险评分就是模型对患者未来健康状况或疾病进展风险的量化评估。

使用predict函数来根据之前构建的随机生存森林模型来预测风险评分，并根据风险评分的中位值区分高低评分组

# 获取训练集的风险评分
dat$riskScore <- predict(RSF_model, newdata = dat)$predicted
dat$risk <- ifelse(dat$riskScore > median(dat$riskScore), 1, 0)

risk <- dat
write.csv(risk, './risk.csv')
write.csv(Importance_df, './gene_importance.csv')

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