對(duì)于生信實(shí)操,，我們認(rèn)為,，在線工具和 R 語(yǔ)言各有千秋,，兩者結(jié)合或許才是最好的策略,。這句話對(duì)于生信學(xué)習(xí),，或者數(shù)據(jù)挖掘都適用。純生信操作簡(jiǎn)單,，適合入門和初步數(shù)據(jù)分析,；R 門檻稍高，個(gè)性化強(qiáng),，適合深度數(shù)據(jù)挖掘,。為此，在純生信教程基礎(chǔ)上,，我們推出在線工具與 R 結(jié)合的新版21天教程（重磅推出 | 新版21天生信實(shí)操教程！）,。新版21天教程將以原純生信教程以及《癌生物學(xué)》,、《R 數(shù)據(jù)科學(xué)》和《詹韋免疫學(xué)》等經(jīng)典權(quán)威教材的核心內(nèi)容為基礎(chǔ)，并結(jié)合腫瘤和非腫瘤疾病最新發(fā)表的論文展開,。在21天學(xué)習(xí)結(jié)束后,，芒果團(tuán)隊(duì)會(huì)集中展開指定文獻(xiàn)中所有圖表的復(fù)現(xiàn)，讓參與學(xué)習(xí)的小伙伴真正實(shí)現(xiàn)“學(xué)了就能會(huì),、會(huì)了就能用”,。這將是新版21天教程的最大特色！盡管本文沒有進(jìn)行Lasso回歸模型構(gòu)建,，但其應(yīng)用廣泛,，值得關(guān)注。本次我們主要分享Lasso回歸,。在與果友交流過(guò)程中,，Lasso回歸經(jīng)常是與Cox回歸同時(shí)被問(wèn)及的話題,。比如①預(yù)后相關(guān)基因是指做單因素Cox回歸有意義的基因，還要深入做Lasso或者多因素Cox嗎,？要?、贚asso預(yù)后建模后得到的風(fēng)險(xiǎn)評(píng)估基因熱圖的差異不顯著該怎么處理呢？AUC明明很高,，但具體到模型內(nèi)的各個(gè)基因在高低風(fēng)險(xiǎn)中的差異表達(dá)無(wú)顯著區(qū)別是怎么回事,？差異分析是前提！③預(yù)測(cè)模型做完單因素后可以直接做Lasso回歸嗎,，就不做多因素了,，因?yàn)長(zhǎng)asso回歸的基因數(shù)只有幾個(gè)了,？④單因素Cox?Lasso篩選預(yù)后相關(guān)DEGs?多因素模型風(fēng)險(xiǎn)評(píng)分,。這種模式能否只做Lasso?多因素呢？因?yàn)槲覇位蚝Y選的基因過(guò)少,，直接做lasso模型效果好一點(diǎn),？群里討論了,，不過(guò)常規(guī)是都要做，或者重置差異分析的條件,。Tibshirani在1996年引入 Lasso (Least Absolute Shrinkage and Selection Operator)的概念,，用于參數(shù)的選擇和收縮。Lasso回歸是以縮小變量集為核心思想的壓縮估計(jì)方法,，通過(guò)構(gòu)造一個(gè)懲罰函數(shù),，將變量系數(shù)進(jìn)行壓縮，并使某些不重要的標(biāo)量回歸系數(shù)變?yōu)?,，從而選擇最重要的變量進(jìn)入下游分析,。Lasso 回歸是一種線性回歸模型，也稱正則化或懲罰回歸模型,，使用L1正則化項(xiàng)來(lái)限制變量的系數(shù),，從而用于變量數(shù)量較多的大數(shù)據(jù)集，而傳統(tǒng)的線性回歸模型無(wú)法處理這類大數(shù)據(jù),。交叉驗(yàn)證用于確證Lasso結(jié)果中的取值是否最佳,。Lasso回歸和嶺回歸都是常用的正則化技術(shù)，用來(lái)減少系數(shù)的大小,，但使用的正則項(xiàng)不同,。Lasso回歸使用L1正則項(xiàng)，嶺回歸使用L2正則項(xiàng),。L1正則項(xiàng)的特點(diǎn)是,，它可以完全移除某些系數(shù)，使得模型中的變量更加稀疏,，即變量較少,，這對(duì)于特征選擇非常有用。而嶺回歸不具有變量的篩選功能,，它只能讓變量之間的差異變得更小,。但L1正則項(xiàng)也有一些缺點(diǎn)，由于L1正則項(xiàng)對(duì)所有變量的懲罰都是絕對(duì)值,，所以對(duì)于變量之間的大小沒有影響,。因此，L1正則項(xiàng)會(huì)使得變量之間存在較大的差異,，導(dǎo)致模型的可解釋性降低,。接下來(lái)，我們進(jìn)行Lasso回歸的實(shí)操,。首先準(zhǔn)備數(shù)據(jù),。
####----11_Lasso_Regression----####rm(list = ls())gc()
load('TCGA_STAD_data.RData')
library(dplyr)library(tidyverse)library(ggplot2)library(survival)
library(glmnet) # 正則化技術(shù)library(car) # 多重共線性檢驗(yàn)library(corrplot) # 相關(guān)性繪圖library(ggpubr) #基于ggplot2的可視化包ggpubrlibrary(patchwork) # 拼圖library(ROCR) # 畫ROC曲線library(caret) # 切割數(shù)據(jù)
Data <- data[,c(6,12:30917)]Data$Status = ifelse(Data$Status =='Dead','1','0')Data <- as.data.frame(lapply(Data,as.numeric))
### 數(shù)據(jù)檢查min(Data$month) # 查看最小值，看是否有0Data2 <- filter(Data, Data$month > 0) # 刪除數(shù)據(jù)框中time值小于0的樣本，有0的話后續(xù)分析會(huì)報(bào)錯(cuò)dataL <- Data2 # 數(shù)據(jù)備份
探究數(shù)據(jù)的特征,， 繪制直方圖和正態(tài)分布概率密度函數(shù)圖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') # 正態(tài)分布的概率密度函數(shù)# shapiro-wilk檢驗(yàn)shapiro.test(dataL$FBLN5) # 適用于小樣本數(shù)據(jù),，N≤50，若p值<0.05,，則數(shù)據(jù)不符合正態(tài)分布# Kolmogorov-Smirnov檢驗(yàn)ks.test(dataL$FBLN5, pnorm, mean = mean(dataL$FBLN5), sd = sd(dataL$FBLN5)) # 適用于大樣本數(shù)據(jù),，N>50，若p值<0.05,，則數(shù)據(jù)不符合正態(tài)分布# 對(duì)所有列的數(shù)據(jù)進(jìn)行Kolmogorov-Smirnov循環(huán)檢驗(yàn)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,，則數(shù)據(jù)不符合正態(tài)分布
繪制相關(guān)矩陣，檢查數(shù)據(jù)的多重共線性問(wèn)題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')
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)# 當(dāng)VIF大于5時(shí),，認(rèn)為存在多重共線性,，當(dāng)VIF大于10時(shí)，認(rèn)為存在嚴(yán)重的多重共線性
構(gòu)建Lasso回歸模型time <- dataL[,'month']status <- dataL[,'Status']x <- as.matrix(dataL[,-c(1,2)]) #x為輸入特征,是矩陣格式y(tǒng) <- as.matrix(dataL$Status)lasso <- glmnet(x = x, y = y, family = 'binomial', alpha = 1, # alpha = 1為L(zhǎng)ASSO回歸,，= 0為嶺回歸,，0和1之間則為彈性網(wǎng)絡(luò) nlambda = 100) # nlambda表示正則化路徑中的個(gè)數(shù)，這個(gè)參數(shù)就可以起到一個(gè)閾值的作用,，決定有多少基因的系數(shù)可以留下來(lái),。默認(rèn)值為100。print(lasso)plot(lasso, xvar = 'lambda', label = TRUE) # 系數(shù)分布圖,，由“l(fā)og-lambda”與變量系數(shù)作圖,，展示根據(jù)lambda的變化情況每一個(gè)特征的系數(shù)變化，展示了Lasso回歸篩選變量的動(dòng)態(tài)過(guò)程plot(lasso, xvar = 'dev', label = TRUE) # 也可以對(duì)%dev繪圖plot(lasso, xvar = 'norm', label = TRUE) # “L1范數(shù)”與變量系數(shù)作圖
交叉驗(yàn)證,，挑選合適的λ值,，變量篩選set.seed(1234) # 設(shè)置種子數(shù)，保證每次得到的交叉驗(yàn)證圖是一致的# 計(jì)算100個(gè)λ,，畫圖,，篩選表現(xiàn)最好的λ值lasso_cv <- cv.glmnet(x = x, y = y, family = 'binomial', alpha = 1, nlambda = 100) # 交叉驗(yàn)證，如果結(jié)果不理想,，可以重新單獨(dú)運(yùn)行這一行代碼,，或者換一下種子數(shù)plot(lasso_cv) # 縱坐標(biāo)：部分似然偏差。上橫坐標(biāo)：當(dāng)前的變量數(shù),。下橫坐標(biāo)：log(λ)
 coef_lasso <- coef(lasso, s = 0.1) # 查看每個(gè)變量的回歸系數(shù),。s為lasso回歸的結(jié)果中的lambda值，選取相應(yīng)的λ值可以得到相應(yīng)數(shù)量的變量的回歸系數(shù)coef_lasso
set.seed(1234) # 設(shè)置種子數(shù),，保證每次得到的交叉驗(yàn)證圖是一致的# 計(jì)算100個(gè)λ,，畫圖，篩選表現(xiàn)最好的λ值lasso_cv <- cv.glmnet(x = x, y = y, family = 'binomial', alpha = 1, nlambda = 100) # 交叉驗(yàn)證,，如果結(jié)果不理想,，可以重新單獨(dú)運(yùn)行這一行代碼,，或者換一下種子數(shù)plot(lasso_cv) # 縱坐標(biāo)：部分似然偏差。上橫坐標(biāo)：當(dāng)前的變量數(shù),。下橫坐標(biāo)：log(λ)
lambda <- lasso_cv$lambda.min# lambda <- lasso_cv$lambda.1selambda
coef_lasso_cv <- coef(lasso, s = lambda) # 查看每個(gè)變量的回歸系數(shù)。s為lasso回歸結(jié)果中的lambda值,，可以得到相應(yīng)數(shù)量的變量的回歸系數(shù)coef_lasso_cv[,1][coef_lasso_cv[,1]!=0]
exp(coef_lasso_cv[,1][coef_lasso_cv[,1]!=0]) # 自然對(duì)數(shù)的指數(shù)函數(shù),，它的定義為exp(x) = e^x，其中e是自然對(duì)數(shù)的底數(shù)（約等于2.718）,，回歸系數(shù)<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)# 這兩個(gè)值體現(xiàn)在參數(shù)lambda上,。有了模型,，可以將篩選的基因挑出來(lái)了。# 所有基因存放于model_lasso_min模型的子集beta中,，用到的基因有一個(gè)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)
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
主要參考資料來(lái)自博客園,，作者小高不高，鏈接： https://www.cnblogs.com/xiaogaobugao/p/17131487.html,。
|
