网络模块内连通度（Zi）和模块间连通度（Pi）及在R中计算

模块内连通度（Zi）和模块间连通度（Pi）及在R中计算

续前文“网络图的凝聚性特征”和“节点和边特征”，本篇简介两种特殊的节点特征，模块内连通度（Zi）和模块间连通度（Pi）。 

模块内连通度（Zi）和模块间连通度（Pi）的定义

有关网络模块的定义可参考前文“网络图的凝聚性特征”，在网络模块的基础上衍生出了两个重要的节点特征，模块内连通度（Within-module connectivity，Zi）和模块间连通度（Among-module connectivity，Pi），计算如下（Guimerà and Nunes 2005）：

对于Zi，K_i是节点i在模块S_i中与其它节点的边数量，`K_Si是模块S_i中所有节点的K值的平均值（K值为该节点在S_i中与其它节点的连接数），σ_KSi是模块S_i中所有节点的K值的标准差。

对于Pi，K_is是节点i与模块S中节点的边数量，k_i是节点i的度，M代表模块，N_M即代表了所有的模块。根据公式可知，如果与某节点有关的边在所有模块中均匀分布，则该节点的Pi值接近于1；如果与某节点有关的所有的边都在其所属的模块内，则该节点的Pi值为0。

通过模块内连通度（Zi）和模块间连通度（Pi）识别关键节点

依据节点的的拓扑特征可将节点属性分为4种类型，包括：Module hubs（模块中心点，在模块内部具有高连通度的节点，Zi > 2.5且Pi < 0.62），Connectors（连接节点，在两个模块之间具有高连通度的节点，Zi < 2.5且Pi > 0.62），Network hubs（网路中心点，在整个网络中具有高连通度的节点，Zi > 2.5且Pi > 0.62）以及Peripherals（外围节点，在模块内部和模块之间均不具有高连通度的节点，Zi < 2.5且Pi <0.62）。通常将除Peripherals外的其余3种类型的节点归为关键节点（Deng et al. 2012），归因于这些节点处于枢纽位置，不难看出这些关键节点的缺失可能会引起模块和网络的分解，它们在维持网络结构的稳定性中发挥重要作用。

因此，基于模块内连通度（Zi）和模块间连通度（Pi）的网络核心节点判别方法得到了广泛的应用。例如，微生物共发生网络一般可以被划分成多个模块，模块是网络中高度连接的区域，模块可能反应了栖息地的异质性、系统发育上亲缘关系较近物种的聚集、生态位的重叠和物种的共进化，被认为是系统发育、进化或功能上独立的单元（Olesen et al. 2007）。在生态网络模块中识别的关键节点，往往代表了在维持微生物群落结构稳定性上可能起重要作用关键物种。Shi等（2016）在野燕麦根际微生物共发生网络中通过模块内连通度（Zi）和模块间连通度（Pi）的概念寻找核心微生物物种就是依据此原理。

R语言计算通过模块内连通度（Zi）和模块间连通度（Pi）

如果提到哪个R包可以计算这两种节点特征，我还真不知道……其实根据上面的公式自己写函数计算就好了，公式理解起来很简单的。在这里就分享一个实验室的祖传代码（流传了约4年了吧），大神师兄写的，我一直在用（师兄肯定不会想到，他造福了多少向我这样的懒人）。

示例数据和R代码链接（提取码：dygt）：

https://pan.baidu.com/s/1X-C5yuYqM_GLXmJJ0oZpgA 

要计算这两种节点特征，首先需要计算网络各节点的度，并对网络划分模块。

节点度的概念及计算可见“节点和边特征”。

模块的概念及划分方法可见“网络图的凝聚性特征”。

##igraph 包计算网络模块
library(igraph)
 
#输入数据示例，邻接矩阵
#这是一个微生物互作网络，数值“1”表示微生物 OTU 之间存在互作，“0”表示无互作
adjacency_unweight <- read.delim('adjacency_unweight.txt', row.names = 1, sep = '\t', check.names = FALSE)
head(adjacency_unweight)[1:6]    #邻接矩阵类型的网络文件
 
#邻接矩阵 -> igraph 的邻接列表，获得非含权的无向网络
igraph <- graph_from_adjacency_matrix(as.matrix(adjacency_unweight), mode = 'undirected', weighted = NULL, diag = FALSE)
igraph    #igraph 的邻接列表
 
#计算节点度
V(igraph)$degree <- degree(igraph)
 
#模块划分，详情 ?cluster_fast_greedy，有多种模型
set.seed(123)
V(igraph)$modularity <- membership(cluster_fast_greedy(igraph))
 
#输出各节点（微生物 OTU）名称、节点度、及其所划分的模块的列表
nodes_list <- data.frame(
    nodes_id = V(igraph)$name, 
        degree = V(igraph)$degree, 
        modularity = V(igraph)$modularity
)
head(nodes_list)    #节点列表，包含节点名称、节点度、及其所划分的模块
 
write.table(nodes_list, 'nodes_list.txt', sep = '\t', row.names = FALSE, quote = FALSE)

然后计算模块内连通度（Zi）和模块间连通度（Pi），计算函数我放在网盘中了。

##计算模块内连通度（Zi）和模块间连通度（Pi）
source('zi_pi.r')

#上述的邻接矩阵类型的网络文件
adjacency_unweight <- read.delim('adjacency_unweight.txt', row.names = 1, sep = '\t', check.names = FALSE)

#节点属性列表，包含节点所划分的模块
nodes_list <- read.delim('nodes_list.txt', row.names = 1, sep = '\t', check.names = FALSE)

#两个文件的节点顺序要一致
nodes_list <- nodes_list[rownames(adjacency_unweight), ]

#计算模块内连通度（Zi）和模块间连通度（Pi）
#指定邻接矩阵、节点列表、节点列表中节点度和模块度的列名称
zi_pi <- zi.pi(nodes_list, adjacency_unweight, degree = 'degree', modularity_class = 'modularity')
head(zi_pi)

write.table(zi_pi, 'zi_pi_result.txt', sep = '\t', row.names = FALSE, quote = FALSE)
##可再根据阈值对节点划分为 4 种类型，并作图展示其分布
library(ggplot2)

zi_pi <- na.omit(zi_pi) #NA 值最好去掉，不要当 0 处理
zi_pi[which(zi_pi$within_module_connectivities < 2.5 & zi_pi$among_module_connectivities < 0.62),'type'] <- 'Peripherals'
zi_pi[which(zi_pi$within_module_connectivities < 2.5 & zi_pi$among_module_connectivities > 0.62),'type'] <- 'Connectors'
zi_pi[which(zi_pi$within_module_connectivities > 2.5 & zi_pi$among_module_connectivities < 0.62),'type'] <- 'Module hubs'
zi_pi[which(zi_pi$within_module_connectivities > 2.5 & zi_pi$among_module_connectivities > 0.62),'type'] <- 'Network hubs'

ggplot(zi_pi, aes(among_module_connectivities, within_module_connectivities)) +
geom_point(aes(color = type), alpha = 0.5, size = 2) +
scale_color_manual(values = c('gray','red','blue','purple'),
limits = c('Peripherals', 'Connectors', 'Module hubs', 'Network hubs'))+
theme(panel.grid = element_blank(), axis.line = element_line(colour = 'black'),
panel.background = element_blank(), legend.key = element_blank()) +
labs(x = 'Among-module connectivities', y = 'Within-module connectivities', color = '') +
geom_vline(xintercept = 0.62) +
geom_hline(yintercept = 2.5)

参考资料

Deng Y, Jiang Y H, Yang Y, et al. Molecular ecological network analyses. Bmc Bioinformatics, 2012, 13(1):113.

Guimerà R, Nunes A, Luís A. Functional cartography of complex metabolic networks. Nature, 2005, 433(7028):895-900.

Olesen J M, Bascompte J, Dupont Y L, et al. The modularity of pollination networks. Proc Natl Acad Sci U S A, 2007, 104(50):19891-19896.

Shi S, Nuccio E E, Shi Z J, et al. The interconnected rhizosphere: High network complexity dominates rhizosphere assemblages. Ecology Letters, 2016, 19(8):926-936.

友情链接

相关性和网络分析基础

网络拓扑结构-网络图的凝聚性特征和R计算 

降维分析

非约束排序（描述性的探索性分析）：

非约束排序中被动添加解释变量：被动添加解释变量

约束排序（将解释变量通过回归方程拟合响应变量的统计模型）：

RDA、CCA的变差分解（VAP） 

对称分析（这类方法意在描述两个或多个矩阵之间的相关性）：

多元因子分析（MFA）

监督降维（带监督的降维方法，常用于分类）：

判别分析（DA）

聚类和分类

层次聚类（无监督，描述性的探索性分析）：

层次聚合：层次聚合分类    层次聚类结果的比较和评估

层次分划：双向指示种分析（TWINSPAN） 

非层次聚类（无监督，描述性的探索性分析）：

划分聚类：k均值划分（k-means）    围绕中心点划分（PAM）

模糊聚类：模糊c均值聚类（FCM）

避免不存在的类

潜变量分类（无监督，潜变量也可视为某种意义上的“降维”）：

潜类别分析（LCA）    潜剖面分析（LPA）约束聚类（无监督，将解释变量通过回归方程“约束”响应变量的模型）：多元回归树 

监督分类（通过已知先验分组构建分类器模型）：

决策树    随机森林分类