其他
听起来很霸气用起来并不难的随机森林
大神一句话，菜鸟跑半年。我不是大神，但我可以缩短你走弯路的半年~
就像歌儿唱的那样，如果你不知道该往哪儿走，就留在这学点生信好不好~
这里有豆豆和花花的学习历程，从新手到进阶，生信路上有你有我！
花花写于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。。。还行。
插个小广告！
再给生信技能树打个call！