github-actions[bot]
commented
3 years ago TCGA文献复现系列|与众不同的lasso回归，三大策略用CNS的方法做lasso回归 by 医学生之学习生信
常规型
Step1.安装加载R包
Step2.数据准备——因人而异
Step.3Lasso回归分析
Bootstrap法
方法由来及作用
代码实现
弹性网络回归
方法由来及作用
代码实现
参考资料
常规型
Step1.安装加载R包
#install.packages("glmnet")
library(glmnet)
library(survival)
library(dplyr)
Step2.数据准备——因人而异
目标
cliDat临床数据(matrix格式)
exprset预测变量数据(matrix格式)
代码实现
LIHC_cliDat <- data.table::fread("TCGA-LIHC.survival.tsv.gz",header = T)
load("LIHC_exprDat.Rdata")
gene <- c("TREM1","EXO1","COLEC10","C7","FCN3","DEFA6","KLK7")
tumorSample <- colnames(LIHC_exprDat)[as.numeric(substr(colnames(LIHC_exprDat),14,15))<10]
tumorSample <- substr(tumorSample,1,16)
#确保表达数据与临床信息样本名的一致性
LIHC_cliDat <- LIHC_cliDat[LIHC_cliDat$sample%in%tumorSample,]
LIHC_cliDat <- filter(LIHC_cliDat,OS.time>30)
colnames(LIHC_exprDat) <- substr(colnames(LIHC_exprDat),1,16)
exprset <- LIHC_exprDat[gene,LIHC_cliDat$sample]
exprset <- t(as.matrix(exprset))
exprset[1:4,1:4]
colnames(LIHC_cliDat) <- c("sample","fustat","_PATIENT","futime")
cliDat <- data.frame(time=LIHC_cliDat$futime,status=LIHC_cliDat$fustat)
#cliDat$status <- ifelse(cliDat$status=="Alive",0,1)
rownames(cliDat) <- LIHC_cliDat$sample
cliDat <- as.matrix(cliDat)
save(exprset,cliDat,file = "lassoData.Rdata")
说明
肝癌患者是文献中的“感兴趣结局事件”，所以必须排除正常患者 为了保证生存时间是疾病导致，而不是其他非疾病因素(如车祸等)导致，所以选择 OS.time>30，具体数值可因癌症类型或者主观判断抉择。(比如某些严重恶性肿瘤，本身生存时间比较短，如果筛选时间过长，很可能导致样本数目严重不足)最好将数据保存为Matrix格式，这是lasso回归要求的输入数据格式。当然，也可以在分析前 as.matrix
Step.3Lasso回归分析
目标：lasso回归两大目的，确定系数和λ值
代码实现
Part1.初步探索
rm(list = ls())
load("lassoData.Rdata")
x <- as.matrix(exprset)
y <- as.matrix(Surv(cliDat[,1],cliDat[,2]))
fit <- glmnet(x,y, family = "cox",alpha = 1,nlambda = 1000)
说明
x，表示预测变量(协变量)矩阵y, 因变量（生存对象矩阵）,注意创建生存对象时time不可以<=0!!!!!!family ，"gaussian"（高斯）意味y为线性数据（连续）,"binomial"("multinomial"),意味y为二分类变量（多分类变量）"cox"意味y为生存资料，y用Surv构建(包括"time"和"status"两列)alpha ，α=1表示lasso回归，α=0表示岭回归
nlambda，表示λ值的数量，默认为100，一般而言，但当%Dev值变化不大时，就会停止
print(fit)
# Df %Dev Lambda
# 1 0 0.0000000 0.14030
# 2 1 0.0002849 0.13900
# 3 1 0.0005652 0.13770
# 4 2 0.0008558 0.13650
# 5 2 0.0013460 0.13520
# 6 2 0.0018270 0.13400
# ...
# 260 7 0.0286300 0.01288
# 261 7 0.0286400 0.01276
# 262 7 0.0286400 0.01265
# 263 7 0.0286500 0.01253
# 264 7 0.0286600 0.01241
# 265 7 0.0286700 0.01230
说明
Df自由度(预示预测变量的数量)，%Dev解析残差百分比(越接近1越好)，Lambda表示λ值
plot(fit,xvar = "lambda",label = T)
plot(fit,xvar = "dev",label = T)
说明
xvar表示以什么作为x轴，label表示是否展示标签
coef(fit,s=0.0125)#预测系数
# 7 x 1 sparse Matrix of class "dgCMatrix"
# 1
# TREM1 0.141542140
# EXO1 0.160796869
# COLEC10 -0.017576378
# C7 -0.016905219
# FCN3 -0.032273735
# DEFA6 0.018806098
# KLK7 0.007648985
Part2.拟合λ值
cvfit <- cv.glmnet(x,y,
family = "cox",
nfolds = 10, #多少折检验，常用10折交叉检验
nlambda = 100,#λ数目，会影响图片竖杠数目
alpha=1)
plot(cvfit)
说明
注意：每次拟合的结果可能略有不同
plot(fit,xvar = "lambda",label = T)
abline(v = log(cvfit$lambda.min),col="#007dc3",lty=3,lwd=2)
text(log(cvfit$lambda.min),0.15,labels = paste0("λ.min=",round(cvfit$lambda.min,4)),
pos = 4,col ="#007dc3",cex = .7 )
abline(v = log(cvfit$lambda.1se),col="#b30838",lty=3,lwd=2)
text(log(cvfit$lambda.1se),0.01,labels = paste0("λ.1se=",round(cvfit$lambda.1se,4)),
pos = 2,col ="#b30838",cex = .7 )lambda.min(最小λ)：表示均方误差最小时的λ值lambda.1se(最小λ+标准误)：表示在不显著改变均方误差时尽可能减少待测特征的数量
Part3.根据拟合λ值，求系数
fitMin <- glmnet(x,y,family = "cox",lambda = cvfit$lambda.min,alpha = 1)
fitLse <- glmnet(x,y,family = "cox",lambda = cvfit$lambda.1se,alpha = 1)
cbind(fitLse$beta,fitMin$beta)#里面的beta等同于coef()的结果。
# 7 x 2 sparse Matrix of class "dgCMatrix"
# s0 s0
# TREM1 0.01138985 0.116068474
# EXO1 0.01896033 0.140059358
# COLEC10 . -0.006842854
# C7 . -0.007127897
# FCN3 . -0.013155785
# DEFA6 . .
# KLK7 . .
Part4.提取预测变量，预测风险评分（相关文章为何predict()函数计算的Riskscore不等于基因的表达量与其系数的乘积的加权呢?）
lassoGene<- rownames(fitMin$beta)[as.numeric(fitMin$beta)!=0]
predict(fit,newx = x,s = cvfit$lambda.min,type = "link")
说明
image-20210204132338317 验证一下：
sum(exprset[1,lassoGene]*fitMin$beta[as.numeric(fitMin$beta)!=0])#验证风险评分
Bootstrap法
方法由来及作用
作用：提高预测变量选择的一致性，使模型更具普适性
https://mp.weixin.qq.com/s/XwNVAWw21zy04lXFqjDIcw