九、随机森林

1.准备输入数据

输入数据是肿瘤样本表达矩阵exprSet、临床信息meta

load("TCGA-KIRC_sur_model.Rdata")
library(randomForest)
library(ROCR)
library(genefilter)
library(Hmisc)

2.构建随机森林模型

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

x=t(exprSet)
y=meta$event
#构建模型，一个叫randomForest的函数,运行时间很长,存Rdata跳过
tmp_rf="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 1.852761e-04     0.1787383
#> hsa-let-7a-2 2.167420e-04     0.1916623
#> hsa-let-7a-3 2.218169e-04     0.1858544
#> hsa-let-7b   7.399404e-05     0.1628646
#> hsa-let-7c   7.658155e-05     0.1635053
#> hsa-let-7d   1.974099e-04     0.2382185
#解释量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.1364447
#> TCGA-A3-3308-01A-02R-1324-13 0 0.1793771
#> TCGA-A3-3311-01A-02R-1324-13 1 0.6709712
#> TCGA-A3-3313-01A-02R-1324-13 1 0.7742376
#> TCGA-A3-3316-01A-01T-0860-13 0 0.2035863
#> TCGA-A3-3317-01A-01T-0860-13 0 0.1619938

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")+ 
  scale_y_continuous(limits = c(0,1)) +
  stat_compare_means()
p1

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

3.3 ROC曲线

library(ROCR)
#library(caret)

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

自己预测自己，auc值竟然是1。

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(train)
y = train_meta$event
tmp_rf="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))
head(choose_gene)
#> [1] "hsa-mir-511-1"  "hsa-mir-155"    "hsa-mir-409"    "hsa-mir-1185-1"
#> [5] "hsa-mir-1277"   "hsa-mir-149"

5.模型预测

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

x = t(test)
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")+ 
  scale_y_continuous(limits = c(0,1)) +
  stat_compare_means()
p1

再看AUC值。

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

*生信技能树课程笔记