听起来很霸气用起来并不难的随机森林

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

1.准备输入数据

输入数据是TCGA的表达矩阵expr、临床信息meta和group_list。保存为forest.Rdata了，在生信星球公众号后台聊天窗口回复“森林”即可获得。

1load("forest.Rdata")
2exprSet = expr[,group_list=="tumor"]
3
4## 必须保证生存资料和表达矩阵，两者一致
5all(substring(colnames(exprSet),1,12)==meta$ID)
6#> [1] TRUE

1library(randomForest)
2library(ROCR)
3library(genefilter)
4library(Hmisc)

2.构建随机森林模型

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

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

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

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

3.模型预测和评估

3.1自己预测自己

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

3.2 箱线图

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

1re=as.data.frame(re)
2colnames(re)=c('event','prob')
3re$event=as.factor(re$event)
4library(ggpubr) 
5p1 = ggboxplot(re, x = "event", y = "prob",
6               color = "event", palette = "jco",
7               add = "jitter")+ stat_compare_means()
8p1

对比lasso回归，这个似乎更好一些？但这个只是自己预测自己，太高也不是好事。还需要看预测集的结果。

3.3 AUC值

1library(ROCR)
2#library(caret)
3
4pred <- prediction(re[,2], re[,1])
5auc = performance(pred,"auc")@y.values[[1]]
6auc
7#> [1] 1

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

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

4.1 切割数据

用R包caret切割数据，生成的结果是一组代表列数的数字，用这些数字来给表达矩阵和meta取子集即可。

1library(caret)
2set.seed(12345679)
3sam<- createDataPartition(meta$event, p = .5,list = FALSE)
4train <- exprSet[,sam]
5test <- exprSet[,-sam]
6train_meta <- meta[sam,]
7test_meta <- meta[-sam,]

4.2 切割后的train数据集建模

和上面的建模方法一样。

 1#建模
 2x = t(log2(train+1))
 3y = train_meta$event
 4tmp_rf=file.path('./Rdata/TCGA_KIRC_miRNA_t_rf_output.Rdata')
 5if(!file.exists(tmp_rf)){
 6  rf_output=randomForest(x=x, y=y,importance = TRUE, ntree = 10001, proximity=TRUE )
 7  save(rf_output,file = tmp_rf)
 8}
 9load(file = tmp_rf)
10
11choose_gene=rownames(tail(rf_importances[order(rf_importances[,2]),],30))

5.模型预测和评估

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

 1x = t(log2(test+1))
 2y = test_meta$event
 3rf.prob <- predict(rf_output, x)
 4re=cbind(y ,rf.prob)
 5re=as.data.frame(re)
 6colnames(re)=c('event','prob')
 7re$event=as.factor(re$event)
 8library(ggpubr) 
 9p1 = ggboxplot(re, x = "event", y = "prob",
10               color = "event", palette = "jco",
11               add = "jitter")+ stat_compare_means()
12p1

再看AUC值。

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

emmm。。。还行。