承接前一篇文章，接下来我们利用复杂网络理论对相关网络数据进行深入的分析。在网络分析中的节点度（node degree）是指和该节点关联的边的条数，或者说连接的个数，又称关联度；显然网络节点越多，节点度越大，为了去除网络规模的影响，使得不同网络可以相互比较，可以使用度中心性（degree centrality）概念。度中心性是在网络分析中刻画节点中心性的最直接度量指标，其值为该节点节点度除以该节点最大可能节点度，也即该节点实际连接数占与其他节点可能连接总数目的比例，如下所示：
其中g为节点总数，度中心性取值范围0-1，为0表明与所有其他节点都无关联，为1表示与其他所有节点有直接关联。在一个网络中节点度越大就意味着这个节点的度中心性越高，该节点在网络中就越重要。此外，还有其他指标例如接近度中心性（closeness centrality，该点与网络中其他点距离之和的倒数，越大说明越在中心）、中介系数中心性（betweennesscentrality，一个点在多大程度上位于图中其他“点对”的“中间”）、特征向量中心性（eigenvectorcentrality，通过相邻点的重要性来衡量该节点的重要性）等，在相关网络中一般使用不到，相关网络中也可以使用加权的节点度（也即相关系数绝对值之和）作为衡量指标。下面以上一篇文章中构建的相关网络为例进行分析：
代码语言：javascript复制#计算节点度
diag(rcorr)=0
degree=numeric(ncol(rcorr))
for (i in 1:ncol(rcorr)) {
degree[i]=length(which(rcorr[i,]!=0))
}
#计算度中心性
centdeg=degree/(ncol(rcorr)-1)
nodedata=data.frame(degree, centdeg, row.names=rownames(rcorr))
nodedata=nodedata[order(nodedata[,2],decreasing=TRUE),]
#节点度分布图
library(ggplot2)
ggplot(nodedata, aes(x=degree))
geom_histogram(position='identity', alpha=0.5, fill="cyan", binwidth=1, aes(y=..density..))
stat_density(geom='line', position='identity', col="cyan4")
xlab("Node Degree") ylab("Density")此外，可以直接使用函数计算节点度等指标：
degree(g)#计算节点度
closeness(g)#计算接近度中心性
betweenness(g)#中介系数中心性
evcent(g)#计算特征向量中心性
节点度分布图是不同节点度范围内的节点数目统计情况，可以反映网络的异质性，也即节点之间的连接状况是否均匀，理论上高关联度节点越多网络结构越复杂，做图结果如下所示：
接下来我们可以筛选出度中心性高的节点，来看那些物种或者环境因子在相关性网络中的影响较大：
代码语言：javascript复制#节点度中心性条形图
nodedata=nodedata[1:20, ]
ggplot(nodedata, aes(x=factor(rownames(nodedata), levels=rev(rownames(nodedata))), y=centdeg))
geom_bar(stat="identity", fill="cyan", alpha=0.6)
xlab("Node") ylab("Node degree centrality") coord_flip()做图结果如下所示：
接下来，我们可以筛选受环境因子直接影响（相关系数之和不为0）的物种，并提取其相对丰度信息以便进行比较分析：
代码语言：javascript复制#提取筛选环境因子与物种相关性
envcor=rcorr[1:m, (m 1):(m n)]
sumcor=numeric(m)
for (i in 1:m) {
sumcor[i]=sum(abs(envcor[i,]))
}
ecocor=data.frame(cbind(envcor, sumcor))
ecocor=ecocor[order(abs(ecocor$sumcor), decreasing=TRUE),]
abund=numeric(m)
for (i in 1:m) {
abund[i]=100*mean(com[, i])
}
names(abund)=colnames(com)
abundance=abund[rownames(ecocor)]
cordata=data.frame(ecocor, abundance)
cordata=cordata[cordata$sumcor!=0,]
write.csv(cordata, file="cordata.csv", row.names=TRUE, quote=FALSE)经整理结果如下所示：
