扩增子分析|microeco包妙用之网络鲁棒性、易损性、内聚力和网络稳定性一行代码获取结果（meconetcomp包）

1、引言

meconetcomp包是microeco包的扩展包。两个包均由福建农林大学姚敏杰教授团队开发，专注于扩增子测序数据的集成分析。Meconetcomp发表在国人创建的顶级生物信息学期刊iMeta中。其中meconetcomp包提供了一种比较微生物共现网络的流程，具有较高的灵活性和可扩展性，可以帮助用户高效地比较不同组样本或不同构建方法构建的网络。（一行代码可实现某个特性分析），因此，掌握了这套代码，可让你快速通关网络分析。本节我们提供了微生物网络分析之网络稳定性评价指标鲁棒性、易损性和内聚力的代码。

⚠️值得注意的是，在周集中老师团队采用的网络相关性计算方法充分考虑了稀有物种，因此计算所耗时间较长。本套代码未使用该方法运行较快。此外，该套代码无论是网络指标计算还是出图均是一行代码可实现，大大简化了操作。不过对于一些个性化的分析要求，因人而异，个人建议可通过github原始代码进行适当修改后使用。

2、示例数据和代码

下述代码和上一期代码可以一起使用也可以分开使用，看个人喜好，上一期代码是：

扩增子分析|microeco包妙用之网络属性（拓扑、边、节点）及多网络比较（子网络拓扑）一行代码获取结果（meconetcomp包）

2.1、数据准备

🌟准备数据

⚠️输入了四个数据分别是otu表（feature_table.csv），样本信息表（sample_table.csv），物种分类表（tax_table.csv）和系统发育树文件（phylo_tree.tre）。

⚠️一般数据输入采用上述4个表格会流畅的获得相关结果，但是由于数据来源的差异，部分数据表格可能没有，因此经测试最少需要otu表（feature_table.csv），样本信息表（sample_table.csv）可满足大部分输出结果。

（1）otu表（feature_table.csv）

（2）样本信息表（sample_table.csv）

（3）物种分类表（tax_table.csv）

（4）系统发育树文件（phylo_tree.tre）

2.2、R包准备和数据读取以及微生物共现性网络构建

🌟基础操作之meconetcomp包等依赖包安装、数据微表创建、共现性网络构建
#############包安装和加载###################
# 安装所需的 R 包
# aplot: microeco 包中 trans_venn 类的依赖包
# agricolae: 用于 Duncan 的新多重范围检验
packages <- c("meconetcomp", "rgexf", "pheatmap", "aplot", "agricolae")
# 检查每个包是否已安装，如果未安装则安装
for(x in packages) {
  if(!require(x, character.only = TRUE)) {
    install.packages(x, dependencies = TRUE) # 自动安装包及其所有依赖项
  }
}
#安装并加载ape包
if(!require("ape")){install.packages("ape")}
library(ape)
##################数据基本操作#######################
# 加载示例数据；16S rRNA基因扩增子测序数据
#读取样本分组信息表
sample_info_16S <- read.csv("sample_table.csv", row.names = 1)
#读取特征表
otu_table_16S <- read.csv("feature_table.csv",row.names = 1)
#读取系统发育树
phylo_tree_16S <- read.tree("phylo_tree.tre", tree.names = NULL)
#读取物种分类表
taxonomy_table_16S <- read.csv("tax_table.csv",row.names = 1)
# 让我们创建需要信息的 microtable 对象
soil_amp <- microtable$new(sample_table = sample_info_16S, otu_table = otu_table_16S, tax_table = taxonomy_table_16S, phylo_tree = phylo_tree_16S)
#########创建三组网络############
# 选择"IW"组的样本
# 使用clone对soil_amp对象深度拷贝
tmp <- clone(soil_amp)
# 直接更改样本表
tmp$sample_table %<>% subset(Group == "IW")
# 清理对象中的所有文件
tmp$tidy_dataset()
# 使用filter_thres参数过滤低丰度特征
tmp <- trans_network$new(dataset = tmp, cor_method = "spearman", filter_thres = 0.0005)
# 设置 p 值阈值
# 设置相关系数阈值
tmp$cal_network(COR_p_thres = 0.01, COR_cut = 0.6)
# 将网络加入列表中
soil_amp_network$IW <- tmp
# 选择"TW"组的样本
tmp <- clone(soil_amp)
tmp$sample_table %<>% subset(Group == "TW")
tmp$tidy_dataset()
tmp <- trans_network$new(dataset = tmp, cor_method = "spearman", filter_thres = 0.0005)
tmp$cal_network(COR_p_thres = 0.01, COR_cut = 0.6)
soil_amp_network$TW <- tmp
# 选择"CW"组的样本
tmp <- clone(soil_amp)
tmp$sample_table %<>% subset(Group == "CW")
tmp$tidy_dataset()
tmp <- trans_network$new(dataset = tmp, cor_method = "spearman", filter_thres = 0.0005)
tmp$cal_network(COR_p_thres = 0.01, COR_cut = 0.6)
soil_amp_network$CW <- tmp
# 现在已创建网络列表 soil_amp_network

2.3、网络鲁棒性指标计算

🌟示例代码

################## 网络鲁棒性 ##################
# 创建一个 'robustness' 对象来评估网络的鲁棒性，即网络在不同移除策略下的抗扰动能力
# remove_strategy 参数指定移除策略，包括随机移除边("edge_rand")、移除最强的边("edge_strong")、随机移除节点("node_rand")和移除度数高的节点("node_degree_high")
# remove_ratio 指定移除比例，从0到0.99，步长为0.1
# measure 参数指定鲁棒性指标，Eff 表示效率，Eigen 表示特征值，Pcr 表示聚类系数
# run = 10 表示每种配置运行10次，便于计算平均鲁棒性
tmp <- robustness$new(soil_amp_network, remove_strategy = c("edge_rand", "edge_strong", "node_rand", "node_degree_high"), 
                      remove_ratio = seq(0, 0.99, 0.1), measure = c("Eff", "Eigen", "Pcr"), run = 10)
# 查看鲁棒性分析的结果表，包括每次实验的详细输出
View(tmp$res_table)
write.csv(tmp$res_table,"robustness_detail.csv")
# 查看鲁棒性分析的总结，包括每种策略的鲁棒性统计结果
View(tmp$res_summary)
write.csv(tmp$res_summary,"robustness_summary.csv")
# 绘制鲁棒性分析图表，显示不同移除策略下网络鲁棒性随移除比例的变化
g1<- tmp$plot(linewidth = 1)
g1
ggsave("robustness.pdf", g1, width = 10, height = 7)

🌟输出结果

robustness_detail.csv为鲁棒性详情输出结果

robustness_summary.csv是鲁棒性汇总输出结果

robustness.pdf是鲁棒性输出图纸

2.4、网络易损性指标计算

🌟示例代码

################节点易损性 #########################
# 使用 'vulnerability' 函数计算每个节点的易损性，即节点被移除时对网络整体稳定性的影响
# 易损性指标有助于识别关键节点
vul_table <- vulnerability(soil_amp_network)
View(vul_table)  # 查看易损性结果表，包含每个节点的易损性评分
#输出所有的节点的易损性值
write.csv(vul_table,"vulnerability_all_otu.csv")

🌟输出结果：

vulnerability_all_otu.csv是所有节点的易损性指标，可全部查看。

# 加载 dplyr 包
library(dplyr)
# 假设 vul_table 是你的数据框
result_table <- vul_table %>%
  group_by(Network) %>%           # 按 network 分组
  slice_max(vulnerability, n = 1) %>%  # 获取 vulnerability 最大值的行
  ungroup()                       # 取消分组
# 输出易损性值
write.csv(result_table,"max_vulnerability.csv")

🌟输出结果：

max_vulnerability.csv是所有网络易损性指标，可用于后续绘图

2.5、网络内聚力及稳定性

🌟网络内聚力计算

############### 内聚力 ####################
# 创建 'cohesionclass' 对象来计算网络的凝聚力，即网络中节点之间的紧密程度
t1 <- cohesionclass$new(soil_amp_network)
# 查看凝聚力计算的结果，包括样本级别的凝聚力信息
View(t1$res_list$sample)
#输出内聚力计算结果
write.csv(t1$res_list$sample,"Cohesion.csv")
# 查看连通性信息，例如基于各特征的连接性分析结果
View(t1$res_list$feature)
#输出连通性计算结果
write.csv(t1$res_list$feature,"连通性.csv")
# 计算不同处理组间凝聚力的差异，使用 ANOVA 检验
t1$cal_diff(method = "anova")
# 绘制内聚力图表，使用正连通性（r_pos）作为指标，显示各组别的凝聚力水平
g1<-t1$plot(measure = "r_pos")
g1
ggsave("r_pos.pdf", g1, width = 7, height = 7)

🌟输出结果

Cohesion.csv为内聚力计算结果。

连通性.csv为连通性计算结果。

图纸为正连通性结果。

🌟基于正负内聚力的网络稳定性计算

#基于正负内聚力的网络稳定性计算
net<- t1$res_list$sample
library(dplyr)
# 创建新列，值为 c_neg 的绝对值除以 c_pos
net <- net %>%
  mutate(stability = abs(c_neg) / c_pos)  # 计算并添加新列
write.csv(net,"network stability.csv")

🌟输出结果

network stability.csv表格包括网络稳定性计算结果。

3、参考文献

[1] Liu, Chi, Chaonan Li, Yanqiong Jiang, Raymond J. Zeng, Minjie Yao, and Xiangzhen Li. 2023. “ A guide for comparing microbial co-occurrence networks.” iMeta 2, e71.

[2] Chi Liu, Yaoming Cui, Xiangzhen Li, Minjie Yao, microeco: an R package for data mining in microbial community ecology, FEMS Microbiology Ecology, Volume 97, Issue 2, February 2021, fiaa255

[2] microeco_tutorial/meconetcomp-package.html#network-modularity-for-all-networks

4、联系作者

！！！有需要的小伙伴可以去评论区自取今天的测试代码和实例数据。

📌示例代码中提供了数据和代码，已经小编测试，可直接运行。

⚠️本账号会不时更新代码，代码失效可联系小编~

以上就是利用R语言microeco包和meconetcomp包解析网络鲁棒性、易损性、内聚力、网络稳定性的全部内容。

💬有任何问题可以后台联系，看文献，学方法，期待大家一起成长！