接下来我们进行实例分析,首先构建两个相关网络,一个区分正、负相关,另一个不区分正、负相关(负的连接在很多算法中是没有现实意义的),具体如下所示:
代码语言:javascript复制#读取物种与环境因子数据
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()函数来实现,如下所示:
代码语言:javascript复制sub1=clusters(g1)
可以提取其结果中的子群成员、大小、子群个数信息,如下所示:
代码语言:javascript复制sub1$membership
sub1$csize
sub1$no
结果如下所示:
可以看到凡是有连接的节点都被归到同一子群,因此在相关性网络分析中较少使用。
⑵随机游走
随机游走的规则是假设一个漫游者在一个网络中随机游走(类似布朗运动),该漫游者很可能被困在连接稠密的区域里,这个稠密区就是漫游者发现的社群。为了度量漫游者的行为,从i点走向j点的概率定义为距离相似度,若i和j在同一子群,显然概率较高,因此可以采取层次聚类的方法建立子群结构。在R中随机游走可以使用walktrap.community()函数,如下所示:
代码语言:javascript复制sub2=walktrap.community(g2, weights=E(g2)$weight, step=4)
其中weights为边的权重,这里即为相关性大小,由于要计算加权概率,负的连接是会有歧义的,所以这里使用g2;step为随机游走的步长,步长越长聚类越粗糙。具体结果如下所示:
⑶基于度中心性的子群分割
首先找出网络中关联度最弱的节点(度中心性最小的点),切断他们之间的边,并以此分裂法对网络进行逐层分割,在适当的时候终止分裂过程,得到划分结果,可以使用edge.betweenness.community()函数实现,如下所示:
代码语言:javascript复制sub3=edge.betweenness.community(g2, weight=E(g2)$weight, directed=F)
此方法同样要求边的权重必须是正数,分析结果如下所示:
⑷标签传播
标签即子群名称,其基本思想是每个节点的群体归属(即标签)由与其直接连接的节点的标签决定,与其直接连接的节点属于哪个群体的最多,那么这个节点就属于那个群体。该算法通过迭代完成,可以预先设置中心节点,可以使用label.propagation.community()函数进行分析,如下所示:
代码语言:javascript复制sub4=label.propagation.community(g2, weights=V(g2)$weight)
分析结果如下所示:
此外还有贪心算法fastgreedy.community()、多层次聚类multilevel.community()、特征值leading.eigenvector.community():
代码语言:javascript复制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)
下面我们进行子群可视化:
代码语言:javascript复制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()函数进行计算,如下所示: