1.準(zhǔn)備輸入數(shù)據(jù)
load("TCGA-KIRC_sur_model.Rdata")
2.構(gòu)建lasso回歸模型
輸入數(shù)據(jù)是表達(dá)矩陣(僅含tumor樣本)和對(duì)應(yīng)的生死。
x=t(exprSet)
y=meta$event
library(glmnet)
cv_fit <- cv.glmnet(x=x, y=y, nlambda = 1000,alpha = 1)
model_lasso_min <- glmnet(x=x, y=y, alpha = 1, lambda=cv_fit$lambda.min)
choose_gene_mi n=rownames(model_lasso_min$beta)[as.numeric(model_lasso_min$beta)!=0]
length(choose_gene_min)
#> [1] 39
3.模型預(yù)測(cè)和評(píng)估
3.1自己預(yù)測(cè)自己
lasso.prob <- predict(cv_fit, newx=x , s=cv_fit$lambda.min )
re=cbind(y ,lasso.prob)
head(re)
#> y 1
#> TCGA-A3-3307-01A-01T-0860-13 0 0.1107463
#> TCGA-A3-3308-01A-02R-1324-13 0 0.3985909
#> TCGA-A3-3311-01A-02R-1324-13 1 0.2875707
#> TCGA-A3-3313-01A-02R-1324-13 1 0.3884061
#> TCGA-A3-3316-01A-01T-0860-13 0 0.3288315
#> TCGA-A3-3317-01A-01T-0860-13 0 0.3793098
3.2 time-ROC
new_dat=meta
library(timeROC)
library(survival)
library(survminer)
new_dat$fp=as.numeric(lasso.prob[,1])
with(new_dat,
ROC <<- timeROC(T=time,#結(jié)局時(shí)間
delta=event,#生存結(jié)局
marker=fp,#預(yù)測(cè)變量
cause=1,#陽(yáng)性結(jié)局賦值悬槽，比如死亡與否
weighting="marginal",#權(quán)重計(jì)算方法，marginal是默認(rèn)值瞬浓，采用km計(jì)算刪失分布
times=c(60,100),#時(shí)間點(diǎn)初婆，選取5年(60個(gè)月)和8年生存率，小的年份寫前面
ROC = TRUE,
iid = TRUE)
)
auc_60 = ROC$AUC[[1]]
auc_100 = ROC$AUC[[2]]
dat = data.frame(tpr_60 = ROC$TP[,1],
fpr_60 = ROC$FP[,1],
tpr_100 = ROC$TP[,2],
fpr_100 = ROC$FP[,2])
library(ggplot2)
ggplot() +
geom_line(data = dat,aes(x = fpr_60, y = tpr_60),color = "blue") +
geom_line(data = dat,aes(x = fpr_100, y = tpr_100),color = "red")+
geom_line(aes(x=c(0,1),y=c(0,1)),color = "grey")+
theme_bw()+
annotate("text",x = .75, y = .25,
label = paste("AUC of 60 = ",round(auc_60,2)),color = "blue")+
annotate("text",x = .75, y = .15,label = paste("AUC of 100 = ",round(auc_100,2)),color = "red")+
scale_x_continuous(name = "fpr")+
scale_y_continuous(name = "tpr")
4.切割數(shù)據(jù)構(gòu)建模型并預(yù)測(cè)
4.1 切割數(shù)據(jù)
用R包c(diǎn)aret切割數(shù)據(jù)猿棉，生成的結(jié)果是一組代表列數(shù)的數(shù)字磅叛，用這些數(shù)字來(lái)給表達(dá)矩陣和meta取子集即可。
library(caret)
set.seed(12345679)
sam<- createDataPartition(meta$event, p = .5,list = FALSE)
train <- exprSet[,sam]
test <- exprSet[,-sam]
train_meta <- meta[sam,]
test_meta <- meta[-sam,]
4.2 切割后的train數(shù)據(jù)集建模
和上面的建模方法一樣萨赁。
#計(jì)算lambda
x = t(train)
y = train_meta$event
cv_fit <- cv.glmnet(x=x, y=y, nlambda = 1000,alpha = 1)
#構(gòu)建模型
model_lasso_min <- glmnet(x=x, y=y, alpha = 1, lambda=cv_fit$lambda.min)
#挑出基因
choose_gene_min=rownames(model_lasso_min$beta)[as.numeric(model_lasso_min$beta)!=0]
length(choose_gene_min)
#> [1] 18
4.3 模型預(yù)測(cè)
用訓(xùn)練集構(gòu)建模型弊琴，預(yù)測(cè)測(cè)試集的生死，注意newx參數(shù)變了杖爽。
lasso.prob <- predict(cv_fit, newx=t(test), s=cv_fit$lambda.min)
re=cbind(test_meta$event ,lasso.prob)
4.4 time-ROC
new_dat = test_meta
library(timeROC)
library(survival)
library(survminer)
new_dat$fp=as.numeric(lasso.prob[,1])
with(new_dat,
ROC <<- timeROC(T=time,#結(jié)局時(shí)間
delta=event,#生存結(jié)局
marker=fp,#預(yù)測(cè)變量
cause=1,#陽(yáng)性結(jié)局賦值敲董，比如死亡與否
weighting="marginal",#權(quán)重計(jì)算方法，marginal是默認(rèn)值慰安，采用km計(jì)算刪失分布
times=c(60,100),#時(shí)間點(diǎn)腋寨，選取5年(60個(gè)月)和8年生存率
ROC = TRUE,
iid = TRUE)
)
auc_60 = ROC$AUC[[1]]
auc_100 = ROC$AUC[[2]]
dat = data.frame(tpr_60 = ROC$TP[,1],
fpr_60 = ROC$FP[,1],
tpr_100 = ROC$TP[,2],
fpr_100 = ROC$FP[,2])
library(ggplot2)
ggplot() +
geom_line(data = dat,aes(x = fpr_60, y = tpr_60),color = "blue") +
geom_line(data = dat,aes(x = fpr_100, y = tpr_100),color = "red")+
geom_line(aes(x=c(0,1),y=c(0,1)),color = "grey")+
theme_bw()+
annotate("text",x = .75, y = .25,
label = paste("AUC of 60 = ",round(auc_60,2)),color = "blue")+
annotate("text",x = .75, y = .15,label = paste("AUC of 100 = ",round(auc_100,2)),color = "red")+
scale_x_continuous(name = "fpr")+
scale_y_continuous(name = "tpr")
*生信技能樹(shù)課程筆記
