欢迎关注R语言数据分析指南
❝本节分享nature communications论文中的随机森林模型分析，论文提供了该分析的数据及代码，均能正常运行。小编根据个人对数据的理解来稍作解读仅供参考，具体的内容请参考论文中的详细介绍。
论文
Prophage-encoded antibiotic resistance genes are enriched in human-impacted environments
论文数据代码
❝https://zenodo.org/records/13301199
下载数据文件夹解压缩后，RF script内即为该分析相关的数据+代码。
结果图
❝该分析通过主成分分析（PCA）对与人类活动相关的变量进行降维，根据Kaiser-Guttman准则（特征值大于1）选取了13个主要的主成分。随后，利用随机森林模型评估这些主成分对细菌抗生素抗性基因（ARGs）总丰度的影响，以量化人类活动对ARGs地理分布的影响程度。模型使用R的randomForest和rfPermute包构建，随机种子设为123，初始参数通过将70%的数据用于训练、30%用于验证进行优化。最终模型在所有数据上构建，参数设定为importance=TRUE、ntree=500、nrep=1000，并通过1000次置换检验评估模型的显著性和交叉验证的R²值。结果显示，人类活动相关的主成分对ARGs的分布具有显著影响。
数据整理
library(vegan) # 加载vegan包，用于生态学和多元统计分析
library(ggplot2) # 加载ggplot2包，用于数据可视化
library(RColorBrewer) # 加载RColorBrewer包，提供颜色调色板
library(tidyverse) # 加载tidyverse包，包含数据处理和可视化工具集
# 读取数据文件，指定分隔符为制表符（\t），首行为列名
data_pca <- read.delim("HH.bacteria_ARG_RF.txt", sep = "\t", header = TRUE)
# 划分数据集
set.seed(123) # 设置随机种子，保证结果可重复
sample_indices <- sample(1:nrow(data_pca), nrow(data_pca) * 0.3) # 从数据集中随机抽取30%的行索引作为测试集
test_data <- data_pca[sample_indices, ] # 根据索引提取测试集数据
train_data <- data_pca[-sample_indices, ] # 剩余的数据作为训练集
使用随机森林模型预测HH.ARG变量
library(randomForest) # 加载randomForest包，用于构建随机森林模型
rf_model <- randomForest(
HH.ARG ~ ., # 公式，表示以HH.ARG为因变量，其他所有变量为自变量
data = train_data, # 指定训练数据集
ntree = 500, # 设置随机森林中决策树的数量为500
importance = TRUE # 计算变量的重要性
)
rf_model # 输出模型的摘要信息
# 使用测试集进行预测
predictions <- predict(rf_model, newdata = test_data) # 基于模型对测试集数据进行预测
# 评估模型性能
mse <- mean((predictions - test_data$HH.ARG)^2) # 计算预测值与实际值之间的均方误差（MSE）
r_squared <- 1 - mse / var(test_data$HH.ARG) # 计算决定系数R²，表示模型的解释能力
print(paste("Mean Squared Error (MSE): ", mse)) # 输出MSE值
print(paste("R-squared (R^2): ", r_squared)) # 输出R²值
进行变量重要性分析
library(rfPermute) # 加载rfPermute包，用于评估变量重要性的显著性
# 使用rfPermute函数重新进行随机森林分析，并计算变量重要性的p值
set.seed(123) # 设置随机种子
factor_rfP <- rfPermute(
HH.ARG ~ ., # 公式，与之前相同
data = data_pca, # 使用完整的数据集
ntree = 500, # 设置决策树数量为500
importance = TRUE, # 计算变量的重要性
nrep = 1000, # 设置置换的重复次数为1000次，以计算p值
num.cores = 6 # 使用6个CPU核心进行并行计算，提高计算效率
)
factor_rfP # 输出rfPermute分析结果
# 提取变量的重要性得分并标准化
importance_factor.scale <- data.frame(
importance(factor_rfP, scale = TRUE), # 提取并标准化变量的重要性
check.names = FALSE # 防止R修改列名
)
importance_factor.scale # 显示变量的重要性得分
# 提取变量的重要性得分的显著性p值
importance_factor.scale.pval <- factor_rfP$pval[ , , 2] # 获取p值矩阵的第二层，即%IncMSE的p值
importance_factor.scale.pval # 显示变量的重要性p值
# 按照%IncMSE（均方误差增加的百分比）对变量进行降序排序
importance_factor.scale <- importance_factor.scale[order(importance_factor.scale$'%IncMSE', decreasing = TRUE), ]
importance_factor.scale # 显示排序后的变量重要性得分
绘制变量的重要性柱状图
# 添加变量名称列，方便绘图时使用
importance_factor.scale$OTU_name <- rownames(importance_factor.scale)
# 将变量名称设置为因子，并按照排序后的顺序
importance_factor.scale$OTU_name <- factor(importance_factor.scale$OTU_name, levels = importance_factor.scale$OTU_name)
p <- ggplot() + geom_col(
data = importance_factor.scale,
aes(x = reorder(OTU_name, `%IncMSE`), # 变量名称，按照%IncMSE排序
y = `%IncMSE`, # y轴为%IncMSE值
fill = `%IncMSE`),colour = 'black',width = 0.8) +
coord_flip() + # 翻转坐标轴，使x轴为纵轴
scale_fill_gradientn(
colors = c(
low = brewer.pal(5, "Blues")[2], # 低值对应的颜色
mid = brewer.pal(5, "YlGnBu")[1], # 中间值对应的颜色
high = brewer.pal(6, "OrRd")[4] # 高值对应的颜色
),
values = scales::rescale(c(0, 15, 50)), # 定义颜色渐变的中断点
guide = FALSE, # 不显示图例
aesthetics = "fill" # 应用于填充颜色
) +
labs(
title = NULL, # 不设置标题
x = NULL, # x轴标签为空
y = 'MSE增加百分比 (%)', # y轴标签
fill = NULL # 填充颜色图例标签为空
) +
theme_test() + # 使用内置的测试主题
theme(
panel.grid = element_blank(), # 去除网格线
panel.background = element_blank(), # 去除背景
axis.line = element_line(colour = 'black') # 设置坐标轴线为黑色
) +
theme(axis.text.x = element_text(angle = 90, hjust = 1)) + # x轴文本旋转90度，右对齐
scale_y_continuous(
expand = c(0, 0), # y轴从0开始，无额外扩展
limit = c(0, 78) # 设置y轴的范围为0到78
)
p
# 标注变量的显著性信息
for (OTU in rownames(importance_factor.scale)) {
# 获取对应变量的p值
importance_factor.scale[OTU, '%IncMSE.pval'] <- importance_factor.scale.pval[OTU, '%IncMSE']
# 根据p值设置显著性符号
if (importance_factor.scale[OTU, '%IncMSE.pval'] >= 0.05) {
importance_factor.scale[OTU, '%IncMSE.sig'] <- '' # 不显著，空白
} else if (importance_factor.scale[OTU, '%IncMSE.pval'] >= 0.01 & importance_factor.scale[OTU, '%IncMSE.pval'] < 0.05) {
importance_factor.scale[OTU, '%IncMSE.sig'] <- '*' # p < 0.05，标记为*
} else if (importance_factor.scale[OTU, '%IncMSE.pval'] >= 0.001 & importance_factor.scale[OTU, '%IncMSE.pval'] < 0.01) {
importance_factor.scale[OTU, '%IncMSE.sig'] <- '**' # p < 0.01，标记为**
} else if (importance_factor.scale[OTU, '%IncMSE.pval'] < 0.001) {
importance_factor.scale[OTU, '%IncMSE.sig'] <- '***' # p < 0.001，标记为***
}
}
# 在图中添加显著性标注和%IncMSE值
p <- p +
annotate(
"text",
x = importance_factor.scale$OTU_name, # 变量名称
y = importance_factor.scale$`%IncMSE`, # %IncMSE值
label = sprintf("%.2f%% %s", importance_factor.scale$`%IncMSE`, importance_factor.scale$`%IncMSE.sig`), # 显示数值和显著性符号
hjust = -0.2, # 水平偏移，使文本靠右一些
size = 3 # 字体大小
)
p #
使用A3包评估模型的整体p值
该步骤非常耗费运行资源
library(A3) # 加载A3包，用于模型的p值计算
# 使用a3函数计算模型的整体p值
set.seed(123) # 设置随机种子
factor_forest.pval <- a3(
HH.ARG ~ ., # 模型公式
data = data_pca, # 数据集
model.fn = randomForest, # 指定模型函数为randomForest
p.acc = 0.001, # p值计算的精度，越小计算次数越多
model.args = list(
importance = TRUE, # 计算变量的重要性
ntree = 500 # 决策树数量为500
)
)
factor_forest.pval # 输出模型p值结果
model.out <- as.data.frame(factor_forest.pval$table) # 提取结果表格
p.value <- as.numeric(gsub("<", "", model.out[1, "p value"])) # 获取模型的p值，并去除可能的"<"符号
# 在图形上添加模型的R²值和p值注释
p <- p +
annotate('text',
label = 'Bacterial ARGs', # 模型名称或标题
x = 3, # x轴位置
y = 40, # y轴位置
size = 3 # 字体大小
) +
annotate(
'text',
label = sprintf('italic(R^2) == %.2f', mean(factor_rfP[["rf"]][["rsq"]])), # 显示平均R²值
x = 2.5, # x轴位置
y = 40, # y轴位置
size = 3,
parse = TRUE # 解析数学表达式
) +
annotate(
'text',
label = sprintf('italic(P) < %.3f', p.value), # 显示模型的p值
x = 2, # x轴位置
y = 40, # y轴位置
size = 3,
parse = TRUE # 解析数学表达式
)
p
关注下方公众号下回更新不迷路
购买介绍
❝本节介绍到此结束，有需要学习R数据可视化的朋友欢迎到淘宝店铺：R语言数据分析指南，购买小编的R语言可视化文档(2024版)，购买将赠送2023年的绘图文档内容。目前此文档(2023+2024)已经更新上传200案例文档，每个案例都附有相应的数据和代码，并配有对应的注释文档，方便大家学习和参考。
2024更新的绘图内容将同时包含数据+代码+注释文档+文档清单，2023无目录仅有数据文件夹，小编只分享案例文档，不额外回答问题，无答疑服务，零基础不推荐买。
案例特点
所选案例图均属于个性化分析图表完全适用于论文发表，内容异常丰富两年累计发布案例图200+，2024年6月起提供html版注释文档更加直观易学。文档累计上千人次购买拥有良好的社群交流体验。R代码结构清晰易懂，为防止中文乱码提供单独的注释文档
R代码结构清晰易懂，为防止中文乱码2024年6月起提供单独html注释文档
