相关性网络的子群划分

SYSU星空

发布于 2022-05-05 12:54:58

5261

发布于 2022-05-05 12:54:58

虽然基于相关性矩阵的聚类热图可以对物种或者环境因子进行聚类，但其原理主要是把行为相似（或者理解为共线性）的聚在一起。而基于网络的聚类手段更加强大，能够将相关网络划分为一个个子群体（community），子群内的个体之间关联度要显著强于与子群之外的个体。

接下来我们进行实例分析，首先构建两个相关网络，一个区分正、负相关，另一个不区分正、负相关（负的连接在很多算法中是没有现实意义的），具体如下所示：

#读取物种与环境因子数据
community=read.table(file="otu_table_L4.txt", header=TRUE)
rownames(community)=community[,1]
community=community[,-1]
com=t(community)
environment=read.table(file="environment.txt", header=TRUE)
rownames(environment)=environment[,1]
env=environment[,-1]
env=env[rownames(com),]
data=as.matrix(cbind(com, env))
#计算相关性矩阵并筛选p<0.05、r>0.45的数据
library(Hmisc)
corr=rcorr(data, type="spearman")
rcorr=corr$r #提取相关系数
pcorr=corr$P #提取检验结果p值
#注意要去掉自相关（对角线上的数据）
for (i in 1:nrow(rcorr)) {
  for (j in 1:nrow(rcorr)) {
    if (i!=j) {
      if (pcorr[i, j]>0.05) {
        rcorr[i, j]=0
      }
    }
    if (rcorr[i, j]>-0.45 & rcorr[i, j]<0.45) {
      rcorr[i, j]=0
    }
  }
}
#构建相关网络
library(igraph)
g1=graph.adjacency(rcorr, mode="undirected", weighted=T, diag=F)
g2=graph.adjacency(abs(rcorr), mode="undirected", weighted=T, diag=F)

接下来进行子群划分，不同算法具体如下：

⑴基于点连接的社群发现

两个节点之间有连接则视为同一个子群，正、负关联度没有影响，可以使用clusters()函数来实现，如下所示：

sub1=clusters(g1)

可以提取其结果中的子群成员、大小、子群个数信息，如下所示：

sub1$membership
sub1$csize
sub1$no

结果如下所示：

可以看到凡是有连接的节点都被归到同一子群，因此在相关性网络分析中较少使用。

⑵随机游走

随机游走的规则是假设一个漫游者在一个网络中随机游走（类似布朗运动），该漫游者很可能被困在连接稠密的区域里，这个稠密区就是漫游者发现的社群。为了度量漫游者的行为，从i点走向j点的概率定义为距离相似度，若i和j在同一子群，显然概率较高，因此可以采取层次聚类的方法建立子群结构。在R中随机游走可以使用walktrap.community()函数，如下所示：

sub2=walktrap.community(g2, weights=E(g2)$weight, step=4)

其中weights为边的权重，这里即为相关性大小，由于要计算加权概率，负的连接是会有歧义的，所以这里使用g2；step为随机游走的步长，步长越长聚类越粗糙。具体结果如下所示：

⑶基于度中心性的子群分割

首先找出网络中关联度最弱的节点（度中心性最小的点），切断他们之间的边，并以此分裂法对网络进行逐层分割，在适当的时候终止分裂过程，得到划分结果，可以使用edge.betweenness.community()函数实现，如下所示：

sub3=edge.betweenness.community(g2, weight=E(g2)$weight, directed=F)

此方法同样要求边的权重必须是正数，分析结果如下所示：

⑷标签传播

标签即子群名称，其基本思想是每个节点的群体归属（即标签）由与其直接连接的节点的标签决定，与其直接连接的节点属于哪个群体的最多，那么这个节点就属于那个群体。该算法通过迭代完成，可以预先设置中心节点，可以使用label.propagation.community()函数进行分析，如下所示：

sub4=label.propagation.community(g2, weights=V(g2)$weight)

分析结果如下所示：

此外还有贪心算法fastgreedy.community()、多层次聚类multilevel.community()、特征值leading.eigenvector.community()：

sub5=leading.eigenvector.community(g2, weights=V(g2)$weight)
sub6=fastgreedy.community(g2, weights=V(g2)$weight)
sub7=multilevel.community(g2, weights=V(g2)$weight)

下面我们进行子群可视化：

m=length(colnames(com))
n=length(colnames(env))
t=length(colnames(rcorr))
size1=numeric(m)
for (i in 1:m) {
  size1[i]=sum(com[,i])
}
size1=(10000*size1)^0.25
size2=numeric(n)
for (i in 1:n) {
  size2[i]=sum(abs(rcorr[,m+i]))-1
}
size2=1.7*sqrt(10*size2)
size=c(size1, size2)
myshape=c(rep("circle", m), rep("square", n))
weight=E(g2)$weight
par(mar=c(0.4,0.4,1,0.4), mfrow=c(2,2))
plot(sub2, g2, vertex.size=size*0.7,vertex.label=NA, vertex.shape=myshape, edge.width=2*weight^2, main="随机游走")
plot(sub3, g2, vertex.size=size*0.7,vertex.label=NA, vertex.shape=myshape, edge.width=2*weight^2, main="度中心性")
plot(sub4, g2, vertex.size=size*0.7,vertex.label=NA, vertex.shape=myshape, edge.width=2*weight^2, main="标签传播")
plot(sub5, g2, vertex.size=size*0.7,vertex.label=NA, vertex.shape=myshape, edge.width=2*weight^2, main="特征值")
plot(sub6, g2, vertex.size=size*0.7,vertex.label=NA, vertex.shape=myshape, edge.width=2*weight^2, main="贪心策略")
plot(sub7, g2, vertex.size=size*0.7,vertex.label=NA, vertex.shape=myshape, edge.width=2*weight^2, main="多层次聚类")

作图结果如下所示：

其中子群之间的连接为红色，子群内部的连接为黑色，不同子群节点以颜色区分，接下来对网络及其子群进行比较衡量分析。

⑴网络聚类系数

网络聚类系数也即聚集系数是对网络中节点聚集程度的衡量，值越大表示网络中节点关联性越强，网络结构越复杂，可以使用transitivity()函数进行计算，如下所示：

⑵网络密度

跟网路聚类系数相似，也是用来形容网络的结构复杂程度。其值越大，说明网络越复杂，可以使用graph.density()函数进行计算，如下所示：

⑶模块化度量值

模块化度量值Q（Modularity)是衡量子群稳定度的指标，其中第i个子群的值计算公式为：Qi=I/E-((2I+O)/2E)2，其中I表示两个节点均在该子群中的边的数目，E为两个节点均不在该子群的边的数目，O表示其中一个端点在该子群中，而另一个端点不在该子群中的边的数目，所有子群的值相加得到Q值。模块化度量值越大说明子群划分越明显，社群之间关联越少，可以使用modularity()函数进行计算，如下所示：