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
*生信技能树课程笔记
