很多時(shí)候遇到的數(shù)據(jù)均為表達(dá)譜數(shù)據(jù)，所以基因數(shù)量很多，需要用lasso回歸來(lái)篩選變量。一般都是使用R來(lái)實(shí)現(xiàn)，下述代碼為篩選變量與篩選出來(lái)的變量做ROC曲線
rm(list = ls())
options(stringsAsFactors = F)
load("sur_model.Rdata")
load("UCECcox_lr.Rdata")
load('UCECcox.Rdata')
load('UCEClog_rank_p.Rdata')
cox1=c('RP4-616B8.5','RP11-389G6.3',"AP000696.2",'CTD-2377D24.6',"LINC01559" ,"LINC00629","AC005062.2","LINC01018")
##A='CTD-2377D24.6',B='RP4-616B8.5',C='RP11-389G6.3'
##'CTD-2377D24.6','RP4-616B8.5','RP11-389G6.3'
##cox1為定義的任何向量，比如里面含有你所感興趣的基因名稱(chēng)
exprSet = exprSet[cox1,]#表達(dá)譜數(shù)據(jù)中提取你所感興趣的基因表達(dá)
exprSet=as.data.frame(exprSet)
dim(exprSet)
exprSet[1:2,1:2]#檢查你的數(shù)據(jù)是否正確
x=t(exprSet)
y=meta$event#提取病例的預(yù)后情況（0與1）列，meta是臨床數(shù)據(jù)
#lasso
library(glmnet)
model_lasso <- glmnet(x, y,nlambda=10, alpha=1)
print(model_lasso)
set.seed(13098)
cv_fit <- cv.glmnet(x=x, y=y, nlambda = 1000,alpha = 1)
plot(cv_fit)
#ROC
#load('for_cox.Rdata')
lasso.prob <- predict(cv_fit, newx=x , s=c(cv_fit$lambda.min,cv_fit$lambda.1se) )
head(lasso.prob)
re=cbind(y ,lasso.prob)
head(re)
re1=as.data.frame(re)
write.csv(re,file='')
#箱線圖（可有無(wú)
re=as.data.frame(re)
colnames(re)=c('event','prob_min','prob_1se')
re$event=as.factor(re$event)
library(ggpubr)
p1 = ggboxplot(re, x = "event", y = "prob_min",
color = "event", palette = "jco",
add = "jitter")+ stat_compare_means()
p2 = ggboxplot(re, x = "event", y = "prob_1se",
color = "event", palette = "jco",
add = "jitter")+ stat_compare_means()
library(patchwork)
p1+p2
#ROC
library(ROCR)
library(caret)
# 自己預(yù)測(cè)自己
#min
pred_min <- prediction(re[,2], re[,1])
auc_min = performance(pred_min,"auc")@y.values[[1]]
perf_min <- performance(pred_min,"tpr","fpr")
#1se
pred_1se <- prediction(re[,3], re[,1])
auc_1se = performance(pred_1se,"auc")@y.values[[1]]
perf_1se <- performance(pred_1se,"tpr","fpr")
#ggpot2
tpr_min = performance(pred_min,"tpr")@y.values[[1]]
tpr_1se = performance(pred_1se,"tpr")@y.values[[1]]
dat = data.frame(tpr_min = perf_min@y.values[[1]],
fpr_min = perf_min@x.values[[1]]
)
library(ggplot2)
ggplot() +
geom_line(data = dat,aes(x = fpr_min, y = tpr_min),color = "blue",size=3) +
geom_line(aes(x=c(0,1),y=c(0,1)),color = "grey",size = 3,linetype=6 )+
theme_bw()+
annotate("text",x = .75, y = .25,
label = paste("AUC of min = ",round(auc_min,2)),color = "blue",size=10)+
scale_x_continuous(name = "1-Specificity")+
scale_y_continuous(name = "Sensitivity")+
theme(plot.title = element_text(hjust = 0.5,size = 24),
axis.text=element_text(size=24),
axis.title.x = element_text(size = 24),
axis.title.y = element_text(size = 24),
legend.title = element_text(size=24),
legend.text = element_text(size=18))
dev.off
有很多時(shí)候需要得到一個(gè)gene對(duì)結(jié)局貢獻(xiàn)性的ROC，及其各個(gè)ROC曲線之間比較，下述代碼為SAS程序代碼
data roc;
input alb tp totscore popind @@;##alb tp totscore 均為自變量，比如一個(gè)gene的表達(dá)量，popind為結(jié)局（0 and 1）
totscore = 10 - totscore;##對(duì)自變量按照需要處理
datalines;
3.0 5.8 10 0 3.2 6.3 5 1 3.9 6.8 3 1 2.8 4.8 6 0
3.2 5.8 3 1 0.9 4.0 5 0 2.5 5.7 8 0 1.6 5.6 5 1
3.8 5.7 5 1 3.7 6.7 6 1 3.2 5.4 4 1 3.8 6.6 6 1
4.1 6.6 5 1 3.6 5.7 5 1 4.3 7.0 4 1 3.6 6.7 4 0
2.3 4.4 6 1 4.2 7.6 4 0 4.0 6.6 6 0 3.5 5.8 6 1
3.8 6.8 7 1 3.0 4.7 8 0 4.5 7.4 5 1 3.7 7.4 5 1
3.1 6.6 6 1 4.1 8.2 6 1 4.3 7.0 5 1 4.3 6.5 4 1
3.2 5.1 5 1 2.6 4.7 6 1 3.3 6.8 6 0 1.7 4.0 7 0
3.7 6.1 5 1 3.3 6.3 7 1 4.2 7.7 6 1 3.5 6.2 5 1
2.9 5.7 9 0 2.1 4.8 7 1 2.8 6.2 8 0 4.0 7.0 7 1
3.3 5.7 6 1 3.7 6.9 5 1 3.6 6.6 5 1
;
ods graphics on;
proc logistic data=roc plots=roc(id=prob);
model popind(event='0') = alb tp totscore / nofit;
roc 'Albumin' alb;
roc 'K-G Score' totscore;
roc 'Total Protein' tp;
roccontrast reference('K-G Score') / estimate e;
run;
ods graphics off;
