零基础入门转录组数据分析—

本文主要是介绍零基础入门转录组数据分析——预后模型的验证，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

零基础入门转录组数据分析——预后模型的验证

零基础入门转录组数据分析——预后模型的验证
- 1. 预后模型的基础知识
- 2. 预后模型的验证（Rstudio）——代码实操
- - 2. 1 数据处理
  - 2. 2 构建多因素cox模型（用输入的全部5个基因）
  - 2. 3 计算风险评分
  - 2. 4 生存分析（验证一）
  - 2. 5 ROC分析（验证二）

1. 预后模型的基础知识

关于预后模型的其他基础知识在此不做过多介绍，前面章节已经有过多次介绍，可从目录进行跳转零基础入门转录组数据分析——目录

1.1 为什么要对预后模型进行验证？

评估预测性能： 预后模型旨在预测患者未来某一时刻发生某事件的概率（如复发、死亡、伤残等）。通过验证，可以评估模型在不同数据集中的预测准确性，确保其在不同情境下都能提供可靠的预测结果。
减少偏差： 模型在开发过程中可能存在过度拟合或欠拟合的情况，导致在特定数据集上表现优异，但在其他数据集中表现不佳，验证过程有助于发现这些问题，并通过调整模型来减少偏差。
提高数据分析可信度： 跟实验重复一样，用别的数据集进行验证，如果能得到相似的结果，那么就说明构建出来的模型具有普适性。

1.2 预后模型的验证通常怎么选择数据集？
通常选择两个数据集，一个作为训练集，一个作为外部验证集，在训练集中筛选基因并构建模型，之后引入一个外部验证集对构建出来的模型性能进行验证。

1.3 预后模型的验证通常包括哪些分析？

生存分析： 要保证训练集和验证集的高低风险组间生存均存在显著差异
ROC分析： 通常选择年份是3个连续的年份，例如：1，2，3年，要求训练集和验证集123年的ROC分析中AUC值都要大于0.6

综上所述： 为了验证构建出来的模型是否在别的数据集中存在普适性，因此要通过生存分析以及ROC分析在另一个数据集（又称为：外部验证集）中验证模型的预测能力。

注意：预后模型的验证要用到前面章节中不同模型最后计算出来的风险评分，如果对于构建模型计算风险评分不了解，可以先看看之前的章节。

本教程将以多因素cox模型展示如何验证模型的普适性，关于其他模型的验证，大家可以参考本贴后自行发挥

2. 预后模型的验证（Rstudio）——代码实操

本项目以TCGA——肺腺癌（LUAD）作为训练集，GSE50081作为外部验证集展开分析
物种：人类（Homo sapiens）
R版本：4.2.2
R包：tidyverse, survival, survminer, timeROC

废话不多说，代码如下：

2. 1 数据处理

设置工作空间：

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

加载包：

library(tidyverse)
library(survival)
library(survminer)
library(timeROC)

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

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

TrainRawData如下图所示，行为基因名（symbol），列为样本名
在这里插入图片描述
导入分组信息表TrainGroup

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

TrainGroup 如下图所示，第一列sample为样本名，第二列为样本对应的分组 （分组为二分类变量：disease和control）
在这里插入图片描述
导入样本对应的生存信息表TrainSurvivalData

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

TrainSurvivalData 如下图所示，第一列sample为样本名，第二列为样本对应的生存状态（0表示存活，1表示死亡），第三列为样本对应的生存时间（单位是天）
在这里插入图片描述

接下来从TrainGroup，TrainRawData 和 TrainSurvivalData 中筛选出疾病样本

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

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

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

HubGene <- data.frame(symbol = gene <- c('VDAC1', 'CYCS', 'ACSL1', 'GLS2', 'MPV17L'))
colnames(HubGene) <- "symbol"

HubGene 如下图所示，只有一列：5个基因的基因名
在这里插入图片描述

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

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

TrainData 如下图所示，行为样本名，第一列OS为生存状态，第二列OS.time为生存时间，后面的列均为基因的表达量。
在这里插入图片描述

准备好训练集数据后，接下来用相同的方式整理外部验证集（GSE50081）的数据

## 验证集
TestRawData <- read.csv("./GSE50081/dat.GSE50081.csv", row.names = 1, check.names = F)  # 行名为全部基因名，每列为样本名TestGroup <- read.csv("./GSE50081/group.GSE50081.csv") # 为每个样本的分组信息（tumor和normal）
colnames(TestGroup) <- c('sample', 'group')TestSurvivalData <- read.csv('./GSE50081/survival.GSE50081.csv', row.names = 1, check.names = F) # 样本生存信息（生存时间，生存状态）TestGroup <- TestGroup[TestGroup$group == 'LUAD', ] ## 筛选疾病样本
TestData <- TestRawData[, TestGroup$sample] ## 从全部表达矩阵中同样筛选出疾病样本
TestSurvivalData <- TestSurvivalData[TestSurvivalData$sample %in% TestGroup$sample, ] ## 从生存信息表中提取出疾病样本的生存信息TestData <- TestData[rownames(TestData) %in% HubGene$symbol, ] %>%t() %>% as.data.frame() # 整理后行为样本名，列为基因名
TestData$sample <- rownames(TestData)
TestData <- merge(TestSurvivalData, TestData, var = 'sample')
TestData <- column_to_rownames(TestData, var = "sample") %>% as.data.frame()

TestData 如下图所示，和训练集TrainData 类似。
在这里插入图片描述

检查训练集（TrainData）和验证集（TestData）列名中的基因是否一致

## 检测两个数据集基因是否相同
identical(colnames(TrainData), colnames(TestData))

注：这里必须要一致，是因为有的数据集中没有定量出同样的基因，出现不一致的情况，只能考虑换别的基因

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

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

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

MultivariableCoxFormula 如下图所示
在这里插入图片描述

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

## 设定种子数
set.seed(123)
## 多因素cox分析
MultivariableCoxModel <- coxph(MultivariableCoxFormula,data = TrainData)

注：在用多因素cox模型往后分析的时候，通常先需要进行PH假定检验，关于做PH假定检验的方式可以参考之前的章节，这里我们假定用于构建模型的5个基因都通过了PH假定检验，直接用于后续分许。

2. 3 计算风险评分

接下来就是用该模型计算不同样本对应的风险评分，来量化评估该模型对患者未来健康状况或疾病进展的风险。并根据风险评分的中位值将不同样本区分为高低风险组。

## 训练集
TrainData$RiskScore <- predict(MultivariableCoxModel, newdata = TrainData, type = "lp")
TrainData$Risk <- ifelse(TrainData$RiskScore > median(TrainData$RiskScore), 1, 0)
## 验证集
TestData$RiskScore <- predict(MultivariableCoxModel, newdata = TestData, type = "lp")
TestData$Risk <- ifelse(TestData$RiskScore > median(TestData$RiskScore), 1, 0)

TrainData如下图所示，可以看到最后多出来两列，一个是RiskScore，这一列就是用多因素cox模型计算出来的风险评分；另一个是Risk，这一列就是根据风险评分的中位值区分的高低风险组，1为高风险组，0为低风险组。
在这里插入图片描述

2. 4 生存分析（验证一）

为了研究区分出来的高低风险组间是否在生存上存在显著差异，用上一步在训练集和验证集划分的高低风险组进行生存分析

使用survfit函数计算生存曲线

Surv(OS.time, OS) ——调用Surv函数的表达式，用于创建生存对象，OS.time对应的就是生存时间，OS对应的生存状态
~ Risk —— 用于指定数据集中的分组，这里0表示低风险组，1表示高风险组

之后使用surv_pvalue函数来获取生存曲线之间的统计显著性，并合并结果

## 训练集生存分析
TrainKm <- survfit(Surv(OS.time, OS) ~ Risk, data =  TrainData)
TrainKmResults <- surv_pvalue(TrainKm, data = TrainData)
TrainKmResults$Type <- 'Train'
## 验证集生存分析
TestKm <- survfit(Surv(OS.time, OS) ~ Risk, data =  TestData)
TestKmResults <- surv_pvalue(TestKm, data = TestData)
TestKmResults$Type <- 'Test'
## 合并训练集和验证集生存分析的结果
FinallyKmResults <- rbind(TrainKmResults, TestKmResults)

FinallyKmResults 如下图所示，method列为Log-rank检验，该方法也称为对数秩检验，是一种非参数检验方法，用于比较两个或多个生存曲线（或称为失败时间分布）之间的差异是否具有统计学意义。通常第二列pval如果小于0.05，则认为两组间生存存在显著差异。
在这里插入图片描述

结果发现在训练集和验证集中计算的风险评分区分的高低风险组患者间生存均存在显著差异

2. 5 ROC分析（验证二）

ROC分析（Receiver Operating Characteristic），是一种将灵敏度和特异度结合起来综合评价诊断准确度或判别效果的方法。ROC分析的核心在于绘制ROC曲线，该曲线以虚惊概率（即1-特异度）为横轴，击中概率（即灵敏度）为纵轴，展示了在不同分类阈值下分类器的性能。简单来说就是：ROC曲线下的面积（AUC）是一个常用的分类器性能度量，AUC值越接近1，表示分类器性能越好。

注意：预后模型ROC曲线的AUC都一般比较低，通常文章认为3个时间节点的AUC只要大于0.6就认为该模型预测性能较好

使用timeROC函数计算1，2，3年的ROC曲线

T = TrainData$OS.time —— 代表每个观察对象的生存时间（以天为单位）。
delta = TrainData$OS —— delta是一个二进制（或类似二进制）的向量，用于指示在每个相应的时间点是否发生了感兴趣的事件（如死亡）。这里，1表示个体死亡，0表示在随访结束时个体仍然存活。
marker = TrainData$RiskScore ——是模型预测的风险评分，用于评估模型在预测生死事件是否发生方面的性能。
cause = 1 —— 参数用于指定感兴趣的事件类型，在这里，1表示只关注个体死亡的事件。
weighting = “marginal” —— 参数指定了如何考虑数据中的权重或样本大小，“marginal” 选项意味着在计算ROC曲线时，不考虑之前时间点的影响，而是基于每个时间点的边缘分布进行独立评估。
times = TimePoint —— 参数是一个向量，指定了计算ROC曲线时要考虑的时间点。在这里，TimePoint 是之前定义的向量，包含了三个时间点（以天为单位），分别是1年、2年和3年的时间点。
ROC = TRUE —— 参数是一个逻辑值，指示函数是否应该计算ROC曲线，在这里，TRUE 表示应该计算ROC曲线及其相关的性能指标（如AUC值）。
iid = TRUE —— 参数用于指定观察对象是否是独立同分布的。在生存分析中，通常假设观察对象不是完全独立的，因为它们之间可能存在时间依赖性或其他类型的依赖关系。将iid设置为FALSE更符合生存数据的实际情况，因为它承认了观察对象之间的潜在依赖性。

## 训练集ROC分析
TimePoint <- c(1, 2, 3)*365
TrainTimeRoc <- timeROC(T = TrainData$OS.time,delta = TrainData$OS,marker = TrainData$RiskScore,cause = 1,weighting = "marginal",times = TimePoint,ROC = TRUE,iid = FALSE
)
TrainTimeRocResults <- TrainTimeRoc$AUC %>% t() %>% as.data.frame()
TrainTimeRocResults$Type <- 'Train'
## 验证集ROC分析
TestTimeRoc <- timeROC(T = TestData$OS.time,delta = TestData$OS,marker = TestData$RiskScore,cause = 1,weighting = "marginal",times = TimePoint,ROC = TRUE,iid = FALSE
)
TestTimeRocResults <- TestTimeRoc$AUC %>% t() %>% as.data.frame()
TestTimeRocResults$Type <- 'Test'
## 合并训练集和验证集ROC分析的结果
FinallyTimeRocResults <- rbind(TrainTimeRocResults, TestTimeRocResults)