零基础入门转录组数据分析——预后模型之lasso模型

1. 预后模型和lasso模型基础知识

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

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

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

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

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

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

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

1.6 lasso是什么？
Lasso（Least Absolute Shrinkage and Selection Operator）回归，由Robert Tibshirani在1996年提出，是一种通过构造惩罚函数来实现变量选择和正则化的回归分析方法。它通过向损失函数中添加L1正则化项（即系数绝对值的和），来压缩一些回归系数，甚至使一些回归系数变为零，从而实现特征选择和降维。

1.7 lasso预后模型和之前章节介绍的lasso部分的区别是什么？
Lasso预后模型是一种结合Lasso回归方法的临床预后评估工具，它利用Lasso回归进行特征选择和系数压缩，从而构建出具有预测患者预后能力的模型。

举个栗子： 有8个基因，先通过lasso算法的特征选择思路从中筛选出4个更加重要的基因（筛选基因的目标是构建一个既简洁又准确的预测模型），之后基于这4个基因构建lasso预后模型，并计算风险评分。

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

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

2. lasso预后模型（Rstudio）——代码实操

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

废话不多说，代码如下：

2. 1 数据处理

设置工作空间：

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

加载包：

library(tidyverse)
library(survival)
library(glmnet)

导入要分析的表达矩阵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 (12个基因)

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

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

从train_data中取出这12个基因对应的表达矩阵，并且与之前准备的生存信息表从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 构建lasso预后模型

先准备要输入的数据x_all 和y_all
x_all 就是dat中去除了OS和OS.time这两列剩余的基因表达矩阵
y_all 就是dat中仅选择OS和OS.time这两列

## 准备输入数据
x_all <- subset(dat, select = -c(OS, OS.time))
y_all <- subset(dat, select = c(OS, OS.time))

接下来构建lasso预后模型：使用cv.glmnet函数通过调整正则化参数（lambda）来寻找最优的lasso模型，以平衡模型的复杂度和拟合度，cv.glmnet函数常用参数介绍如下：

x参数——是要输入的基因表达矩阵（也称为特征或自变量），as.matrix函数将其转换为矩阵格式，因为glmnet函数需要输入数据为矩阵形式
y 参数——这部分使用survival包中的Surv函数创建生存对象。OS.time是生存时间，OS是一个指示变量，表示每个观测在观测结束时是否发生了事件（比如：死亡）。Surv函数将这两个变量结合，生成一个包含生存时间和事件状态的对象，供Cox模型使用。
family参数——这个参数指定了模型的类型，"cox"表示指定了模型家族为Cox比例风险模型
nfolds 参数——这个参数指定了交叉验证的折数（folds）。在这个例子中，数据集被分为15个大致相等的部分，模型在其中的14部分上进行训练，并在剩下的1部分上进行测试。这个过程重复15次，每次选择不同的部分作为测试集。通过这种方式，可以对整个数据集的不同部分评估模型的性能，并得到一个平均的性能指标，从而帮助选择最佳的lambda值

（cv.glmnet函数中比较关注的参数就是上述的这些，当然还有其他参数，如果想深入了解可自行查看官方说明文档）

set.seed(754) ## 设置种子
cvfit <- cv.glmnet(x = as.matrix(x_all), y = Surv(y_all$OS.time, y_all$OS), nfolds=15,family = "cox")

注：在构建模型的时候切记要设置种子（设置随机种子是为了确保结果的可重复性。由于交叉验证涉及随机分割数据，因此设置种子可以确保每次运行代码时，数据的分割方式都是相同的，从而得到相同的模型结果）

接下来从构建的最优模型中提取每个基因对应的系数（简单说就是提取基因）

## 提取指定lambda时特征的系数
coef.min <- coef(cvfit, s = "lambda.min") ## lambda.min & lambda.1se 取一个

注：lambda.min 或 lambda.1se 取一个就行，看个人需求，我这边选择lambda.min旨在考虑最佳预测性能，不考虑稳健性

lambda.min——是交叉验证过程中使得平均偏差（或其他指定的损失函数）最小的lambda值。换句话说，它是根据交叉验证结果，直接选择出的最优lambda值，因为它在平均意义上给出了最好的预测性能
lambda.1se——lambda.1se是lambda.min之后，使得平均偏差增加不超过一倍标准误差的最大的lambda值。这个选择标准旨在在模型的预测性能和复杂度之间找到一个折中。通过选择一个稍微更大的lambda值（即更强的正则化），lambda.1se可以帮助避免过拟合，在接受一些预测性能的损失的同时来换取更好的稳健性

coef.min中就对应着每个基因的系数，如下图所示，这些系数在后面就参与了风险评分的计算
在这里插入图片描述

2. 3 提取Lasso预后基因

提取出最优Lasso模型中筛选出来的那些基因及其对应的系数并保存（筛选出来的基因被认为相对更加重要一些）

# 找出那些回归系数没有被惩罚为0的
active.min <- which(coef.min@i != 0)# 提取基因名称
lasso_geneids <- coef.min@Dimnames[[1]][coef.min@i+1]df.coef <- cbind(gene = rownames(coef.min), coefficient = coef.min[,1]) %>% as.data.frame()
df.coef <- subset(df.coef, coefficient != 0) %>% as.data.framewrite.csv(df.coef, "Lasso_Coefficients.csv")
write.csv(lasso_geneids, './lasso_genes.csv')

df.coef如下图所示，这里面就是筛选出来的8个预后基因，可以看到相较于一开始输入的12个基因，剔除了4个与预后不相干的基因。
在这里插入图片描述

2. 4 计算风险评分

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

使用glmnet包的predict函数来根据之前通过交叉验证得到的最佳模型（cvfit）来预测风险评分

riskScore <- predict(cvfit, newx = as.matrix(x_all), s = cvfit$lambda.min) %>% as.data.frame() # 行为每个样本名，第一列为其对应的风险评分
colnames(riskScore) <- 'riskscore'

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

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

riskScore$sample <- rownames(riskScore)lasso_data <- dplyr::select(dat, OS, OS.time, lasso_geneids)
lasso_data$sample <- rownames(lasso_data)risk <- merge(lasso_data, riskScore, by = 'sample')
risk$risk <- ifelse(risk$riskscore > median(risk$riskscore), 1, 0)write.csv(risk, './risk.csv')

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