零基础入门转录组数据分析——机器学习算法之lasso

1. Lasso基础知识

1.1 lasso是什么？
Lasso回归算法（Least Absolute Shrinkage and Selection Operator，最小绝对收缩和选择算子回归）是一种用于特征选择和参数估计的线性模型构建方法。

1.2 lasso是用来干什么的？
主要目的是为了解决传统线性回归在处理高维数据时遇到的问题，如变量选择困难、模型过拟合等。
通俗来说：就是筛选关键特征

1.3 lasso基本原理是什么？
Lasso回归通过在损失函数中加入一个L1正则化项（即变量系数的绝对值之和）来压缩模型，使某些系数缩减至零，从而实现了对变量的选择。

1.4 lasso优点是什么？
（1）自动进行变量选择，简化了模型，提高了模型的解释性。
（2）能够有效处理参数的多重共线性问题，可以控制模型的复杂度，避免过拟合。

1.5 lasso的本质是什么？
筛选出一些关键特征，这些关键特征相对于其他特征来说，区分样本的能力更加精确。

举个栗子： 基因A在疾病组中普遍表达为9-10， 基因A在对照组中普遍表达为5-7，这样明显就是 基因A在疾病组中高表达， 基因A表达量对于疾病和对照区分度非常高； 基因B在疾病组中普遍表达9-10， 基因B在对照组中普遍表达8-10，这样 基因B就无法非常好的区分疾病和对照，那么 基因B表达量对于疾病和对照区分度比较低。 基因A 相较于 基因B来说就是关键特征。

综上所述： lasso就是一种用来筛选 关键特征 的方法，这个关键特征可以是临床指标，也可以是重要基因等，并且在关键特征选择的时候避免了人为的选择，提高了可解释性。

2. Lasso（Rstudio）——代码实操

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

废话不多说，代码如下：

2. 1 数据处理

设置工作空间：

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

加载包：

library(tidyverse)
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和normal）
colnames(group) <- c('sample', 'group')

group 如下图所示，第一列sample为样本名，第二列为样本对应的分组 （分组为二分类变量：disease和control）
在这里插入图片描述
导入要筛选的基因hub_gene (8个基因)

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

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

从全部的基因表达矩阵中取出这8个基因对应的表达矩阵，并且与之前准备的分组信息表进行合并

dat <- train_data[rownames(train_data) %in% hub_gene$symbol, ] %>%
  t() %>%
  as.data.frame() # 整理后行为样本名，列为基因名
dat$sample <- rownames(dat)
dat <- merge(dat, group, var = "sample") # 筛选后的表达矩阵与分组信息表合并
dat <- column_to_rownames(dat, var = "sample") %>% as.data.frame() # 列名转为行名

table(dat$group)
dat$group <- factor(dat$group, levels = c('disease', 'control'))

dat 如下图所示，行为基因名，前8列为基因对应的表达矩阵，第9列为合并的分组信息表
在这里插入图片描述

2. 2 构建Lasso模型

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

x参数——是要输入的基因表达矩阵（也称为特征或自变量）
y 参数——这是要区分的目标变量。这里指向的是dat$group，根据前面得知这里是二分类变量分组——disease和control
family参数——这个参数指定了模型的类型，"binomial"表示二项分布，适用于二分类问题
type.measure参数——这个参数指定了交叉验证时用于评估模型性能的度量标准。“default"意味着对于"binomial”，默认使用偏差（deviance）作为度量标准，偏差是模型拟合优度的一种度量，较低的偏差表示更好的拟合。
nfolds 参数——这个参数指定了交叉验证的折数（folds）。在这个例子中，数据集被分为10个大致相等的部分，模型在其中的9部分上进行训练，并在剩下的1部分上进行测试。这个过程重复10次，每次选择不同的部分作为测试集。通过这种方式，可以对整个数据集的不同部分评估模型的性能，并得到一个平均的性能指标，从而帮助选择最佳的lambda值

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

set.seed(381) # 设置种子
res.lasso <- cv.glmnet(x = as.matrix(dat[-ncol(dat)]), 
                       y = dat$group, 
                       family = "binomial", 
                       type.measure = "default",
                       nfolds = 10)

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

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

coef.min <- coef(res.lasso, 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]
lasso_geneids <- lasso_geneids[-1] %>% as.data.frame()
colnames(lasso_geneids) <- 'symbol'

write.csv(lasso_geneids, file = 'lasso_gene.csv')

lasso_geneids如下图所示，这里面就是筛选出来的关键基因，可以看到相较于一开始输入的8个基因，剔除了一些相对来说不重要的基因（CTSB）
在这里插入图片描述
注：从上面的分析也可以看到最后筛选出来的基因是根据算法自动得出来的，这样就避免了人为选择的情况，提高了模型的可解释性

2. 4 Lasso结果简单可视化

接下来一步就是要对lasso结果进行简单可视化，毕竟文章里是要放图的!!!

# lasso结果简单可视化
par(mfrow = c(1, 2))
plot(res.lasso$glmnet.fit, 
     xvar = 'lambda', 
     label = TRUE, 
     las = 1)
abline(v=log(c(res.lasso$lambda.min, res.lasso$lambda.1se)), lty="dashed")
plot(res.lasso, 
     las =1)