TCGA-11.听起来很霸气用起来并不难的随机森林

花花写于2020-02-05 今天豆豆点了个外卖，以前可以送到房间门口，现在只能去大门口取了，而且回来还需要例行检查。感谢疫情中守护我们的工作人员，愿一切快点好起来。
今天继续学习机器学习算法做生存模型，高大上的随机森林，听起来需要学两年，但初步应用的代码很简单，并不是搞懂了算法、原理才能去用。这个运行时间很长，有点考验电脑哦。

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

library(randomForest)
library(ROCR)
library(genefilter)
library(Hmisc)

2.构建随机森林模型

输入数据是表达矩阵(仅含tumor样本)和对应的生死。

高能：此处可能需要运行十几分钟，取决于你电脑的配置，配置很不好就不要跑这一步了，容易卡死，看看代码理解一下意思就好。

x=t(log2(exprSet+1))
y=meta$event
#构建模型，一个叫randomForest的函数,运行时间很长,存Rdata跳过
tmp_rf=file.path('./Rdata/TCGA_KIRC_miRNA_rf_output.Rdata')
if(!file.exists(tmp_rf)){
  rf_output=randomForest(x=x, y=y,importance = TRUE, ntree = 10001, proximity=TRUE )
  save(rf_output,file = tmp_rf)
}
load(file = tmp_rf)
#top30的基因
varImpPlot(rf_output, type=2, n.var=30, scale=FALSE, 
           main="Variable Importance (Gini) for top 30 predictors",cex =.7)

rf_importances=importance(rf_output, scale=FALSE)
head(rf_importances)
#>                   %IncMSE IncNodePurity
#> hsa-let-7a-1 0.0005839108     0.3952803
#> hsa-let-7a-2 0.0006337582     0.4783033
#> hsa-let-7a-3 0.0005867823     0.4106104
#> hsa-let-7b   0.0005974527     0.2871573
#> hsa-let-7c   0.0003302193     0.2344674
#> hsa-let-7d   0.0003187660     0.1929435
#解释量top30的基因,和图上是一样的，从下往上看。
choose_gene=rownames(tail(rf_importances[order(rf_importances[,2]),],30))

3.模型预测和评估

3.1自己预测自己

rf.prob <- predict(rf_output, x)
re=cbind(y ,rf.prob)
head(re)
#>                              y   rf.prob
#> TCGA-A3-3307-01A-01T-0860-13 0 0.1094791
#> TCGA-A3-3308-01A-02R-1324-13 0 0.1719678
#> TCGA-A3-3311-01A-02R-1324-13 1 0.6914492
#> TCGA-A3-3313-01A-02R-1324-13 1 0.7408259
#> TCGA-A3-3316-01A-01T-0860-13 0 0.1635503
#> TCGA-A3-3317-01A-01T-0860-13 0 0.1204913

3.2 箱线图

对预测结果进行可视化。以实际的生死作为分组，画箱线图整体上查看预测结果。

re=as.data.frame(re)
colnames(re)=c('event','prob')
re$event=as.factor(re$event)
library(ggpubr) 
p1 = ggboxplot(re, x = "event", y = "prob",
               color = "event", palette = "jco",
               add = "jitter")+ stat_compare_means()
p1

对比lasso回归，这个似乎更好一些？

3.3 AUC值

library(ROCR)
#library(caret)

pred <- prediction(re[,2], re[,1])
auc = performance(pred,"auc")@y.values[[1]]
auc
#> [1] 1

自己预测自己，auc值竟然是1。ROC曲线我就不画了，和前面的几个模型画的代码一样。

4.切割数据构建模型并预测

4.1 切割数据

用R包caret切割数据，生成的结果是一组代表列数的数字，用这些数字来给表达矩阵和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数据集建模

和上面的建模方法一样。

#建模
x = t(log2(train+1))
y = train_meta$event
tmp_rf=file.path('./Rdata/TCGA_KIRC_miRNA_t_rf_output.Rdata')
if(!file.exists(tmp_rf)){
  rf_output=randomForest(x=x, y=y,importance = TRUE, ntree = 10001, proximity=TRUE )
  save(rf_output,file = tmp_rf)
}
load(file = tmp_rf)

choose_gene=rownames(tail(rf_importances[order(rf_importances[,2]),],30))
length(choose_gene)
#> [1] 30

5.模型评估和预测

用训练集构建模型，预测测试集的生死。

x = t(log2(test+1))
y = test_meta$event
rf.prob <- predict(rf_output, x)
re=cbind(y ,rf.prob)
re=as.data.frame(re)
colnames(re)=c('event','prob')
re$event=as.factor(re$event)
library(ggpubr) 
p1 = ggboxplot(re, x = "event", y = "prob",
               color = "event", palette = "jco",
               add = "jitter")+ stat_compare_means()
p1

再看AUC值。

library(ROCR)
# 训练集模型预测测试集
pred <- prediction(re[,2], re[,1])
auc= performance(pred,"auc")@y.values[[1]]
auc
#> [1] 0.7118081

还行。

微信公众号生信星球同步更新我的文章，欢迎大家扫码关注！

我们有为生信初学者准备的学习小组,点击查看◀️
想要参加我的线上线下课程，也可加好友咨询🔼
如果需要提问，请先看生信星球答疑公告