简介面向微生物网络分析中节点模块内连通度与模块间连通度的量化需求这份资源提供了基于R语言的完整计算方案适用于生态学、生物信息学等领域研究者。压缩包内共2个文件包含1个R脚本和1个graphml网络文件脚本可直接运行并输出两类连通度指标graphml文件既可作为示例网络导入分析工具也便于替换为自有数据整体体积仅18KB轻量易用。目前已有3244人学习下载适合希望掌握网络模块化指标计算方法的入门至进阶用户。通过该资料读者可获得可复用的连通度计算代码、标准图数据格式样例并理解模块内/间连通度在微生物群落功能解析中的实际应用逻辑为后续开展网络关键节点识别、群落稳定性评估等研究提供工具与思路参考。1. 从微生物网络到关键物种识别这份资源到底在算什么微生物网络里连边最多的节点不一定是让整个群落在扰动后还能撑住的节点。真正决定稳定性的往往是那些在模块内部连接紧密、又在模块间充当桥的物种。这正是模块内连通度Zi和模块间连通度Pi想回答的问题。这份资源给出一套可直接复现的R脚本zi_pi_test.R和一个GraphML格式网络文件mygraph.graphml把网络导入 igraph、用社区检测切好模块、逐节点统计模块内与模块间连接情况最终得到每个网络节点的 Zi、Pi 取值和角色分类。适合已经在做 OTU/ASV 共现网络、手里有网络文件和 R 环境、但卡在“找关键物种”这一步的研究者。不需要从零推公式把输入输出路径一改就能跑。2. 模块内连通度与模块间连通度公式拆解和算法选型2.1 为什么微生物网络要同时看模块内和模块间两个维度对于微生物网络节点是 OTU/ASV 或物种边是共现或互斥关系。最常见的错误是用度中心性直接排序然后宣布一堆高连接物种是关键节点。但度中心性只看连了几条边不看边分布在哪一个腐殖质降解菌连到十几个不同模块的不同物种度很高但对其中任意一个模块内部的互助结构参与得都很浅另一个菌只在自己的模块内部和十来个个邻居打成一片度不算高却可能是维持该功能的承重墙。Zi-Pi 两个指标把“连接数量”拆成“模块内贡献”和“模块间跨接”两个坐标能避免这类误判。在生态网络的经典工作里这两个坐标被用来做角色分类模块内连通度高的节点被称为模块中心Module hub模块间连通度高的节点被称为连接者Connector两者都高的则是网络中心Network hub都低的属于外围节点Peripheral。在微生物群落里这种分类有直观的生物学含义高 Zi 往往对应同一生态功能组内部的活跃成员高 Pi 往往对应参与多个生物地球化学过程的跨界物种。因此 Zi-Pi 分析的结果常被当作“候选关键物种”名单。2.2 公式拆解什么在影响 Zi 和 Pi这两个量在文献里最常见的代码命名就是Zi和Pi后面脚本里的变量也沿用这个约定。先说模块内连通度。用 k_i 表示节点 i 的总度数κ_i 表示 i 与其所属模块内部节点的连边数。设节点 i 所属模块记为 S模块内所有节点的 κ 值算均值为 μ_S、标准差为 σ_S那么模块内连通度定义为z_i (κ_i - μ_S) / σ_S这个式子本质是把“该节点在模块内部的连接数”相对同模块平均水平做一个标准化。z_i 大于 2.5相当于该节点在模块内部的活跃程度显著高于同模块同伴z_i 为负则说明它偏居模块边缘。需要特别说明的是当 σ_S0 时公式会除零脚本里会强制取 0否则整列全是 NaN。模块间连通度用的是参与系数participation coefficient的框架。记 κ_is 为节点 i 与模块 s 之间的连边数s 遍历网络中全部模块那么p_i 1 - Σ (κ_is / k_i)^2这个式子的观察点很有意思。一个节点如果把所有边都连在自身模块内那么它到其他模块的 κ_is 全是 0只有自身模块那一项为 1p_i0如果它的边平摊到许多模块上平方和会变小p_i 会趋近 1。所以 Pi 衡量的是“跨模块连接在总连接中的占比”而不是“跨模块连接的总条数”。一个只有 4 条边的边缘节点可能 Pi 高达 0.75而一个有 50 条边但 49 条都连在自身模块内的中枢节点 Pi 只有 0.04。这两个值放在一起才能把“局部核心”和“全局桥梁”两种角色分开。指标代码变量取值范围微生物网络中的解释模块内连通度 ZiZi无固定上下界常见 -35高值同模块内互动远超平均水平疑似模块功能核心低值模块边缘或不稳定成员模块间连通度 PiPi01高值连接对象分散在多个模块疑似跨功能桥梁低值连接几乎全部局限在自身模块内2.3 模块划分算法怎么选Louvain、Infomap 与 ModuLandZi-Pi 的可靠性完全建立在模块划分的基础上。zi_pi_test.R 默认走 igraph 的cluster_louvain路径原因有三参数少对加权网络有原生支持几千到几万节点的微生物网络在普通笔记本上都能跑完。Louvain 优化的是模块度 Q倾向于让模块内部边尽可能密、模块间边尽可能稀这个目标正好与“共现网络找生态功能组”的需求一致。Infomap 用随机游走信息编码做社区发现对边权重比例比 Louvain 更敏感带权、带方向的网络下结果常有差异。如果网络是时间序列挖掘出来的且权重跨度大可以换 Infomap 验证。ModuLand 系列方法支持重叠模块能处理一个节点同时属于多个功能模块的情形这个特性很吸引做多功能微生物的研究者但代价是输出结构复杂计算 Zi-Pi 前还得自行决定重叠归属怎么折算。常见做法是把重叠模块合并为唯一归属或者一个节点产出多组 Zi-Pi 再取均值。我的建议是默认用 Louvain 出主结果用 Infomap 和重叠模块方法做敏感性验证而不把宝押在单一的算法上。3. 把 Zi-Pi 跑出来读取 GraphML、计算指标并绘制散点图3.1 准备工作用 igraph 读取 mygraph.graphml 并核对属性在跑任何计算前先把网络读取进来并检查属性名。GraphML 是比较通用的网络交换格式mygraph.graphml里通常同时带节点属性和边属性但不同工具导出的 GraphML 属性命名习惯不同有的边权重叫weight有的叫value、correlation有的根本不带权重。直接用E(g)$weight取值之前先列出全部边属性名这一步能省掉后面大半的报错。library(igraph) g - read_graph(mygraph.graphml, format graphml) cat(节点数:, vcount(g), 边数:, ecount(g), \n) print(edge_attr_names(g)) if (weight %in% edge_attr_names(g)) { cat(权重范围:, range(E(g)$weight), \n) } else { warning(没有找到weight属性将按无权网络计算) }GraphML 格式本身有命名空间read_graph不指定format也能自动识别但显式写明format graphml更稳妥。edge_attr_names的返回结果决定后面计算按有权还是无权处理。如果网络中含有负权重优先考虑把它清洗掉或转成正相关强度再做模块划分否则模块结构会被明显拉歪这一点第 4 章会细说。到这一步数据读取和属性核验就完成了。3.2 核心环节模块识别与 Zi-Pi 计算读取网络后先用 Louvain 方法做模块划分然后把模块结果写回节点属性逐节点统计模块内邻居数、总度数和到各模块的连接数最后套用上一章的公式。下面这段是按zi_pi_test.R的思路组织的最简版本输出一张包含节点名、所属模块、模块内度、总度、Zi、Pi 的表。g - induced_subgraph(g, V(g)[degree(g) 0]) set.seed(123) com - cluster_louvain(g, weights E(g)$weight) V(g)$module - membership(com) intra_deg - sapply(V(g), function(v) { nbr - neighbors(g, v) sum(V(g)$module[nbr] V(g)$module[v]) }) total_deg - degree(g, loops FALSE) deg_to_module - sapply(V(g), function(v) { nbr - neighbors(g, v) tabulate(factor(V(g)$module[nbr], levels seq_len(max(V(g)$module))), nbins max(V(g)$module)) }) p_i - 1 - rowSums((deg_to_module / total_deg)^2) z_i - numeric(vcount(g)) for (i in seq_len(vcount(g))) { mod - V(g)$module[i] idx - which(V(g)$module mod) mu - mean(intra_deg[idx]) sdv - sd(intra_deg[idx]) z_i[i] - ifelse(sdv 0, 0, (intra_deg[i] - mu) / sdv) } result - data.frame(node V(g)$name, module V(g)$module, intra_deg intra_deg, total_deg total_deg, Zi round(z_i, 4), Pi round(p_i, 4)) write.csv(result, zi_pi_output.csv, row.names FALSE)先解释过滤那行induced_subgraph(g, V(g)[degree(g) 0])把孤立节点剔掉。度数为 0 的节点不参与任何模块内部统计但在p_i计算里会让total_deg变成 0直接除零出 NaN所以先过滤最省心。这样会改变节点数写报告时要注明“剔除孤立节点 N 个”。逐段看逻辑。set.seed(123)用来锁住 Louvain 的随机初始化cluster_louvain的weights参数接收边权重如果第 3.1 步发现没有weight属性这里参数要改成weights NULL。intra_deg的sapply里neighbors(g, v)返回节点 v 的所有邻居顶点 idV(g)$module[nbr]取出邻居的模块编号与 v 自身模块比较后求和得到模块内连接边数。total_deg用degree(g, loops FALSE)把自环排除掉因为微生物共现网络中一个节点到自己的边通常没有语义。deg_to_module里用tabulate统计 v 到每个模块的连接数得到公式里的 κ_is 矩阵。p_i的计算里以deg_to_module / total_deg得到每个模块连接占总连接的比例平方后按行求和再用 1 减正好是参与系数。z_i的循环里最关键的是sdv 0的判断单节点模块或者模块内所有节点模块内度都相同的情况下方差为 0此时直接让 Zi 取 0避免产出 NaN。参数上需要根据实际情况调整的地方有两个。第一max(V(g)$module)必须在tabulate前拿到因为 Louvain 返回的模块编号不一定是连续的如果硬编码模块数量遇到空编号会报错或错位。第二Louvain 的分辨率可以直接传resolution参数默认 1.0网络明显被切得过碎或过粗时再调不要一上来就动。跑完把zi_pi_output.csv打开扫一眼如果 Pi 全为 0说明整张网络被划分成了一个模块需要回过去看网络的连通性或者分辨率设置。注意如果结果里 Pi 全为 0先检查最大模块数量是否等于 1Louvain 把全网划成一个模块时Pi 失去区分度需要从分辨率或网络连通性入手排查。3.3 可视化Zi-Pi 散点图与生态角色着色计算完成后把 Zi 放纵轴、Pi 放横轴画散点图。这里用经典阈值 2.5 和 0.62 把节点划成四类并用 ggrepel 给非外围节点加标签避免全图标签堆成黑疙瘩。library(ggplot2) library(ggrepel) result$role - ifelse(result$Zi 2.5 result$Pi 0.62, Network hub, ifelse(result$Zi 2.5, Module hub, ifelse(result$Pi 0.62, Connector, Peripheral))) ggplot(result, aes(Pi, Zi, color role)) geom_point(size 2.2, alpha 0.75) geom_hline(yintercept 2.5, linetype dashed) geom_vline(xintercept 0.62, linetype dashed) geom_text_repel(aes(label ifelse(role ! Peripheral, node, )), max.overlaps 15) labs(x 模块间连通度 Pi, y 模块内连通度 Zi, color 生态角色) theme_minimal() ggsave(zi_pi_plot.png, width 8, height 6, dpi 300)画图的逻辑对应生态网络的经典分类两者都高的是 Network hub全局和局部都很活跃Zi 高 Pi 低是 Module hub模块内部核心Zi 低 Pi 高是 Connector跨模块桥梁两者都低就是 Peripheral 外围节点。linetype dashed两条阈值线只是参照不是坐标轴。geom_text_repel只标非外围节点max.overlaps 15限制最多展示 15 个标签节点多时省得图糊成一团。dpi 300保证放进论文或汇报里不虚。4. 避坑指南Zi-Pi 计算中五个容易翻车的地方4.1 现象Pi 算出来几乎全大于 0.9所有节点都成了 Connector这是模块划分太碎的表现。Louvain 在分辨率默认 1.0 时倾向于把网络切得比较细尤其是大型稀疏共现网络模块数量一多节点的连接只要稍微分散就很容易让 Pi 冲高。具体到结果上x 轴的 0.62 阈值线失去区分意义。原因模块数量过多导致 κ_is 的分布被摊薄平方和缩小Pi 整体抬高。解决先跑一下table(V(g)$module)看模块大小分布如果大量模块只有两三个节点说明划分过碎。常见做法是降低分辨率例如cluster_louvain(g, resolution 0.6)让模块合并得更粗或者直接用 Infomap 算法交叉验证看 Pi 分布是否回归到合理的 0 到 1 散布范围。注意改分辨率后要重新跑一遍计算不能只改画图阈值。4.2 现象Zi 输出列出现成片 NaN出现 NaN 的原因基本集中在一个模块内的 κ 值没有波动标准差为 0。单节点模块、只有两条自环的模块或者模块内所有节点模块内度恰好相等的极端数据都会触发除零。原因z_i (κ_i - μ_S) / σ_S在 σ_S0 时未做保护除零后 R 返回 NaN。解决计算前先过滤孤立节点g - induced_subgraph(g, V(g)[degree(g) 0])计算中在sdv 0时让 Zi 取 0而不是保留 Inf 或 NaN。上面第 3.2 节的实现里已经写了ifelse(sdv 0, 0, ...)但如果数据里有特别怪异的模块结构仍建议在写 CSV 之前对 NaN 做一次显式检查result[is.na(result$Zi), ]看到底是哪些节点出了问题。4.3 现象同一份网络两次结果里的关键节点对不上Louvain 是启发式优化算法初始化顺序和随机种子都会影响模块划分结果。微生物网络动辄上千节点不同次运行可能给出不同的模块边界进而让同一节点的 Zi/Pi 发生漂移。原因没有固定随机种子或者模块本身在真实结构中就不稳定算法在不同局部最优解之间切换。解决set.seed()在运行前固定更稳妥的做法是多次运行并只把模块归属一致性高的节点纳入最终关键节点名单。对每次运行都计算一次 Zi/Pi然后看不同种子之间排名稳定程度如果某个节点在三四次运行里一会儿是 Module hub 一会儿是 Peripheral那它就不该被写进结论。报告里可以给“模块归属不确定”单独留一列备注不硬下结论。4.4 现象导入 GraphML 后脚本报“没有 weight 属性”或节点名称丢失GraphML 是一个容器格式不同工具写出来的属性键名差异很大。Cytoscape 或 Gephi 导出的文件里边权重可能叫value、叫co-expression节点名称可能不在name属性里。原因属性键名不标准脚本写死了取E(g)$weight和V(g)$name遇到不匹配的键名就失效。解决先用edge_attr_names(g)和vertex_attr_names(g)列出所有属性名再决定映射。如果原有键名不是weight常见做法是E(g)$weight - E(g)$value或者用set_edge_attr显式改名。节点名同理V(g)$name不存在时用expr或label列手工赋给name。属性核对应当作为流程第一步而不是报错后再来。4.5 现象把负相关权重直接喂给 Louvain模块结构全错微生物共现网络经常保留正相关与负相关构建网络时 SPIEC-EASI、SparCC 的输出可能既有正权重又有负权重。直接把这些权重作为 Louvain 的weights传入等于把强负相关当成强吸引模块会被拉向完全错误的边界。原因Louvain 的模块度优化隐含“权重高亲缘近”的语义负权重不满足这个前提会把负相关配对强行聚成一个模块。解决计算 Zi-Pi 前先把边表清洗到正权重语义上。常见做法是只保留正相关边把原边表过滤到相关系数大于 0 的边后再建图如果非要保留负相关结构就把负相关单独建一张负网络分开分析不要和正网络混在一个模块划分里跑。这个决策要在数据清洗阶段定好并写进方法部分因为评审对负相关边如何处理非常敏感。5. 进阶用法阈值怎么改、模块怎么动态观察、结果怎么用5.1 阈值 2.5 和 0.62 不是铁律如何根据网络特性调整2.5 这个值近似正态分布约 99% 置信区间0.62 则是早期生态网络研究里参与系数分布的观察阈值它们并非从微生物网络本身推导出来的。实际项目里网络规模、模块数量会严重影响这两个阈值的位置。比如一个只有 200 个节点的共现网络Zi 很少能超过 2.5会漏掉不少实际重要的模块成员。常见做法是用分位数替换固定阈值取全网络 Zi 和 Pi 的 95% 或 99% 分位数作为门槛再重新划分角色。qz - quantile(result$Zi, 0.95, na.rm TRUE) qp - quantile(result$Pi, 0.95, na.rm TRUE) result$role_q - ifelse(result$Zi qz result$Pi qp, Network hub, ifelse(result$Zi qz, Module hub, ifelse(result$Pi qp, Connector, Peripheral)))用分位数阈值的优点是适应不同尺度的网络缺点是没有统一参照跨研究比较时不方便。我一般的习惯是正文用经典 2.5/0.62 画主图补充材料里给一版按分位数划分的结果两者一致的地方才是真正值得写进结论的关键物种。阈值调整时还要记得把平均路径长度和聚类系数一起算出来看因为这些网络级指标可以反映模块划分的整体质量如果聚类系数在阈值调整后显著下降说明模块切得不自然。5.2 动态模块化观察改分辨率追踪单个节点的角色漂移项目摘要里提到可以通过改变模块划分的分辨率或动态模块化分析观察连通度随网络结构的变化。这个思路在实操里就是循环遍历resolution对同一个网络做多次 Louvain 划分然后把某个关注节点的 Zi/Pi 轨迹画出来。如果一个节点在粗尺度下是 Module hub、在细尺度下变成 Connector那它的生态角色就依赖于观察粒度这样的节点需要回到功能注释或实验里找证据不能直接靠一次划分下结论。res_range - seq(0.5, 2.0, by 0.1) track_list - lapply(res_range, function(res) { set.seed(42) com2 - cluster_louvain(g, resolution res) V(g)$module2 - membership(com2) intra2 - sapply(V(g), function(v) { nbr - neighbors(g, v) sum(V(g)$module2[nbr] V(g)$module2[v]) }) deg2 - degree(g, loops FALSE) deg_to_mod2 - sapply(V(g), function(v) { nbr - neighbors(g, v) tabulate(factor(V(g)$module2[nbr], levels seq_len(max(V(g)$module2))), nbins max(V(g)$module2)) }) p2 - 1 - rowSums((deg_to_mod2 / deg2)^2) z2 - numeric(vcount(g)) for (i in seq_len(vcount(g))) { mod - V(g)$module2[i] idx - which(V(g)$module2 mod) sdv - sd(intra2[idx]) z2[i] - ifelse(sdv 0, 0, (intra2[i] - mean(intra2[idx])) / sdv) } data.frame(resolution res, node V(g)$name, Zi z2, Pi p2) }) track_df - do.call(rbind, track_list)resolution大于 1 会让模块切得更细小于 1 会让模块合并得更粗。track_df输出后挑出那些角色变化最大的节点单独画线图就能直观看到“粒度敏感性”。这类动态观察结论常见于生态网络的方法学讨论但在微生物组里很少有人展示写进补充材料会显得分析很扎实。注意这里计算 Pi 时如果deg2为 0 会除零建议循环前把孤点过滤掉这一步和第 3.2 节一样。5.3 把 Zi-Pi 的结果接到鲁棒性分析和功能注释里Zi-Pi 算出来不只是为了画一张散点图。比较常见的下游用法是把不同角色节点作为移除对象做网络扰动模拟观察最大连通片的变化。若先移除所有 Connector、再移除所有 Module hub看哪个组对网络连通性的打击更大就能判断跨模块节点是否真的是网络的承重结构。robustness_curve - function(g, remove_by) { g2 - g max_comp - numeric(length(remove_by)) for (i in seq_along(remove_by)) { g2 - delete_vertices(g2, remove_by[i]) if (vcount(g2) 0) break comp - components(g2) max_comp[i] - max(comp$csize) / vcount(g2) } max_comp } connectors - result$node[result$role Connector] hubs - result$node[result$role Module hub] curve_c - robustness_curve(g, connectors) curve_h - robustness_curve(g, hubs)这条曲线可以在鲁棒性分析里画出两条下降曲线比较哪个角色被删后最大连通片下降更快。再把高 Zi、高 Pi 节点的功能注释KEGG 通路或 COG 类别列出来看是否与已知功能模块对应。Zi-Pi 在这里的角色更像“筛选器”负责把候选节点压到小名单功能验证仍然要靠注释或实验。我一般会在最终报告里同时呈现角色表、扰动曲线和通路富集三样东西因为单独一张散点图的生态学说服力有限。6. 验证你的 Zi-Pi 结果三种交叉检验方法6.1 人工构造网络做正对照用一个已知结构的图做正向校验是最快的排错方式。比如构造三个全连接模块模块间只留几条桥边跑完脚本预期桥节点的 Pi 应当显著偏高而模块内部全连接节点的 Zi 才会高。如果这个人工例子的输出不符合预期问题多半在数据处理环节。tg - make_full_graph(10) make_full_graph(10) make_full_graph(10) tg - tg edges(c(1, 15, 2, 23, 5, 19)) set.seed(123) com_t - cluster_louvain(tg)这个例子虽然简单却能立刻暴露tabulate长度、除零这类基础问题。正控通过后再跑真实网络才有意义。6.2 用第二、第三套算法交叉验证Louvain 的结果和 Infomap、标签传播算法给出的 Zi/Pi 做相关性分析。保留在至少两套算法下都落入同一角色的节点算法之间互相矛盾的节点标记为不稳定节点不放进关键物种讨论范围。com_infomap - cluster_infomap(g, e.weights E(g)$weight)Infomap 的输入参数是e.weights注意和 Louvain 的weights参数名不同。交叉验证的核心目的不是判断哪个算法更好而是找交集。生态学结论应该建立在交集之上而不是一次 Louvain 的随机结果。6.3 把计算结果映射回原始网络把高 Pi 节点在原图中高亮逐一检查它是否真的连接了多个子图区域把高 Zi 节点检查是否为所在模块的连接稠密区域成员。这个人工目检虽然主观但能拦住不少因为模块划分边界错误造成的假阳性。再看平均路径长度和聚类系数是否落在合理区间如果聚类系数异常低说明模块边界取得太碎优先回去调分辨率而不是继续解读结果。从那以后我每次跑 Zi-Pi 都强制走一遍正控、交叉算法、映射回原图三步尤其是报告要给合作方之前宁可多花半天在这三步上也不愿被数据和结论不一致坑一次。希望帮到你。本文还有配套的精品资源点击获取
