小云昨天下班后，想著看看有啥生信好文章，看看能不能發(fā)現(xiàn)一些好套路來給小伙伴們謀福利，功夫不負(fù)有心人，還真的發(fā)現(xiàn)了一個(gè)。小云發(fā)現(xiàn)很多文章都喜歡用lasso回歸來進(jìn)行特征選擇。這是為什么呢？今天小云就以一篇Journal of Clinical Oncology雜志（IF=26.3）中的文章為例來說明這個(gè)原因。

這篇文章采用影響組學(xué)方法在結(jié)直腸癌淋巴結(jié)轉(zhuǎn)移方面進(jìn)行預(yù)測研究，目的是為了預(yù)測結(jié)直腸癌術(shù)前淋巴結(jié)轉(zhuǎn)移的風(fēng)險(xiǎn)，進(jìn)而為臨床醫(yī)生提供決策參考，具有重要的臨床價(jià)值。好文章一定要有好數(shù)據(jù)！這篇文章便意識到了這一點(diǎn)，于是作者在提取原始圖片數(shù)據(jù)的特征后，便考慮到了特征選擇。那么用什么方法最好呢？Lasso回歸！作者利用Lasso回歸從150個(gè)特征中篩選出24個(gè)回歸系數(shù)非零的具有預(yù)測作用的特征，在此基礎(chǔ)上完成后續(xù)實(shí)驗(yàn)。今天小云就用一個(gè)案例分析來帶大家認(rèn)識這個(gè)高分生信必備算法——Lasso回歸。

今天使用的數(shù)據(jù)集是一個(gè)前列腺癌數(shù)據(jù)集，該數(shù)據(jù)集有97個(gè)樣本，每個(gè)樣本包含9個(gè)特征。數(shù)據(jù)樣例如下圖所示。

其中train中的TRUE和FALSE表示該行樣本是否為訓(xùn)練集，lpsa為結(jié)果變量，PSA表示病人前列腺癌特異性抗原水平，lpsa即為其對數(shù)。本次分析即為利用lasso回歸提取出可以表現(xiàn)lpsa表達(dá)水平的特征。

分析開始

加載R包及數(shù)據(jù)

library(caret)

library(glmnet)

prostate1<-read.csv("prostate.csv")

head(prostate)

#測試集訓(xùn)練集劃分

train<-subset(prostate, train ==TRUE)[, 1:9]

str(train)

test <-subset(prostate, train==FALSE)[,1:9]

str(test)

#構(gòu)建Lasso回歸模型

x <-as.matrix(train[, 1:8])

y <-train[, 9]

#Lasso回歸，指定模型類型，這里采用gaussian類型

lasso<-?glmnet(x, y, family =“gaussian”, alpha =1)

print(lasso)

#查看系數(shù)和L1范數(shù)的關(guān)系

plot(ridge, label =TRUE)

#查看系數(shù)隨著lambda值變化的趨勢

plot(lasso, xvar =“l(fā)ambda”, label =TRUE)

#十折交叉驗(yàn)證，選擇合適的lambda值

cvfit <- cv.glmnet(x, y,nfolds = 10)

plot(cvfit)

#這里選擇lambda.lse，得到參數(shù)最少且性能優(yōu)良的模型

print(cvfit$lambda.1se) #0.3158532

#可以查看模型參數(shù)估計(jì)結(jié)果

coefficients<-coef(lasso,s= cvfit$lambda.1se)

Active.Index<-which(coefficients!=0) #系數(shù)不為0的特征索引

Active.coefficients<-coefficients[Active.Index]?#系數(shù)不為0的特征系數(shù)值

結(jié)果可以看到根據(jù)lasso回歸構(gòu)建模型， 發(fā)現(xiàn)lcavol、lweight、svi三個(gè)特征系數(shù)不為0，可以作為有效特征來代替病人樣本PSA表達(dá)水平，進(jìn)而展開后續(xù)分析。小伙伴們也可以將樣例數(shù)據(jù)換成自己的數(shù)據(jù)，利用lasso回歸提取出能夠替代病人臨床結(jié)果的有限特征，進(jìn)而展開更有意義的分析。高分文章近在咫尺！心動(dòng)的小伙伴一定不要錯(cuò)過！