花花写于2020-02-02，看似雄赳赳气昂昂，其实不难。我们是R语言的应用者，不是开发者，所以诸如差异分析，富集分析，生存分析，绘图等等让初学者听起来遥不可及的操作，都只是一个R包，几个函数。打好基础，找到规律，任何函数都不难,不要被吓退。
今天是千年一遇的对称日，去年各地民政局都表示宠网民，表示今天上半天班给登记结婚（今天周日本来不用上班的），不料赶上疫情，又纷纷取消了今天的登记活动，真的是造化弄人啊。湖北省更是全面暂停了婚姻登记业务，不知道该说什么，武汉加油，希望一切快点好起来。
1.准备输入数据
输入数据是TCGA的表达矩阵expr、临床信息meta和group_list。保存为forest.Rdata了，在生信星球公众号后台聊天窗口回复“森林”即可获得。（为什么是森林呢，因为和随机森林图用的同一个数据，我又又又懒得改了。）
load("forest.Rdata")
exprSet = expr[,group_list=="tumor"]
## 必须保证生存资料和表达矩阵，两者一致
all(substring(colnames(exprSet),1,12)==meta$ID)
#> [1] TRUE
2.构建lasso回归模型
lasso回归的输入数据是表达矩阵(仅含tumor样本)和对应的生死。
x=t(log2(exprSet+1))
y=meta$event
library(glmnet)
model_lasso <- glmnet(x, y,nlambda=10, alpha=1)
print(model_lasso)
#>
#> Call: glmnet(x = x, y = y, alpha = 1, nlambda = 10)
#>
#> Df %Dev Lambda
#> [1,] 0 0.00000 0.127800
#> [2,] 7 0.08689 0.076610
#> [3,] 22 0.19490 0.045930
#> [4,] 49 0.28660 0.027530
#> [5,] 97 0.42080 0.016510
#> [6,] 164 0.55460 0.009895
#> [7,] 217 0.66240 0.005932
#> [8,] 275 0.74850 0.003556
#> [9,] 345 0.81660 0.002132
#> [10,] 403 0.86840 0.001278
这里是举例子，所以只计算了10个λ值，解释一下输出结果三列的意思。
- Df 是自由度
- 列%Dev代表了由模型解释的残差的比例，对于线性模型来说就是模型拟合的R^2(R-squred)。它在0和1之间，越接近1说明模型的表现越好，如果是0，说明模型的预测结果还不如直接把因变量的均值作为预测值来的有效。
- Lambda 是构建模型的重要参数。
解释的残差百分比越高越好，但是构建模型使用的基因的数量也不能太多，需要取一个折中值。
2.1挑选合适的λ值
计算1000个，画图，筛选表现最好的λ值
cv_fit <- cv.glmnet(x=x, y=y, nlambda = 1000,alpha = 1)
plot.cv.glmnet(cv_fit)
两条虚线分别指示了两个特殊的λ值,一个是lambda.min,一个是lambda.1se,这两个值之间的lambda都认为是合适的。lambda.1se构建的模型最简单，即使用的基因数量少，而lambda.min则准确率更高一点，使用的基因数量更多一点。
2.2 用这两个λ值重新建模
model_lasso_min <- glmnet(x=x, y=y, alpha = 1, lambda=cv_fit$lambda.min)
model_lasso_1se <- glmnet(x=x, y=y, alpha = 1, lambda=cv_fit$lambda.1se)
这两个值体现在参数lambda上。有了模型，可以将筛选的基因挑出来了。所有基因存放于模型的子集beta中，用到的基因有一个s0值，没用的基因只记录了“.”，所以可以用下面代码挑出用到的基因。
head(model_lasso_min$beta)
#> 6 x 1 sparse Matrix of class "dgCMatrix"
#> s0
#> hsa-let-7a-1 -0.0003827849
#> hsa-let-7a-2 .
#> hsa-let-7a-3 .
#> hsa-let-7b .
#> hsa-let-7c .
#> hsa-let-7d .
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)
#> [1] 33
length(choose_gene_1se)
#> [1] 12
3.模型预测和评估
3.1自己预测自己
newx参数是预测对象。输出结果lasso.prob是一个矩阵，第一列是min的预测结果，第二列是1se的预测结果，预测结果是概率，或者说百分比，不是绝对的0和1。
将每个样本的生死和预测结果放在一起，直接cbind即可。
lasso.prob <- predict(cv_fit, newx=x , s=c(cv_fit$lambda.min,cv_fit$lambda.1se) )
re=cbind(y ,lasso.prob)
head(re)
#> y 1 2
#> TCGA-A3-3307-01A-01T-0860-13 0 0.1345889 0.2609906
#> TCGA-A3-3308-01A-02R-1324-13 0 0.3984848 0.3547232
#> TCGA-A3-3311-01A-02R-1324-13 1 0.3165679 0.3206726
#> TCGA-A3-3313-01A-02R-1324-13 1 0.3996740 0.3487908
#> TCGA-A3-3316-01A-01T-0860-13 0 0.3221988 0.2928797
#> TCGA-A3-3317-01A-01T-0860-13 0 0.3740788 0.3733633
3.2 箱线图
对预测结果进行可视化。以实际的生死作为分组，画箱线图整体上查看预测结果。
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
可以看到，真实结果是死(0)的样本，预测的值就小一点(靠近0)，真实结果是活着(1)的样本，预测的值就大一点(靠近1)，整体上趋势是对的，但不是完全准确，模型是可用的。
对比两个λ值构建的模型，差别不大，1se的预测值不论生死都整体偏大一点点。
3.3 ROC曲线
计算AUC取值范围在0.5-1之间，越接近于1越好。可以根据预测结果绘制ROC曲线。
library(ROCR)
library(caret)
# 自己预测自己
#min
pred_min <- prediction(re[,2], re[,1])
auc_min = performance(pred_min,"auc")@y.values[[1]]
perf_min <- performance(pred_min,"tpr","fpr")
plot(perf_min,colorize=FALSE, col="blue")
lines(c(0,1),c(0,1),col = "gray", lty = 4 )
text(0.8,0.2, labels = paste0("AUC = ",round(auc_min,3)))
#1se
pred_1se <- prediction(re[,3], re[,1])
auc_1se = performance(pred_1se,"auc")@y.values[[1]]
perf_1se <- performance(pred_1se,"tpr","fpr")
plot(perf_1se,colorize=FALSE, col="red")
lines(c(0,1),c(0,1),col = "gray", lty = 4 )
text(0.8,0.2, labels = paste0("AUC = ",round(auc_min,3)))
- 强迫症选项，想把两个模型画一起。
plot(perf_min,colorize=FALSE, col="blue")
plot(perf_1se,colorize=FALSE, col="red",add = T)
lines(c(0,1),c(0,1),col = "gray", lty = 4 )
text(0.8,0.3, labels = paste0("AUC = ",round(auc_min,3)),col = "blue")
text(0.8,0.2, labels = paste0("AUC = ",round(auc_1se,3)),col = "red")
-还想用ggplot2画的更好看一点
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]],
tpr_1se = perf_1se@y.values[[1]],
fpr_1se = perf_1se@x.values[[1]])
library(ggplot2)
ggplot() +
geom_line(data = dat,aes(x = fpr_min, y = tpr_min),color = "blue") +
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 min = ",round(auc_min,2)),color = "blue")+
scale_x_continuous(name = "fpr")+
scale_y_continuous(name = "tpr")
ggplot() +
geom_line(data = dat,aes(x = fpr_min, y = tpr_min),color = "blue") +
geom_line(data = dat,aes(x = fpr_1se, y = tpr_1se),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 min = ",round(auc_min,2)),color = "blue")+
annotate("text",x = .75, y = .15,label = paste("AUC of 1se = ",round(auc_1se,2)),color = "red")+
scale_x_continuous(name = "fpr")+
scale_y_continuous(name = "tpr")
5.切割数据构建模型并预测
5.1 切割数据
用R包caret切割数据，生成的结果是一组代表列数的数字，用这些数字来给表达矩阵和meta取子集即可。
library(caret)
set.seed(12345679)
sam<- createDataPartition(meta$event, p = .5,list = FALSE)
head(sam)
#> Resample1
#> [1,] 5
#> [2,] 9
#> [3,] 13
#> [4,] 17
#> [5,] 19
#> [6,] 22
可查看两组一些临床参数切割比例
train <- exprSet[,sam]
test <- exprSet[,-sam]
train_meta <- meta[sam,]
test_meta <- meta[-sam,]
prop.table(table(train_meta$stage))
#>
#> i ii iii iv
#> 0.4636015 0.1072797 0.2796935 0.1494253
prop.table(table(test_meta$stage))
#>
#> i ii iii iv
#> 0.5249042 0.1111111 0.1954023 0.1685824
prop.table(table(test_meta$race))
#>
#> asian black or african american
#> 0.01171875 0.08593750
#> white
#> 0.90234375
prop.table(table(train_meta$race))
#>
#> asian black or african american
#> 0.01937984 0.13953488
#> white
#> 0.84108527
5.2 切割后的train数据集建模
和上面的建模方法一样。
#计算lambda
x = t(log2(train+1))
y = train_meta$event
cv_fit <- cv.glmnet(x=x, y=y, nlambda = 1000,alpha = 1)
plot.cv.glmnet(cv_fit)
#构建模型
model_lasso_min <- glmnet(x=x, y=y, alpha = 1, lambda=cv_fit$lambda.min)
model_lasso_1se <- glmnet(x=x, y=y, alpha = 1, lambda=cv_fit$lambda.1se)
#挑出基因
head(model_lasso_min$beta)
#> 6 x 1 sparse Matrix of class "dgCMatrix"
#> s0
#> hsa-let-7a-1 .
#> hsa-let-7a-2 .
#> hsa-let-7a-3 .
#> hsa-let-7b .
#> hsa-let-7c .
#> hsa-let-7d .
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)
#> [1] 48
length(choose_gene_1se)
#> [1] 3
4.模型预测
用训练集构建模型，预测测试集的生死，注意newx参数变了。
lasso.prob <- predict(cv_fit, newx=t(log2(test+1)), s=c(cv_fit$lambda.min,cv_fit$lambda.1se) )
re=cbind(test_meta$event ,lasso.prob)
head(re)
#> 1 2
#> TCGA-A3-3307-01A-01T-0860-13 0 0.2760721 0.3084088
#> TCGA-A3-3308-01A-02R-1324-13 0 0.3646723 0.3370767
#> TCGA-A3-3311-01A-02R-1324-13 1 0.3190559 0.3292243
#> TCGA-A3-3313-01A-02R-1324-13 1 0.3427763 0.3433195
#> TCGA-A3-3317-01A-01T-0860-13 0 0.3489759 0.3589717
#> TCGA-A3-3319-01A-02R-1324-13 0 0.7132020 0.3882291
再次进行可视化
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
目测还行哈？再看AUC值。
library(ROCR)
library(caret)
# 训练集模型预测测试集
#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")
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]],
tpr_1se = perf_1se@y.values[[1]],
fpr_1se = perf_1se@x.values[[1]])
ggplot() +
geom_line(data = dat,aes(x = fpr_min, y = tpr_min),color = "blue") +
geom_line(data = dat,aes(x = fpr_1se, y = tpr_1se),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 min = ",round(auc_min,2)),color = "blue")+
annotate("text",x = .75, y = .15,label = paste("AUC of 1se = ",round(auc_1se,2)),color = "red")+
scale_x_continuous(name = "fpr")+
scale_y_continuous(name = "tpr")
主要看AUC值，可以说一般，也可以说还行。
网友评论