当前位置：首页 > news >正文

2025 - 生信信息学 - GEO数据分析 - RF分析（随机森林）

news 2025/6/29 5:49:19

GEO数据分析 - RF分析（随机森林）

在这里插入图片描述

01 准备数据文件

在这里插入图片描述

#install.packages("randomForest")#引用包
library(randomForest)
set.seed(123456)inputFile="diffGeneExp.txt"       #输入文件
setwd("/Users/wangyang/Desktop/BCBM/02randomForest-ANN")      #设置工作目录#读取输入文件
data=read.table(inputFile, header=T, sep="\t", check.names=F, row.names=1)
data=t(data)
group=gsub("(.*)\\_(.*)", "\\2", row.names(data))#随机森林树
rf=randomForest(as.factor(group)~., data=data, ntree=500)
pdf(file="forest.pdf", width=6, height=6)
plot(rf, main="Random forest", lwd=2)
dev.off()#找出误差最小的点
optionTrees=which.min(rf$err.rate[,1])
optionTrees
rf2=randomForest(as.factor(group)~., data=data, ntree=optionTrees)#查看基因的重要性
importance=importance(x=rf2)#绘制基因的重要性图
pdf(file="geneImportance.pdf", width=6.2, height=7.8)
varImpPlot(rf2, main="")
dev.off()#挑选疾病特征基因
rfGenes=importance[order(importance[,"MeanDecreaseGini"], decreasing = TRUE),]
rfGenes=names(rfGenes[rfGenes>1])     #挑选重要性评分大于1的基因
#rfGenes=names(rfGenes[1:30])         #挑选重要性评分最高的30个基因
write.table(rfGenes, file="rfGenes.txt", sep="\t", quote=F, col.names=F, row.names=F)#输出重要基因的表达量
sigExp=t(data[,rfGenes])
sigExpOut=rbind(ID=colnames(sigExp),sigExp)
write.table(sigExpOut, file="rfGeneExp.txt", sep="\t", quote=F, col.names=F)