本文介绍了如何使用R语言复现环境因子与物种相关性的网络图，包括计算相关性、构建简单图、调整节点布局以及绘制成环径向网络图。通过手动计算节点坐标或使用create_layout()生成布局数据，实现物种与环境因子分别成环的展示效果。
https://doi.org/10.1016/j.soilbio.2020.107782
有师妹想要绘制一张类似上图的环境因子-物种相关性网络图。这张图其实还挺好复现的，将环境因子与物种都作为变量一起计算相关性指数，然后绘图时更改网络图的布局参数，将环境因子与物种分内外环放置。下面使用虚构数据绘制一幅类似的图。
一、 导入数据
物种名中的空格、-等特殊字符需要注意，导入语R中会变成".",可能导致后面合并表格，计算节点出现NA值。所以导入数据时设置check.names = FALSE，使行名和列名原样导入。计算相关性系数后，行名中的的特殊字符也会变成.,可以自行重命名一下。最好导入数据中不要存在特殊字符。
# 1.1 设置工作路径、加载R包
setwd("D:\\EnvStat\\network\\layout")
library(tidyverse)
library(igraph)
library(ggraph)
library(tidygraph)
library(reshape2)
library(ggplot2)
library(psych)
# 1.1 导入数据：观测-变量数据表
data<- read.csv("data.csv",header = TRUE,
row.names = 1,
check.names=FALSE,
comment.char = "")
head(data)
# 1.2 变量分类表
type = read.csv("type.csv",header = TRUE,check.names = FALSE)
head(type) # 变量分类信息
图1|环境因子-物种丰度数据，data.csv。行为样本，列为变量。
图2|变量分类信息表，type.csv。
二、 确定相关性关系
2.1 计算相关性系数
# 2.1 计算相关性系数，建议使用spearman系数。
cor <- psych::corr.test(data, use = "pairwise",
method="spearman",
adjust="fdr",
alpha=.05,ci=FALSE) # ci=FALSE,不进行置信区间计算，数据量较大时，可以加快计算速度。
cor.r <- data.frame(cor$r) # 提取R值
cor.p <- data.frame(cor$p) # 提取p值
colnames(cor.r) = rownames(cor.r)
colnames(cor.p) = rownames(cor.p) # 变量名称中存在特殊字符，为了防止矩阵行名与列名不一致，必须运行此代码。
write.csv(cor.r,"cor.r.csv",quote = FALSE,col.names = NA,row.names = TRUE) # 保存结果到本地
write.csv(cor.p,"cor.p.csv",quote = FALSE,col.names = NA,row.names = TRUE)
head(cor.r)
head(cor.p)2.2 确定相关性关系
这里以p<0.05作为确定相关性系数的阈值，不对r值进行筛选。也可以使用MENA(http://129.15.40.240/mena/)筛选r值。
# 2.2 确定阈值：保留p<0.05的物种间相关关系
## 2.2.1 将数据转换为long format进行合并，过滤，并添加链接属性
cor.r$from = rownames(cor.r)
cor.p$from = rownames(cor.p)
p = cor.p %>%
gather(key = "to", value = "p", -from) %>%
data.frame()
cor.data = cor.r %>%
gather(key = "to", value = "r", -from) %>%
data.frame() %>%
left_join(p, by=c("from","to")) %>%
filter(p <= 0.05, from != to) %>%
mutate(
linecolor = ifelse(r > 0,"positive","negative"), # 设置链接线属性，可用于设置线型和颜色。
linesize = abs(r) # 设置链接线宽度。
) # 此输出仍有重复链接，后面需进一步去除。
head(cor.data)
图3|相关性链接表，cor.data。data存在重复的链接(例AK-pH，pH-AK)，后面还需要过滤。
三、 构建网络
上步骤生成的cor.data中，存在重复的链接(例AK-pH，pH-AK)，用此数据构建的网络图，不是简单图。所以依此计算的节点degree有误，后面需要转为简单图后，重新计算节点degree。
3.1 准备节点属性文件
# 3.1 准备节点属性文件
## 3.1.1 计算每个节点具有的链接数(degree)
c(as.character(cor.data$from),as.character(cor.data$to)) %>%
as_tibble() %>%
group_by(value) %>%
summarize(n=n()) -> vertices
colnames(vertices) <- c("name", "n")
## 3.1.2 添加变量分类属性
vertices <- vertices %>%
select(-n) %>% # 因为此处的n不准确，所以先删除。
left_join(type,by="name")
## 3.1.3 对节点属性表进行排序
#网络图中节点会按照节点属性文件的顺序依次绘制，
##为了使同类型变量位置靠近，按照节点属性对节点进行排序。
##环境变量最好置于表格末尾，否则后续绘图可能会效果不好
##此处使用type进行排序，依次是细菌、真菌和环境因子,
##然后按照group排序，使同一属的物种彼此靠近。
vertices$type <- factor(vertices$type,
levels = c("bacteria","Fungus","Soil physi-chemical factors" ))
vertices <- vertices %>%
arrange(type,group)
dim(vertices)
head(vertices)
图4|节点属性表，vertices。节点按照type，group依次排序，环境因子排在最后。
3.2 构建简单图(graph)数据结构
## 3.2.1 构建graph数据结构
graph <- graph_from_data_frame(cor.data, vertices = vertices, directed = FALSE )
graph # 非简单图
图5|非简单图，graph。
## 3.2.2 转为简单图，去除cor.data数据表中的重复链接
is.simple(graph) # 非简单图，链接数会偏高，所以需要转换为简单图。
E(graph)$weight <- 1 # 将链接权重赋值为1
graph <- igraph::simplify(graph,edge.attr.comb = "first")# 转为简单图
is.simple(graph)
E(graph)$weight <- 1 # 将链接权重赋值为1
is.weighted(graph)
graph
图6|简单图，graph。去除自身链接形成的loops和两个节点互相链接形成的多边。
3.3 计算每个节点具有的链接数(degree)
# 3.3 重新计算每个节点具有的链接数(degree)
V(graph)$degree <- degree(graph)3.4 保存绘图数据到本地
此步从graph提取的节点和链接属性表，后面绘图分析需要用到。
# 3.4 保存绘图数据到本地
write.graph(graph,file = "graph.gml",format="gml") # 直接保存graph结构，gml能保存的graph信息最多。
net.data <- igraph::as_data_frame(graph, what = "both")$edges # 提取链接属性
write.csv(net.data,"net.data.csv",quote = FALSE,col.names = NA,row.names = FALSE) # 保存结果到本地。
head(net.data)
vertices <- igraph::as_data_frame(graph, what = "both")$vertices # 提取节点属性
write.csv(vertices,"vertices.csv",quote = FALSE,col.names = NA,row.names = FALSE)
head(vertices) # 直接读入此步保存的链接和节点属性文件，之后可直接生成graph或用于其他绘图软件绘图。
图7|提取的简单图节点属性表，vertices。包含计算的degree。
四、 绘制网络图
复现的径向网络图中，物种与环境因子分别单独成环，因此需要分别计算节点标签的位置，然后再重构网络图。此处使用ggraph绘制物种-环境因子分别成环相关性网络图，则需要更改ggraph()默认的节点布局参数，即节点的(x,y)坐标。
有两种方式可以更改布局: 1）手动计算节点坐标，然后绘图前使用create_layout()生成布局数据表，更改其中的x,y列，然后绘图时直接使用此布局文件。2)提取物种或环境因子子图，分别生成布局数据表，获取坐标后，更改create_layout()为全部网络节点生成的(x，y)或更改绘图后的data数据(例，net.cir$data)中的(x,y)。
两种方法中的重点就是为每个节点生成理想的坐标。下面介绍两种生成物种-环境因子分别成环相关性网络图节点坐标的方法。
4.1 设置颜色
# 4.1 设置颜色
## 4.1.1 节点分组颜色
unique(vertices$group) # 7个分组信息
#??ggsci # 查看ggsci包帮助文件，选择颜色。
library(ggsci)
library(scales)
mycolor = pal_d3("category20",alpha = 1)(20)
mycolor
cols = mycolor[c(1:7,9:13,14)]
show_col(cols) # 以图形形式展示颜色
## 4.1.2 链接线颜色
#linetypes = c("positive" ="solid","negative" ="dashed") # 设置线类型
color = c("positive" ="#D62728FF","negative" ="#2CA02CFF")
图8|cols颜色集。
4.2 方法1: 手动计算节点坐标绘图
此图的两个环其实就是物种或环境因子节点均匀的分布在两个半径不同的圆的周长上。即外环，每个物种所占圆的角度是1/物种总数，半径设为1，则排序第一的物种节点的坐标应该为(cos(1/物种总数)，sin(1/物种总数)),起始分子可以更改，但后面物种的角度则都是分子依次+1。添加外环则可增大半径，添加内环则减小半径。依据此原理，可分别手动计算环境因子或物种的网络图坐标。
# 4.2 方法1: 手动计算节点坐标
##添加标签角度数据
## seq()生成一组1-数据表行数的数值。
#number_of_bar<-nrow(vertices)
#vertices$id = seq(1, nrow(vertices))
#angle= 360 * (vertices$id-0.5) /number_of_bar
#vertices$hjust<-ifelse(angle>180, 1, 0)
#vertices$angle<-ifelse(angle>180, 90-angle+180, 90-angle)
##物种与环境因子位于一个环上，则运行上面的代码。
## 4.2.1 为物种和环境因子标签分别生成角度
### 计算各类型变量总数
bar1 = nrow(vertices[vertices$group == "env",])
bar1 # 9个环境因子
bar2 = nrow(vertices[vertices$group != "env",])
bar2 # 17个物种
### 为各变量添加序号
vertices$id =1
vertices$id[vertices$group == "env"] = seq(1, bar1)
vertices$id[vertices$group != "env"] = seq(1, bar2)
head(vertices)
### 计算变量与轴的角度
vertices$var_angel <- NA
vertices$var_angel[vertices$group == "env"] <- 360 * (vertices$id[vertices$group == "env"]-0.5)/bar1
vertices$var_angel[vertices$group != "env"] <- 360 * (vertices$id[vertices$group != "env"]-0.5)/bar2
### 计算物种和环境因子坐标-物种环半径设为1，环境因子环半径为0.5。
vertices$x <- ifelse(vertices$group == "env",cos(vertices$var_angel)/2,cos(vertices$var_angel))
vertices$y <- ifelse(vertices$group == "env",sin(vertices$var_angel)/2,sin(vertices$var_angel))
head(vertices)
## 4.2.2 为物种和环境因子标签添加hjust和angle属性
vertices$hjust <- ifelse(vertices$var_angel > 180, 1, 0)
vertices$angle <- ifelse(vertices$var_angel > 180, 90-vertices$var_angel+180, 90-vertices$var_angel)
## 4.2.3 内环的环境因子标签，设置为在节点中心水平放置
vertices$hjust[vertices$group == "env"] = "center" # 或0.5
vertices$angle[vertices$group == "env"] = 0
write.csv(vertices,"vertices.csv",quote = FALSE,row.names = FALSE) #保存结果到本地
head(vertices)
图9|新vertices节点属性表。新增节点标签位置属性和节点(x，y)。
## 4.2.4 将group属性加入链接属性表
net.data <- net.data %>%
left_join(vertices,by=c('from'='name'))
write.csv(net.data,"net.data.csv",quote = FALSE,row.names = FALSE) # 保存结果到本地
head(net.data)
## 4.2.5 重新构建graph结构数据
graph <- graph_from_data_frame(net.data, vertices = vertices, directed = FALSE)
graph
图10|最终graph结构文件，graph。
## 4.2.6 create_layout()创建布局数据，并将其中的(x,y)替换为上面计算的坐标。
layout <- create_layout(graph,layout = 'linear', circular = TRUE)
head(layout)
layout$x <- vertices$x
layout$y <- vertices$y
head(layout)
图11|更改坐标前的graph布局信息，layout。
图12|更改坐标后的graph布局信息，layout。
下面直接使用创建的layout绘图即可。
## 4.2.7 绘制物种-环境因子分别成环径向网络图1
### 使用上步生成的layout布局绘图
##根据物种属设置节点颜色；
##节点degree设置节点大小；
##除环境因子(黑色)外，节点标签使用节点颜色设置；
##相关性正负设置链接线颜色；
##使用r绝对值大小设置链接线粗细。
net.cir1 <- ggraph(layout) +
geom_edge_link(aes(edge_colour = as.factor(linecolor),
edge_width = abs(r)
#edge_linetype = linetype # 想用不同线形表示相关性正负，取消此句注释
),
#edge_width=0.5,
edge_alpha=0.9 ) + # 使用直线作为链接线。
###theme()中的legend参数不利于多个图例设置，这里在scale*函数中直接设置图例参数###
scale_edge_colour_manual(values=color,
breaks = c("positive","negative"),
guide = guide_legend(title="Correlation",
direction="vertical",
order=2,# 图例排列序号
ncol=1,
byrow=FALSE,
title.theme = element_text(
size = 14,
face = "bold",
colour = "black"),
label.theme = element_text(
size = 12,colour = "black")
)) +
scale_edge_width(
breaks = seq(0.2,1,0.2),
label = seq(0.2,1,0.2),
range = c(0.2,1),
guide = guide_legend(title="|r|",
direction="vertical",
order=1,# 图例排列序号
ncol=1,
byrow=FALSE,
title.theme = element_text(
size = 14,
face = "bold",
colour = "black"),
label.theme = element_text(
size = 12,colour = "black")
))+
###想用不同线形表示相关性正负，取消此句注释###
#scale_edge_linetype_manual(values = linetypes,
#name = "Correlation",
#breaks = c("positive","negative"))+
geom_node_point(aes(size=degree, fill=as.factor(group)),
shape=21,alpha=1) +
#scale_size_continuous(name = "Degree",range=c(10,15)) +
scale_size(
breaks = seq(3,max(vertices$degree),3),
label = seq(3,max(vertices$degree),3),
range = c(5, 15),
guide = guide_legend(title="Degree",
direction="vertical",
order=3,# 图例排列序号
ncol=2, # 图例排两列
byrow=FALSE,
title.theme = element_text(
size = 14,
face = "bold",
colour = "black"),
label.theme = element_text(
size = 12,colour = "black")
))+
###可使用cols = c(...,"env"="black")为每个节点设置指定颜色，适合节点分类少的网络图###
scale_fill_manual(values=cols,
guide = guide_legend(title="Genus",
#keywidth=grid::unit(2,"cm"),
#keyheight=grid::unit(2,"cm"),
direction="vertical",
order=4,# 图例排列序号
ncol=2,
byrow=FALSE,
title.theme = element_text(
size = 14,
face = "bold",
colour = "black"),
label.theme = element_text(
size = 12,colour = "black",face = "italic")
)) +
###展示物种+环境因子节点标签，不展示节点标签，则注释掉geom_node_text()函数###
#geom_node_text(aes(x = x*1.1, y=y*1.1,
#label=as.character(name),
#angle=angle,hjust=hjust,
#color=as.factor(group)),size=4.5,
#show.legend = FALSE) +
###只展示环境因子节点标签###
geom_node_text(aes(x = x, y=y,
label=ifelse(V(graph)$group == "env",as.character(name),""),
angle=angle,hjust=hjust,
color = as.factor(group)
),
size=4.75,# 14pt
show.legend = FALSE) +
###将env分类对应的颜色改为“black”,更改节点标签颜色###
scale_color_manual(values = c(cols[1:6],"black"))+
#expand_limits(x = c(-1, 1), y = c(-1, 1))+ # 设置网络图区域坐标范围
coord_fixed()+
theme_minimal()+
theme(
legend.title.align = 0,# 图例标题左对齐方式(0,左对齐；1，右对齐；0.5，居中对齐)
legend.text.align = 0,# 图例文字左侧对齐方式(0,左对齐；1，右对齐；0.5，居中对齐)
###图例的左上角位于图中(1,0.9)坐标处###
#legend.position = c(1,0.9), # 坐标设置图例位置左→右：0-1；下→上：0-1.
#legend.justification = c(0,1), # 更改图例锚点(左→右：0-1；下→上：0-1,默认(0,1))
legend.spacing.x = unit(0.1,"cm"), # 图例中各个元素(key与key标签)的距离
legend.spacing.y = unit(0.1,"cm"),
legend.box.margin = margin(t=0.2,r=0.2,b=0.2,l=0.2,"cm"), # 设置图例中心与图例框边缘的距离
plot.margin=unit(rep(0,4), "null"), # 设置图片区与边界的距离,值越大图与图片边界距离越远。
panel.spacing=unit(c(0,0,0,0), "null"),
###去除图中坐标轴和网格线###
panel.grid = element_blank(),
axis.line = element_blank(),
axis.ticks =element_blank(),
axis.text =element_blank(),
axis.title = element_blank()
)
### 输出pdf格式图片到本地
cairo_pdf("net.cir1.pdf",height=8,width=10,family="Times")
print(net.cir1) # 图中的节点分布不是特别均匀，而且计算过程有些复杂，推荐使用第二种方法更省事和更美观。
dev.off()
图13|环境因子-物种分别成环径向网络图，net.cir1.pdf。连接线设置为直线，此方法绘制的点不太均匀，推荐方法2.
4.3 方法二：creat_layout()分别计算物种和环境因子的坐标
# 4.3 方法2：creat_layout()分别计算物种和环境因子的坐标
## 4.3.1 计算环境因子、物种子网络径向布局数据表
###使用induced.subgraph()提取环境因子、物种子网络，
###create_layout()分别计算两者的径向布局(x,y)坐标。
### 此种方法环境因子必须全部在节点的末尾(多次尝试发现，内环的变量不位于末尾，画出来效果就不好，所以前面步骤应该提前排序)。
### 此方法要添加节点标签，也需要计算节点水平位置和角度，只是不需要计算(x,y),前面方法1已经计算过了，这里不再重复。
### 计算物种子网络节点布局数据表
layout1 <- create_layout(
induced.subgraph(graph,c(V(graph)$name %in% vertices$name[vertices$group != "env"])),
layout = 'linear', circular = TRUE)
layout1
### 计算环境因子子网络节点布局数据表
layout2 <- create_layout(
induced.subgraph(graph,c(V(graph)$name %in% vertices$name[vertices$group == "env"])),
layout = 'linear', circular = TRUE)
layout2[1:2] <- layout2[1:2]/2 # 内环变量的坐标同时除2。
## 4.3.2 create_layout()创建布局数据，并将其中的(x,y)替换为上面计算的坐标。
layout <- create_layout(graph,layout = 'linear', circular = TRUE)
head(layout)
## 4.3.3 替换坐标-要求name的顺序一致。
###rbind(layout2,layout1)，不能用于设置ggraph(graph,layout = layout),绘图时会缺少物种链接线布局信息。
layout$x <- rbind(layout1,layout2)$x
layout$y <- rbind(layout1,layout2)$y
head(layout)
图14|方法2更改坐标后的graph布局信息，layout。
# 4.3.4 绘制径向布局网络图2
##根据物种属设置节点颜色；
##节点degree设置节点大小；
##相关性正负设置链接线颜色；
##使用r绝对值大小设置链接线粗细。
Correlation <- guide_legend(title="Correlation",
direction="vertical",
order=2,# 图例排列序号
ncol=1,
byrow=FALSE,
title.theme = element_text(
size = 14,
face = "bold",
colour = "black"),
label.theme = element_text(
size = 12,colour = "black")
)
width_legend <- guide_legend(title="|r|",
direction="vertical",
order=1,# 图例排列序号
ncol=1,
byrow=FALSE,
title.theme = element_text(
size = 14,
face = "bold",
colour = "black"),
label.theme = element_text(
size = 12,colour = "black")
)
size_legend <- guide_legend(title="Degree",
direction="vertical",
order=3,# 图例排列序号
ncol=2, # 图例排两列
byrow=FALSE,
title.theme = element_text(
size = 14,
face = "bold",
colour = "black"),
label.theme = element_text(
size = 12,colour = "black")
)
fill_legend <- guide_legend(title="Genus",
direction="vertical",
order=4,# 图例排列序号
ncol=2,
byrow=FALSE,
title.theme = element_text(
size = 14,
face = "bold",
colour = "black"),
label.theme = element_text(
size = 12,
colour = "black",face = "italic")
)
set_graph_style(plot_margin = margin(0,0,1,0))
net.cir2 <- ggraph(layout) +
geom_edge_arc(aes(edge_colour = as.factor(linecolor),
edge_width = abs(r)
),
edge_alpha=0.9 ) + # 改用曲线链接。
scale_edge_colour_manual(values=color,
breaks = c("positive","negative"),
guide = Correlation) +
scale_edge_width(
breaks = seq(0.2,1,0.2),
label = seq(0.2,1,0.2),
range = c(0.2,1),
guide = width_legend )+
geom_node_point(aes(size=degree, fill=as.factor(group)),
shape=21,alpha=1) +
scale_fill_manual(values=cols,
guide = fill_legend) +
scale_size(
breaks = seq(3,max(vertices$degree),3),
label = seq(3,max(vertices$degree),3),
range = c(5, 15),
guide = size_legend)+
geom_node_text(aes(x = x, y=y,
label= ifelse(V(graph)$group == "env",as.character(name),""),
angle=angle,hjust=hjust,
#color = as.factor(group)
),
color = "black",# 所有节点标签设置为黑色。
size=4.75,# 14pt
show.legend = FALSE) +
##添加物种标签，取消下面代码的注释即可##
#geom_node_text(aes(x = x*1.4, y=y*1.4, # 物种与环境因子标签的放置方式不同，最好分别设置，
#label= ifelse(V(graph)$group != "env",as.character(name),""),
#angle=angle,hjust=hjust,
# ),
#color = "black",# 所有节点标签设置为黑色。
#size=4.75,# 14pt
#show.legend = FALSE) +
coord_fixed()+
theme_graph()
### 输出pdf格式图片到本地-输出到本地之后，可用AI或pdf编辑器删除物种标签，初期探索数据的时候，把标签加上比较好。
cairo_pdf("net.cir2.pdf",height=8,width=10,family="Times")
print(net.cir2) # 节点标签过长，可只展示种名。
dev.off()
图15|环境因子-物种分别成环径向网络图，net.cir2.pdf。连接线设置为曲线，此为未经AI修改的原图。
参考资料：
https://www.statology.org/ggplot2-legend-size/
单个图例位置修改：https://blog.csdn.net/lsxxx2011/article/details/98764286
图例设置：http://www.360doc.com/content/19/0805/18/52645714_853164588.shtml
绘图区拓展：https://www.5axxw.com/questions/content/w63e7x
http://127.0.0.1:12026/library/ggraph/doc/Layouts.html
EcoEvoPhylo公众号后台发送“network_layout”或QQ群文件获取数据和代码。
公众号文章链接：R绘图-物种、环境因子相关性网络图(简单图、提取子图、修改图布局参数、物种-环境因子分别成环径向网络图)
EcoEvoPhylo ：主要分享微生物生态和phylogenomics的数据分析教程。扫描上方二维码，关注EcoEvoPhylo。让我们大家一起学习，互相交流，共同进步。
知乎 | 生态媛
学术交流QQ群 | 438942456
学术交流微信群 | jingmorensheng412
加好友时，请备注姓名-单位-研究方向。
推荐阅读
R统计绘图-分子生态相关性网络分析(拓扑属性计算，ggraph绘图)
R统计-PCA/PCoA/db-RDA/NMDS/CA/CCA/DCA等排序分析教程
R统计绘图-环境因子相关性+mantel检验组合图(linkET包介绍1)
机器学习-分类随机森林分析(randomForest模型构建、参数调优、特征变量筛选、模型评估和基础理论等)
R统计-单因素ANOVA/Kruskal-Wallis置换检验
R统计绘图-单、双、三因素重复测量方差分析[Translation]
R统计绘图-多变量相关性散点矩阵图(GGally::ggpairs())
5217
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
