Day06.大讨论！Lasso回归模型！

对于生信实操，我们认为，在线工具和 R 语言各有千秋，两者结合或许才是最好的策略。这句话对于生信学习，或者数据挖掘都适用。纯生信操作简单，适合入门和初步数据分析；R 门槛稍高，个性化强，适合深度数据挖掘。为此，在纯生信教程基础上，我们推出在线工具与 R 结合的新版21天教程（重磅推出 | 新版21天生信实操教程！）。

生信数据挖掘的策略

新版21天教程将以原纯生信教程以及《癌生物学》、《R 数据科学》和《詹韦免疫学》等经典权威教材的核心内容为基础，并结合肿瘤和非肿瘤疾病最新发表的论文展开。在21天学习结束后，芒果团队会集中展开指定文献中所有图表的复现，让参与学习的小伙伴真正实现“学了就能会、会了就能用”。这将是新版21天教程的最大特色！尽管本文没有进行Lasso回归模型构建，但其应用广泛，值得关注。本次我们主要分享Lasso回归。

往期推荐

Day01.最新生信思路 Fig1复现！差异分析和生存分析！尝鲜版！

Day02.Fig2复现！热图、富集分析和互作分析！尝鲜版！

Day03.预测模型，从批量生存分析谈起！

Day04.预测模型，聊聊ROC曲线！尝鲜版！

Day05.单因素、多因素COX回归分析！尝鲜版！

在与果友交流过程中，Lasso回归经常是与Cox回归同时被问及的话题。比如①预后相关基因是指做单因素Cox回归有意义的基因，还要深入做Lasso或者多因素Cox吗？要！②Lasso预后建模后得到的风险评估基因热图的差异不显著该怎么处理呢？AUC明明很高，但具体到模型内的各个基因在高低风险中的差异表达无显著区别是怎么回事？差异分析是前提！③预测模型做完单因素后可以直接做Lasso回归吗，就不做多因素了，因为Lasso回归的基因数只有几个了？④单因素Cox➕Lasso筛选预后相关DEGs➕多因素模型风险评分。这种模式能否只做Lasso➕多因素呢？因为我单基因筛选的基因过少，直接做lasso模型效果好一点？群里讨论了，不过常规是都要做，或者重置差异分析的条件。

Tibshirani在1996年引入 Lasso (Least Absolute Shrinkage and Selection Operator)的概念，用于参数的选择和收缩。Lasso回归是以缩小变量集为核心思想的压缩估计方法，通过构造一个惩罚函数，将变量系数进行压缩，并使某些不重要的标量回归系数变为0，从而选择最重要的变量进入下游分析。Lasso 回归是一种线性回归模型，也称正则化或惩罚回归模型，使用L1正则化项来限制变量的系数，从而用于变量数量较多的大数据集，而传统的线性回归模型无法处理这类大数据。交叉验证用于确证Lasso结果中的取值是否最佳。

Lasso回归和岭回归都是常用的正则化技术，用来减少系数的大小，但使用的正则项不同。Lasso回归使用L1正则项，岭回归使用L2正则项。L1正则项的特点是，它可以完全移除某些系数，使得模型中的变量更加稀疏，即变量较少，这对于特征选择非常有用。而岭回归不具有变量的筛选功能，它只能让变量之间的差异变得更小。但L1正则项也有一些缺点，由于L1正则项对所有变量的惩罚都是绝对值，所以对于变量之间的大小没有影响。因此，L1正则项会使得变量之间存在较大的差异，导致模型的可解释性降低。

接下来，我们进行Lasso回归的实操。首先准备数据。

####----11_Lasso_Regression----####rm(list = ls())gc()
load('TCGA_STAD_data.RData')
library(dplyr)library(tidyverse)library(ggplot2)library(survival)
library(glmnet) # 正则化技术library(car) # 多重共线性检验library(corrplot) # 相关性绘图library(ggpubr) #基于ggplot2的可视化包ggpubrlibrary(patchwork) # 拼图library(ROCR) # 画ROC曲线library(caret) # 切割数据
Data <- data[,c(6,12:30917)]Data$Status = ifelse(Data$Status =='Dead','1','0')Data <- as.data.frame(lapply(Data,as.numeric))
### 数据检查min(Data$month) # 查看最小值，看是否有0Data2 <- filter(Data, Data$month > 0) # 删除数据框中time值小于0的样本，有0的话后续分析会报错dataL <- Data2  # 数据备份

探究数据的特征， 绘制直方图和正态分布概率密度函数图

hist(dataL$FBLN5, main = 'Distribution of gene expression levels',     xlab = 'FBLN5', ylab = 'gene expression', col = 'skyblue')  # 直方图curve(dnorm(x, mean = mean(dataL$FBLN5), sd = sd(dataL$FBLN5)), add = TRUE, col = 'red')  # 正态分布的概率密度函数# shapiro-wilk检验shapiro.test(dataL$FBLN5)  # 适用于小样本数据，N≤50，若p值<0.05，则数据不符合正态分布# Kolmogorov-Smirnov检验ks.test(dataL$FBLN5, pnorm, mean = mean(dataL$FBLN5), sd = sd(dataL$FBLN5)) # 适用于大样本数据，N>50，若p值<0.05，则数据不符合正态分布# 对所有列的数据进行Kolmogorov-Smirnov循环检验ks_result <- data.frame('Gene' = character(), 'p.value' = numeric())for (i in 3:length(colnames(dataL))){  p.value <- ks.test(dataL[,i], pnorm, mean = mean(dataL[,i]), sd = sd(dataL[,i]))$p.value  Gene <- colnames(dataL)[i]  tmp <- data.frame(Gene = Gene,                    p.value = p.value)  ks_result <- rbind(ks_result, tmp)}max(ks_result$p.value)  # 若p值<0.05，则数据不符合正态分布

绘制相关矩阵，检查数据的多重共线性问题

cor_matrix <- cor(dataL[,3:12], use = 'everything', method = c('spearman'))corrplot(cor_matrix, type = 'upper', method = 'color', tl.cex = 0.8, tl.col = 'black', order = 'AOE')

多重共线性检查之方差膨胀因子（VIF）

dataL2 <- dataL[,1:12]model <- coxph(Surv(month, Status) ~ ., data = dataL2)vif(model)plot(vif(model),     xlim = c(0,10), ylim = c(0,10),     xlab = 'variables', ylab = 'VIF',     cex = 3, pch = 1)# 当VIF大于5时，认为存在多重共线性，当VIF大于10时，认为存在严重的多重共线性

构建Lasso回归模型

time <- dataL[,'month']status <- dataL[,'Status']x <- as.matrix(dataL[,-c(1,2)]) #x为输入特征,是矩阵格式y <- as.matrix(dataL$Status)lasso <- glmnet(x = x, y = y, family = 'binomial', alpha = 1, # alpha = 1为LASSO回归，= 0为岭回归，0和1之间则为弹性网络 nlambda = 100) # nlambda表示正则化路径中的个数，这个参数就可以起到一个阈值的作用，决定有多少基因的系数可以留下来。默认值为100。print(lasso)plot(lasso, xvar = 'lambda', label = TRUE) # 系数分布图，由“log-lambda”与变量系数作图，展示根据lambda的变化情况每一个特征的系数变化，展示了Lasso回归筛选变量的动态过程plot(lasso, xvar = 'dev', label = TRUE) # 也可以对%dev绘图plot(lasso, xvar = 'norm', label = TRUE) # “L1范数”与变量系数作图

交叉验证，挑选合适的λ值，变量筛选

set.seed(1234)  # 设置种子数，保证每次得到的交叉验证图是一致的# 计算100个λ，画图，筛选表现最好的λ值lasso_cv <- cv.glmnet(x = x, y = y, family = 'binomial', alpha = 1, nlambda = 100)  # 交叉验证，如果结果不理想，可以重新单独运行这一行代码，或者换一下种子数plot(lasso_cv)  # 纵坐标：部分似然偏差。上横坐标：当前的变量数。下横坐标：log(λ)

coef_lasso <- coef(lasso, s = 0.1) # 查看每个变量的回归系数。s为lasso回归的结果中的lambda值，选取相应的λ值可以得到相应数量的变量的回归系数coef_lasso
set.seed(1234) # 设置种子数，保证每次得到的交叉验证图是一致的# 计算100个λ，画图，筛选表现最好的λ值lasso_cv <- cv.glmnet(x = x, y = y, family = 'binomial', alpha = 1, nlambda = 100) # 交叉验证，如果结果不理想，可以重新单独运行这一行代码，或者换一下种子数plot(lasso_cv) # 纵坐标：部分似然偏差。上横坐标：当前的变量数。下横坐标：log(λ)
lambda <- lasso_cv$lambda.min# lambda <- lasso_cv$lambda.1selambda
coef_lasso_cv <- coef(lasso, s = lambda) # 查看每个变量的回归系数。s为lasso回归结果中的lambda值，可以得到相应数量的变量的回归系数coef_lasso_cv[,1][coef_lasso_cv[,1]!=0]
exp(coef_lasso_cv[,1][coef_lasso_cv[,1]!=0]) # 自然对数的指数函数，它的定义为exp(x) = e^x，其中e是自然对数的底数（约等于2.718），回归系数<0的exp>1，反之<1。
model_lasso_min <- glmnet(x = x, y = y, alpha = 1, lambda = lasso_cv$lambda.min)# model_lasso_1se <- glmnet(x = x, y = y, alpha = 1, lambda = lasso_cv$lambda.1se)# 这两个值体现在参数lambda上。有了模型，可以将筛选的基因挑出来了。# 所有基因存放于model_lasso_min模型的子集beta中，用到的基因有一个s0值，没用的基因只记录了“.”，所以可以用下面代码挑出用到的基因。model_lasso_min$beta[,1][model_lasso_min$beta[,1]!=0]choose_gene_min = rownames(model_lasso_min$beta)[as.numeric(model_lasso_min$beta) != 0]# choose_gene_1se = rownames(model_lasso_1se$beta)[as.numeric(model_lasso_1se$beta) != 0]length(choose_gene_min)# length(choose_gene_1se)

最终筛选得到25个基因

DCBLD2        NECAB2          ECH1      SERPINE1  0.0001176182  0.0050609354  0.0013001252  0.0076438492        PAPPA2           EGF        IGFBP1          GPC3  0.0058767537  0.0017086554  0.0018483925  0.0031675436          SYN2          PRTG         H3C14            U3  0.0046942426  0.0172503772 -0.0108994643 -0.0100155497      CTAGE11P    AL355001.1    AL022316.1       PRDX3P2  0.0059144276 -0.0078457384 -0.0172823521 -0.0105412767      MTHFD1P1    AC002480.1         SCAT8          PWP2 -0.0035888050  0.0047360137  0.0010923467  0.0090971061     PTPRJ.AS1           HPR    AL049839.2      SOWAHCP2 -0.0011330955  0.0002045581  0.0045504295 -0.0013408436    AC016583.2  0.0054207717      

主要参考资料来自博客园，作者小高不高，链接：

https://www.cnblogs.com/xiaogaobugao/p/17131487.html。