【原】ggClusterNet包文档更新-增加30%以上网络分析相关内容（2023年2月）

微生信生物 2023-02-13 发布于北京

展开全文

R包安装

导入R包

Hide


#--导入所需R包#-------
library(phyloseq)
library(igraph)
library(network)
library(sna)
library(tidyverse)
library(ggClusterNet)

输入方式

输入方式一

直接输入phyloseq格式的数据。

Hide

#-----导入数据#-------
data(ps)

输入方式二

可以从https://github.com/taowenmicro/R-_function下载数据，构造phylsoeq文件。自己的数据也按照网址示例数据进行准备。虽然phylsoeq对象不易用常规手段处理，但是组学数据由于数据量比较大，数据注释内容多样化，所以往往使用诸如phyloseq这类对象进行处理，并简化数据处理过程。ggClusterNet同样使用了phyloseq对象作为微生物网络的分析。

phyloseq对象构建过程如下，网络分析主要用到otu表格，后续pipeline流程可能用到分组文件metadata，如果按照分类水平山色或者区分模块则需要taxonomy。这几个部分并不是都必须加入phyloseq对象中，可以用到那个加那个。

或者直接从网站读取，但是由于github在国外，所以容易失败

library(phyloseq)library(ggClusterNet)library(tidyverse)library(Biostrings)
metadata = read.delim("./metadata.tsv",row.names = 1)otutab = read.delim("./otutab.txt", row.names=1)taxonomy = read.table("./taxonomy.txt", row.names=1,header = T)head(taxonomy)# tree  = read_tree("./otus.tree")# rep = readDNAStringSet("./otus.fa")
ps = phyloseq(sample_data(metadata),              otu_table(as.matrix(otutab), taxa_are_rows=TRUE),              tax_table(as.matrix(taxonomy))#,              # phy_tree(tree),              # refseq(rep)              )
psrank.names(ps)

#或者直接从网站读取，但是由于github在国外，所以容易失败
metadata = read.delim("https://raw./taowenmicro/R-_function/main/metadata.tsv",row.names = 1)otutab = read.delim("https://raw./taowenmicro/R-_function/main/otutab.txt", row.names=1)taxonomy = read.table("https://raw./taowenmicro/R-_function/main/taxonomy.txt", row.names=1)# tree  = read_tree("https://raw./taowenmicro/R-_function/main/otus.tree")# rep = readDNAStringSet("https://raw./taowenmicro/R-_function/main/otus.fa")
ps = phyloseq(sample_data(metadata),              otu_table(as.matrix(otutab), taxa_are_rows=TRUE),              tax_table(as.matrix(taxonomy))#,              # phy_tree(tree),              # refseq(rep)              )

使用ggCLusterNet进行网络分析的过程

corMicro函数用于计算网络相关

按照丰度过滤微生物表格，并却计算相关矩阵，按照指定的阈值挑选矩阵中展示的数值。调用了psych包中的corr.test函数，使用三种相关方法。N参数提取丰度最高的150个OTU；method.scale参数确定微生物组数据的标准化方式，这里我们选用TMM方法标准化微生物数据。

Hide


#-提取丰度最高的指定数量的otu进行构建网络


#----------计算相关#----
result = corMicro(ps = ps,
                   N = 150,
                   method.scale = "TMM",
                   r.threshold=0.8,
                   p.threshold=0.05,
                   method = "spearman"
                   
                   
                   )

#--提取相关矩阵
cor = result[[1]]
cor[1:6,1:6]
#>         ASV_1 ASV_28 ASV_125 ASV_135 ASV_8 ASV_106
#> ASV_1       1      0       0       0     0       0
#> ASV_28      0      1       0       0     0       0
#> ASV_125     0      0       1       0     0       0
#> ASV_135     0      0       0       1     0       0
#> ASV_8       0      0       0       0     1       0
#> ASV_106     0      0       0       0     0       1

#-网络中包含的OTU的phyloseq文件提取
ps_net = result[[3]]

#-导出otu表格
otu_table = ps_net %>% 
  vegan_otu() %>%
  t() %>%
  as.data.frame()

制作分组,我们模拟五个分组

这是网络布局的基础，无论是什么聚类布局，都需要制作一个分组文件，这个文件有两列，一列是节点，一列是分组信息，这个分组信息名称为：group。这个文件信息就是用于对节点进行分组，然后按照分组对节点归类，使用包中可视化函数计算节点位置。

注意分组文件的格式，分为两列，第一列是网络中包含的OTU的名字，第二列是分组信息，同样的分组标记同样的字符。

Hide

#--人工构造分组信息：将网络中全部OTU分为五个部分，等分
netClu = data.frame(ID = row.names(otu_table),group =rep(1:5,length(row.names(otu_table)))[1:length(row.names(otu_table))] )
netClu$group = as.factor(netClu$group)
head(netClu)

	ID <chr>	group <fct>
1	ASV_1	1
2	ASV_28	2
3	ASV_125	3
4	ASV_135	4
5	ASV_8	5
6	ASV_106	1

6 rows

PolygonClusterG 根据分组，计算布局位置坐标

不同的模块按照分组聚集成不同的圆，并且圆形的大小一样。如果一个分组只有一个点，则这个点坐落在圆心位置。

Hide

#--------计算布局#---------
result2 = PolygonClusterG (cor = cor,nodeGroup =netClu  )
node = result2[[1]]
head(node)

	X1 <dbl>	X2 <dbl>	elements <chr>
ASV_1	0.000000e+00	10.5	ASV_1
ASV_14	2.598076e+00	9.0	ASV_14
ASV_72	2.598076e+00	6.0	ASV_72
ASV_109	3.673819e-16	4.5	ASV_109
ASV_71	-2.598076e+00	6.0	ASV_71
ASV_88	-2.598076e+00	9.0	ASV_88

6 rows

nodeadd 节点注释的：用otu表格和分组文件进行注释

nodeadd函数只是提供了简单的用注释函数，用户可以自己在node的表格后面添加各种注释信息。

Hide


tax_table = ps_net %>%
  vegan_tax() %>%
  as.data.frame()
#---node节点注释#-----------
nodes = nodeadd(plotcord =node,otu_table = otu_table,tax_table = tax_table)
nodes[1:6,1:6]

	X1 <dbl>	X2 <dbl>	elements <chr>	Kingdom <chr>	Phylum <chr>	Class <chr>
ASV_1	0.000000	10.5000000	ASV_1	Bacteria	Actinobacteria	Actinobacteria
ASV_10	3.792382	-8.9964485	ASV_10	Bacteria	Proteobacteria	Alphaproteobacteria
ASV_100	7.269286	-5.1349547	ASV_100	Bacteria	Proteobacteria	Betaproteobacteria
ASV_101	5.911236	5.0628075	ASV_101	Bacteria	Proteobacteria	Unassigned
ASV_102	4.278276	1.3951201	ASV_102	Bacteria	Proteobacteria	Gammaproteobacteria
ASV_103	-8.359017	-0.4411929	ASV_103	Bacteria	Proteobacteria	Alphaproteobacteria

6 rows

计算边

Hide

#-----计算边#--------
edge = edgeBuild(cor = cor,node = node)
head(edge)

	X2 <dbl>	Y2 <dbl>	OTU_2 <chr>	OTU_1 <chr>	weight <dbl>	X1 <dbl>	Y1 <dbl>	cor <chr>
1	2.598076	9.0000000	ASV_14	ASV_172	-0.8020673	4.278276	3.2492221	-
2	7.131446	5.3221711	ASV_28	ASV_20	0.8041345	7.755181	-0.6122717	+
3	4.533370	0.8221711	ASV_140	ASV_64	0.8484250	1.818041	-4.5620057	+
4	7.014193	-4.5620057	ASV_166	ASV_43	-0.8064018	-2.853170	8.4270510	-
5	7.014193	-4.5620057	ASV_166	ASV_172	-0.8016529	4.278276	3.2492221	-
6	1.818041	-7.5620057	ASV_29	ASV_43	-0.8415076	-2.853170	8.4270510	-

6 rows

出图

Hide

p <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = as.factor(cor)),
                                data = edge, size = 0.5) +
  geom_point(aes(X1, X2,fill = Phylum,size = mean),pch = 21, data = nodes) +
  scale_colour_brewer(palette = "Set1") +
  scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) +
  theme(panel.background = element_blank()) +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
  theme(legend.background = element_rect(colour = NA)) +
  theme(panel.background = element_rect(fill = "white",  colour = NA)) +
  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
p

Hide


# ggsave("cs1.pdf",p,width = 8,height = 6)

模拟不同的分组–可视化

模拟不同分组效果展示：3个分组

这是网络布局的基础，无论是什么聚类布局，都需要制作一个分组文件，这个文件有两列，一列是节点，一列是分组信息，这个分组信息名称必须设定为：group。

Hide


netClu = data.frame(ID = row.names(tax_table),group =rep(1,length(row.names(tax_table)))[1:length(row.names(tax_table))] )
netClu$group = as.factor(netClu$group)
result2 = PolygonClusterG (cor = cor,nodeGroup =netClu )
node = result2[[1]]
# ---node节点注释#-----------
nodes = nodeadd(plotcord =node,otu_table = otu_table,tax_table = tax_table)

#-----计算边#--------
edge = edgeBuild(cor = cor,node = node)

### 出图
pnet <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = as.factor(cor)),
                                data = edge, size = 0.5) +
  geom_point(aes(X1, X2,fill = Phylum,size = mean),pch = 21, data = nodes) +
  scale_colour_brewer(palette = "Set1") +
  scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) +
  # labs( title = paste(layout,"network",sep = "_"))+
  # geom_text_repel(aes(X1, X2,label=Phylum),size=4, data = plotcord)+
  # discard default grid + titles in ggplot2
  theme(panel.background = element_blank()) +
  # theme(legend.position = "none") +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
  theme(legend.background = element_rect(colour = NA)) +
  theme(panel.background = element_rect(fill = "white",  colour = NA)) +
  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
pnet

Hide


# ggsave("cs2.pdf",pnet,width = 7,height = 5.5)

模拟不同的分组查看效果：8个分组

Hide

netClu = data.frame(ID = row.names(cor),group =rep(1:8,length(row.names(cor)))[1:length(row.names(cor))] )
netClu$group = as.factor(netClu$group)
result2 = PolygonClusterG (cor = cor,nodeGroup =netClu )
node = result2[[1]]
# ---node节点注释#-----------
nodes = nodeadd(plotcord =node,otu_table = otu_table,tax_table = tax_table)
#-----计算边#--------
edge = edgeBuild(cor = cor,node = node)
### 出图
pnet <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = as.factor(cor)),
                                data = edge, size = 0.5) +
  geom_point(aes(X1, X2,fill = Phylum,size = mean),pch = 21, data = nodes) +
  scale_colour_brewer(palette = "Set1") +
  scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) +
  # labs( title = paste(layout,"network",sep = "_"))+
  # geom_text_repel(aes(X1, X2,label=Phylum),size=4, data = plotcord)+
  # discard default grid + titles in ggplot2
  theme(panel.background = element_blank()) +
  # theme(legend.position = "none") +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
  theme(legend.background = element_rect(colour = NA)) +
  theme(panel.background = element_rect(fill = "white",  colour = NA)) +
  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
pnet

Hide


# ggsave("cs3.pdf",pnet,width = 6,height = 5)

微生物分类分组可视化

Hide


#----------计算相关#----
result = corMicro (ps = ps,
                   N = 200,
                   method.scale = "TMM",
                   r.threshold=0.8,
                   p.threshold=0.05,
                   method = "spearman"
                   )
#--提取相关矩阵
cor = result[[1]]
# head(cor)
#-网络中包含的OTU的phyloseq文件提取
ps_net = result[[3]]
#-导出otu表格
otu_table = ps_net %>% 
  vegan_otu() %>%
  t() %>%
  as.data.frame()
tax = ps_net %>% vegan_tax() %>%
  as.data.frame()
tax$filed = tax$Phylum
group2 <- data.frame(ID = row.names(tax),group = tax$Phylum)
group2$group  =as.factor(group2$group)
result2 = PolygonClusterG (cor = cor,nodeGroup =group2 )
node = result2[[1]]
# ---node节点注释#-----------
nodes = nodeadd(plotcord =node,otu_table = otu_table,tax_table = tax_table)
#-----计算边#--------
edge = edgeBuild(cor = cor,node = node)
### 出图
pnet <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = as.factor(cor)),
                                data = edge, size = 0.5) +
  geom_point(aes(X1, X2,fill = Phylum,size = mean),pch = 21, data = nodes) +
  scale_colour_brewer(palette = "Set1") +
  scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) +
  # labs( title = paste(layout,"network",sep = "_"))+
  # geom_text_repel(aes(X1, X2,label=Phylum),size=4, data = plotcord)+
  # discard default grid + titles in ggplot2
  theme(panel.background = element_blank()) +
  # theme(legend.position = "none") +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
  theme(legend.background = element_rect(colour = NA)) +
  theme(panel.background = element_rect(fill = "white",  colour = NA)) +
  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
pnet

Hide


# ggsave("cs3.pdf",pnet,width = 6,height = 5)

微生物分类可视化布局优化1-圆环半径调整PolygonRrClusterG

结果发现这些高丰度OTU大部分属于放线菌门和变形菌门，其他比较少。所以下面我们按照OTU数量的多少，对每个模块的大小进行重新调整。

Hide


result2 = PolygonRrClusterG(cor = cor,nodeGroup =group2 )
node = result2[[1]]
# ---node节点注释#-----------
nodes = nodeadd(plotcord =node,otu_table = otu_table,tax_table = tax_table)
#-----计算边#--------
edge = edgeBuild(cor = cor,node = node)
### 出图
pnet <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = as.factor(cor)),
                                data = edge, size = 0.5) +
  geom_point(aes(X1, X2,fill = Phylum,size = mean),pch = 21, data = nodes) +
  scale_colour_brewer(palette = "Set1") +
  scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) +
  # labs( title = paste(layout,"network",sep = "_"))+
  # geom_text_repel(aes(X1, X2,label=Phylum),size=4, data = plotcord)+
  # discard default grid + titles in ggplot2
  theme(panel.background = element_blank()) +
  # theme(legend.position = "none") +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
  theme(legend.background = element_rect(colour = NA)) +
  theme(panel.background = element_rect(fill = "white",  colour = NA)) +
  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
pnet

Hide


# ggsave("cs4.pdf",pnet,width = 6,height = 5)

微生物分类可视化布局优化2-model_filled_circle

用实心点作为每个模块的布局方式

Hide


set.seed(12)
#-实心圆2
result2 = model_filled_circle(cor = cor,
                              culxy =TRUE,
                              da = NULL,# 数据框，包含x,和y列
                              nodeGroup = group2,
                              mi.size = 1,# 最小圆圈的半径，越大半径越大
                              zoom = 0.3# 不同模块之间距离
                              )
#> [1] 0
#> [1] 1
#> [1] 2
#> [1] 0
#> [1] 1
#> [1] 2
#> [1] 3
#> [1] 4

# result2 = PolygonRrClusterG (cor = cor,nodeGroup =group2 )
node = result2[[1]]
# ---node节点注释#-----------
nodes = nodeadd(plotcord =node,otu_table = otu_table,tax_table = tax_table)
#-----计算边#--------
edge = edgeBuild(cor = cor,node = node)
### 出图
pnet <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = as.factor(cor)),
                                data = edge, size = 0.5) +
  geom_point(aes(X1, X2,fill = Phylum,size = mean),pch = 21, data = nodes) +
  scale_colour_brewer(palette = "Set1") +
  scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) +
  # labs( title = paste(layout,"network",sep = "_"))+
  # geom_text_repel(aes(X1, X2,label=Phylum),size=4, data = plotcord)+
  # discard default grid + titles in ggplot2
  theme(panel.background = element_blank()) +
  # theme(legend.position = "none") +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
  theme(legend.background = element_rect(colour = NA)) +
  theme(panel.background = element_rect(fill = "white",  colour = NA)) +
  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
pnet

Hide

# ggsave("cs5.pdf",pnet,width = 6,height = 5)

微生物分类可视化布局优化3 model_maptree_group

用实心点作为每个模块布局方式2：model_maptree_group，智能布局不同分组之间的距离，在美学上特征更明显一点。

Hide


set.seed(12)
#-实心圆2
result2 = model_maptree_group(cor = cor,
                              nodeGroup = group2,
                              )

# result2 = PolygonRrClusterG (cor = cor,nodeGroup =group2 )
node = result2[[1]]
# ---node节点注释#-----------
nodes = nodeadd(plotcord =node,otu_table = otu_table,tax_table = tax_table)
#-----计算边#--------
edge = edgeBuild(cor = cor,node = node)
### 出图
pnet <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = as.factor(cor)),
                                data = edge, size = 0.5) +
  geom_point(aes(X1, X2,fill = Phylum,size = mean),pch = 21, data = nodes) +
  scale_colour_brewer(palette = "Set1") +
  scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) +
  # labs( title = paste(layout,"network",sep = "_"))+
  # geom_text_repel(aes(X1, X2,label=Phylum),size=4, data = plotcord)+
  # discard default grid + titles in ggplot2
  theme(panel.background = element_blank()) +
  # theme(legend.position = "none") +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
  theme(legend.background = element_rect(colour = NA)) +
  theme(panel.background = element_rect(fill = "white",  colour = NA)) +
  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
pnet

Hide


# ggsave("cs6.pdf",pnet,width = 6,height = 5)

按照网络模块化分组

模块布局算法 model_maptree_group

按照物种组成分类完成网络分析其实并不常用，更多的是按照模块分组，进行网络可视化。

Hide

#--modulGroup函数用于计算模块并整理成分组信息
netClu  = modulGroup( cor = cor,cut = NULL,method = "cluster_fast_greedy" )

result2 = model_maptree_group(cor = cor,
                              nodeGroup = group2,
                              )

# result2 = PolygonRrClusterG (cor = cor,nodeGroup =group2 )
node = result2[[1]]

# ---node节点注释#-----------
nodes = nodeadd(plotcord =node,otu_table = otu_table,tax_table = tax_table)
# head(nodes)

nodes2 = nodes %>% inner_join(netClu,by = c("elements" = "ID"))
nodes2$group = paste("Model_",nodes2$group,sep = "")

#-----计算边#--------
edge = edgeBuild(cor = cor,node = node)

### 出图
pnet <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = as.factor(cor)),
                                data = edge, size = 0.5) +
  geom_point(aes(X1, X2,fill = group,size = mean),pch = 21, data = nodes2) +
  scale_colour_brewer(palette = "Set1") +
  scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) +
  # labs( title = paste(layout,"network",sep = "_"))+
  # geom_text_repel(aes(X1, X2,label=Phylum),size=4, data = plotcord)+
  # discard default grid + titles in ggplot2
  theme(panel.background = element_blank()) +
  # theme(legend.position = "none") +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
  theme(legend.background = element_rect(colour = NA)) +
  theme(panel.background = element_rect(fill = "white",  colour = NA)) +
  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
pnet

Hide


# ggsave("cs7.pdf",pnet,width = 6,height = 5)

模块布局算法 model_maptree2

使用升级的model_maptree2：不在可以将每个模块独立区分，而是将模块聚拢，并在整体布局上将离散的点同这些模块一同绘制到同心圆内。控制了图形的整体布局为圆形。

Hide


result2 = model_maptree2(cor = cor,
                              method = "cluster_fast_greedy"
                              )

# result2 = PolygonRrClusterG (cor = cor,nodeGroup =group2 )
node = result2[[1]]

# ---node节点注释#-----------
nodes = nodeadd(plotcord =node,otu_table = otu_table,tax_table = tax_table)
head(nodes)

	X1 <dbl>	X2 <dbl>	elements <chr>	Kingdom <chr>	Phylum <chr>	Class <chr>
ASV_1	-13.8901865	-0.03852081	ASV_1	Bacteria	Actinobacteria	Actinobacteria
ASV_10	-0.4097123	3.33964715	ASV_10	Bacteria	Proteobacteria	Alphaproteobacteria
ASV_100	-1.1227948	-12.61245151	ASV_100	Bacteria	Proteobacteria	Betaproteobacteria
ASV_101	5.7161404	20.40310410	ASV_101	Bacteria	Proteobacteria	Unassigned
ASV_102	-18.4899174	-10.64072413	ASV_102	Bacteria	Proteobacteria	Gammaproteobacteria
ASV_103	16.2557654	10.33992920	ASV_103	Bacteria	Proteobacteria	Alphaproteobacteria

6 rows | 1-7 of 12 columns

Hide


nodes2 = nodes %>% inner_join(netClu,by = c("elements" = "ID"))
nodes2$group = paste("Model_",nodes2$group,sep = "")

#-----计算边#--------
edge = edgeBuild(cor = cor,node = node)

### 出图
pnet <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = as.factor(cor)),
                                data = edge, size = 0.5) +
  geom_point(aes(X1, X2,fill = group,size = mean),pch = 21, data = nodes2) +
  scale_colour_brewer(palette = "Set1") +
  scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) +
  # labs( title = paste(layout,"network",sep = "_"))+
  # geom_text_repel(aes(X1, X2,label=Phylum),size=4, data = plotcord)+
  # discard default grid + titles in ggplot2
  theme(panel.background = element_blank()) +
  # theme(legend.position = "none") +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
  theme(legend.background = element_rect(colour = NA)) +
  theme(panel.background = element_rect(fill = "white",  colour = NA)) +
  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
pnet

Hide



# ggsave("cs8.pdf",pnet,width = 6,height = 5)

上千OTU相关性计算测试model_igraph2布局

这个布局最近几年文章上使用非常多。

Hide

# cor_Big_micro2 增加了标准化方法和p值矫正方法
result = cor_Big_micro2(ps = ps,
                       N = 1000,
                       r.threshold=0.85,
                       p.threshold=0.05,
                       method = "pearson",
                       scale = FALSE
)

#--提取相关矩阵
cor = result[[1]]
dim(cor)
#> [1] 1000 1000

# model_igraph2
result2 <- model_igraph2(cor = cor,
                         method = "cluster_fast_greedy",
                         seed = 12
)
node = result2[[1]]
dim(node)
#> [1] 293   3


dat = result2[[2]]
head(dat)

	orig_model <dbl>	model <chr>	color <chr>	OTU <chr>	X1 <dbl>	X2 <dbl>
1	45	mini_model	#C1C1C1	ASV_1591	7.835825	12.339063
2	15	model_15	#50A9AF	ASV_688	-1.523618	-9.873582
3	34	model_34	#FEF6B1	ASV_8	4.877237	13.710152
4	6	model_6	#F99655	ASV_174	2.811909	-10.288624
5	6	model_6	#F99655	ASV_716	4.314573	-10.310663
6	34	model_34	#FEF6B1	ASV_106	3.768050	13.839837

6 rows

Hide

tem = data.frame(mod = dat$model,col = dat$color) %>%  
  dplyr::distinct( mod, .keep_all = TRUE)  
col = tem$col
names(col) = tem$mod

#---node节点注释#-----------
otu_table = as.data.frame(t(vegan_otu(ps)))
tax_table = as.data.frame(vegan_tax(ps))
nodes = nodeadd(plotcord =node,otu_table = otu_table,tax_table = tax_table)
head(nodes)

	X1 <dbl>	X2 <dbl>	elements <chr>	Kingdom <chr>	Phylum <chr>	Class <chr>
ASV_10	-2.259642	10.577060	ASV_10	Bacteria	Proteobacteria	Alphaproteobacteria
ASV_100	-7.994668	8.125110	ASV_100	Bacteria	Proteobacteria	Betaproteobacteria
ASV_102	6.730077	-13.330438	ASV_102	Bacteria	Proteobacteria	Gammaproteobacteria
ASV_104	-2.682913	12.612126	ASV_104	Bacteria	Proteobacteria	Betaproteobacteria
ASV_1043	3.463727	-9.298511	ASV_1043	Bacteria	Acidobacteria	Acidobacteria_Gp10
ASV_1048	12.307434	-7.316681	ASV_1048	Bacteria	Actinobacteria	Actinobacteria

6 rows | 1-7 of 12 columns

Hide

#-----计算边#--------
edge = edgeBuild(cor = cor,node = node)
colnames(edge)[8] = "cor"
head(edge)

	X2 <dbl>	Y2 <dbl>	OTU_2 <chr>	OTU_1 <chr>	weight <dbl>	X1 <dbl>	Y1 <dbl>	cor <chr>
1	7.835825	12.339063	ASV_1591	ASV_710	0.9060991	7.958023	13.239291	+
2	-1.523618	-9.873582	ASV_688	ASV_596	0.9262611	-2.153803	-10.584342	+
3	-1.523618	-9.873582	ASV_688	ASV_780	0.8884699	-1.365654	-10.865238	+
4	-1.523618	-9.873582	ASV_688	ASV_20	0.8557562	-1.207825	-8.872789	+
5	4.877237	13.710152	ASV_8	ASV_106	0.8734179	3.768050	13.839837	+
6	4.877237	13.710152	ASV_8	ASV_334	0.9389280	4.640294	14.523675	+

6 rows

Hide



tem2 = dat %>% 
  dplyr::select(OTU,model,color) %>%
  dplyr::right_join(edge,by =c("OTU" = "OTU_1" ) ) %>%
  dplyr::rename(OTU_1 = OTU,model1 = model,color1 = color)
head(tem2)

	OTU_1 <chr>	model1 <chr>	color1 <chr>	X2 <dbl>	Y2 <dbl>	OTU_2 <chr>	weight <dbl>	X1 <dbl>	Y1 <dbl>
1	ASV_1591	mini_model	#C1C1C1	7.958023	13.239291	ASV_710	0.9060991	7.835825	12.339063
2	ASV_688	model_15	#50A9AF	-2.153803	-10.584342	ASV_596	0.9262611	-1.523618	-9.873582
3	ASV_688	model_15	#50A9AF	-1.365654	-10.865238	ASV_780	0.8884699	-1.523618	-9.873582
4	ASV_688	model_15	#50A9AF	-1.207825	-8.872789	ASV_20	0.8557562	-1.523618	-9.873582
5	ASV_8	model_34	#FEF6B1	3.768050	13.839837	ASV_106	0.8734179	4.877237	13.710152
6	ASV_8	model_34	#FEF6B1	4.640294	14.523675	ASV_334	0.9389280	4.877237	13.710152

6 rows | 1-10 of 11 columns

Hide


tem3 = dat %>% 
  dplyr::select(OTU,model,color) %>%
  dplyr::right_join(edge,by =c("OTU" = "OTU_2" ) ) %>%
  dplyr::rename(OTU_2 = OTU,model2 = model,color2 = color)
head(tem3)

	OTU_2 <chr>	model2 <chr>	color2 <chr>	X2 <dbl>	Y2 <dbl>	OTU_1 <chr>	weight <dbl>	X1 <dbl>	Y1 <dbl>
1	ASV_1591	mini_model	#C1C1C1	7.835825	12.339063	ASV_710	0.9060991	7.958023	13.239291
2	ASV_688	model_15	#50A9AF	-1.523618	-9.873582	ASV_596	0.9262611	-2.153803	-10.584342
3	ASV_688	model_15	#50A9AF	-1.523618	-9.873582	ASV_780	0.8884699	-1.365654	-10.865238
4	ASV_688	model_15	#50A9AF	-1.523618	-9.873582	ASV_20	0.8557562	-1.207825	-8.872789
5	ASV_8	model_34	#FEF6B1	4.877237	13.710152	ASV_106	0.8734179	3.768050	13.839837
6	ASV_8	model_34	#FEF6B1	4.877237	13.710152	ASV_334	0.9389280	4.640294	14.523675

6 rows | 1-10 of 11 columns

Hide


tem4 = tem2 %>%inner_join(tem3)
head(tem4)

	OTU_1 <chr>	model1 <chr>	color1 <chr>	X2 <dbl>	Y2 <dbl>	OTU_2 <chr>	weight <dbl>	X1 <dbl>	Y1 <dbl>
1	ASV_1591	mini_model	#C1C1C1	7.958023	13.239291	ASV_710	0.9060991	7.835825	12.339063
2	ASV_688	model_15	#50A9AF	-2.153803	-10.584342	ASV_596	0.9262611	-1.523618	-9.873582
3	ASV_688	model_15	#50A9AF	-1.365654	-10.865238	ASV_780	0.8884699	-1.523618	-9.873582
4	ASV_688	model_15	#50A9AF	-1.207825	-8.872789	ASV_20	0.8557562	-1.523618	-9.873582
5	ASV_8	model_34	#FEF6B1	3.768050	13.839837	ASV_106	0.8734179	4.877237	13.710152
6	ASV_8	model_34	#FEF6B1	4.640294	14.523675	ASV_334	0.9389280	4.877237	13.710152

6 rows | 1-10 of 13 columns

Hide


edge2 = tem4 %>% mutate(color = ifelse(model1 == model2,as.character(model1),"across"),
                        manual = ifelse(model1 == model2,as.character(color1),"#C1C1C1")
                        )


col_edge = edge2 %>% dplyr::distinct(color, .keep_all = TRUE)  %>% 
  select(color,manual)
col0 = col_edge$manual
names(col0) = col_edge$color

library(ggnewscale)

p1 <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = color),
                              data = edge2, size = 1) +
  scale_colour_manual(values = col0) 

# ggsave("./cs1.pdf",p1,width = 16,height = 14)
p2 = p1 +
   new_scale_color() +
  geom_point(aes(X1, X2,color =model), data = dat,size = 4) +
  scale_colour_manual(values = col) +
  scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) +
  theme(panel.background = element_blank()) +
  theme(axis.title.x = element_blank(), axis.title.y = element_blank()) +
  theme(legend.background = element_rect(colour = NA)) +
  theme(panel.background = element_rect(fill = "white",  colour = NA)) +
  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
p2

Hide

# ggsave("./cs2.pdf",p2,width = 12,height = 10)

网络属性和节点属性

网络性质计算

22年6月升级后版本包括了16项网络属性，包括周集中老师21年NCC文章中全部属性

Hide

dat = net_properties(igraph)
head(dat)
#>                      value
#> num.edges      122.0000000
#> num.pos.edges  101.0000000
#> num.neg.edges   21.0000000
#> num.vertices    49.0000000
#> connectance      0.1037415
#> average.degree   4.9795918


# 升级后包含的网络属性更多
dat = net_properties.2(igraph,n.hub = T)
head(dat,n = 16)
#>                                               value
#> num.edges(L)                            122.0000000
#> num.pos.edges                           101.0000000
#> num.neg.edges                            21.0000000
#> num.vertices(n)                          49.0000000
#> Connectance(edge_density)                 0.1037415
#> average.degree(Average K)                 4.9795918
#> average.path.length                       2.1710403
#> diameter                                  4.7234210
#> edge.connectivity                         0.0000000
#> mean.clustering.coefficient(Average.CC)   0.5690439
#> no.clusters                               5.0000000
#> centralization.degree                     0.1879252
#> centralization.betweenness                0.1533413
#> centralization.closeness                  1.1468429
#> RM(relative.modularity)                   0.2573947
#> the.number.of.keystone.nodes              1.0000000

dat = net_properties.3(igraph,n.hub = T)
head(dat,n = 16)
#>                                               value
#> num.edges(L)                            122.0000000
#> num.pos.edges                           101.0000000
#> num.neg.edges                            21.0000000
#> num.vertices(n)                          49.0000000
#> Connectance(edge_density)                 0.1037415
#> average.degree(Average K)                 4.9795918
#> average.path.length                       2.1710403
#> diameter                                  4.7234210
#> edge.connectivity                         0.0000000
#> mean.clustering.coefficient(Average.CC)   0.5690439
#> no.clusters                               5.0000000
#> centralization.degree                     0.1879252
#> centralization.betweenness                0.1533413
#> centralization.closeness                  1.1468429
#> RM(relative.modularity)                  -8.4867310
#> the.number.of.keystone.nodes              1.0000000

# 增加了网络模块性（modularity.net ）和随机网络模块性（modularity_random ）
# dat = net_properties.4(igraph,n.hub = F)
# head(dat,n = 16)

节点性质计算

Hide

nodepro = node_properties(igraph)
head(nodepro)
#>        igraph.degree igraph.closeness igraph.betweenness igraph.cen.degree
#> ASV_28             3        0.2569444           2.333333                 3
#> ASV_48             1        0.2032967           0.000000                 1
#> ASV_34             5        0.3274336         140.216667                 5
#> ASV_18             4        0.2846154          20.916667                 4
#> ASV_20             3        0.2569444          36.000000                 3
#> ASV_44            11        0.3627451           7.874009                11

Zipi基于模块对OTU进行分类

Hide


result = cor_Big_micro2(ps = ps,
                       N = 500,
                       r.threshold=0.6,
                       p.threshold=0.05,
                       # method = "pearson",
                       scale = FALSE
)

#--提取相关矩阵
cor = result[[1]]

result4 = nodeEdge(cor = cor)
#提取变文件
edge = result4[[1]]
#--提取节点文件
node = result4[[2]]
igraph  = igraph::graph_from_data_frame(edge, directed = FALSE, vertices = node)
res = ZiPiPlot(igraph = igraph,method = "cluster_fast_greedy")
p <- res[[1]]
p

Hide

# ggsave("./cs2.pdf",p,width = 8,height = 7)

扩展-关键OTU挑选

Hub节点是在网络中与其他节点连接较多的节点，Hub微生物就是与其他微生物联系较为紧密的微生物，可以称之为关键微生物（keystone）

Hide


hub = hub_score(igraph)$vector %>%
  sort(decreasing = TRUE) %>%
  head(5) %>%
  as.data.frame()

colnames(hub) = "hub_sca"

ggplot(hub) +
  geom_bar(aes(x = hub_sca,y = reorder(row.names(hub),hub_sca)),stat = "identity",fill = "#4DAF4A")

Hide



# ggsave("./cs2.pdf",width = 5,height = 4)

对应随机网络构建和网络参数比对

Hide


result = random_Net_compate(igraph = igraph, type = "gnm", step = 100, netName = layout)
p1 = result[[1]]
sum_net = result[[4]]
p1

Hide

head(sum_net)
#>                          KO        Means SD
#> num.edges      1.078000e+03 1.078000e+03  0
#> num.pos.edges  7.570000e+02 0.000000e+00  0
#> num.neg.edges  3.210000e+02 0.000000e+00  0
#> num.vertices   3.350000e+02 3.350000e+02  0
#> connectance    1.926892e-02 1.926892e-02  0
#> average.degree 6.435821e+00 6.435821e+00  0
# ggsave("./cs3.pdf",p1,width = 5,height = 4)

微生物网络流程

微生物组小网络：model_Gephi.2

使用network函数运行微生物网络全套分析：

使用OTU数量建议少于250个，如果OTU数量为250个，同时计算zipi，整个运算过程为3-5min。

Hide


data("ps16s")
path = "./result_micro_200/"
dir.create(path)

result = network(ps = ps16s,
                 N = 200,
                layout_net = "model_Gephi.2",
                 r.threshold=0.6,
                 p.threshold=0.05,
                 label = FALSE,
                 path = path,
                 zipi = TRUE)
#> [1] "Group1"
#> [1] "cor matrix culculating over"
#> [1] "1"
#> [1] 0
#> [1] 1
#> [1] 2
#> [1] 3
#> [1] 4
#> [1] 5
#> [1] 6
#> [1] "Group2"
#> [1] "cor matrix culculating over"
#> [1] "1"
#> [1] 0
#> [1] 1
#> [1] 2
#> [1] 3
#> [1] 4
#> [1] 5
#> [1] 6
#> [1] "Group3"
#> [1] "cor matrix culculating over"
#> [1] "1"
#> [1] 0
#> [1] 1
#> [1] 2
#> [1] 3
#> [1] 4
#> [1] 5
#> [1] 6

# 多组网络绘制到一个面板
p = result[[1]]
p

Hide

# 全部样本网络参数比对
data = result[[2]]


plotname1 = paste(path,"/network_all.jpg",sep = "")
ggsave(plotname1, p,width = 48,height = 16,dpi = 72)

plotname1 = paste(path,"/network_all.pdf",sep = "")
ggsave(plotname1, p,width = 48,height = 16)


tablename <- paste(path,"/co-occurrence_Grobel_net",".csv",sep = "")
write.csv(data,tablename)

微生物大网络：model_maptree2

大网络运算时间会比较长，这里我没有计算zipi，用时5min完成全部运行。N=0，代表用全部的OTU进行计算。

3000个OTU不计算zipi全套需要18min。

Hide


path = "./result_big_1000/"
dir.create(path)

result = network.2(ps = ps16s,
                 N = 1000,
                 big = TRUE,
                 maxnode = 5,
                 select_layout = TRUE,
                  layout_net = "model_maptree2",
                 r.threshold=0.4,
                 p.threshold=0.01,
                 label = FALSE,
                 path = path,
                 zipi = FALSE)
#> [1] "Group1"
#> [1] "cor matrix culculating over"
#> [1] "Group2"
#> [1] "cor matrix culculating over"
#> [1] "Group3"
#> [1] "cor matrix culculating over"

# 多组网络绘制到一个面板
p = result[[1]]
p

Hide

# 全部样本网络参数比对
data = result[[2]]
num= 3
# plotname1 = paste(path,"/network_all.jpg",sep = "")
# ggsave(plotname1, p,width = 16*num,height = 16,dpi = 72)

plotname1 = paste(path,"/network_all.pdf",sep = "")
ggsave(plotname1, p,width = 10*num,height = 10,limitsize = FALSE)

tablename <- paste(path,"/co-occurrence_Grobel_net",".csv",sep = "")
write.csv(data,tablename)

微生物大网络：model_igraph

Hide


path = "./result_1000_igraph/"
dir.create(path)
# map = sample_data(ps)
# map$Group = "one"
# sample_data(ps16s) = map
result = network.2(ps = ps16s,
                 N = 1000,
                 big = TRUE,
                 maxnode = 5,
                 select_layout = TRUE,
                  layout_net = "model_igraph",
                 r.threshold=0.4,
                 p.threshold=0.01,
                 label = FALSE,
                 path = path,
                 zipi = FALSE)
#> [1] "Group1"
#> [1] "cor matrix culculating over"
#> [1] "Group2"
#> [1] "cor matrix culculating over"
#> [1] "Group3"
#> [1] "cor matrix culculating over"

# 多组网络绘制到一个面板
p = result[[1]]
# 全部样本网络参数比对
data = result[[2]]
num= 3
# plotname1 = paste(path,"/network_all.jpg",sep = "")
# ggsave(plotname1, p,width = 16*num,height = 16,dpi = 72)

plotname1 = paste(path,"/network_all.pdf",sep = "")
ggsave(plotname1, p,width = 16*num,height = 16,limitsize = FALSE)

tablename <- paste(path,"/co-occurrence_Grobel_net",".csv",sep = "")
write.csv(data,tablename)

微生物大网络：model_igraph2-network.i

model_igraph2在model_igraph的基础上去除了离散点，这对于重复数量少的研究比较有利，可以更加清楚的展示小样本网络。

network.i函数专门为model_igraph2算法而写，可以额外输出模块上色的网络。

Hide



path = "./result_1000_igraph2/"
dir.create(path)

# map = sample_data(ps)
# map$Group = "one"
# sample_data(ps16s) = map

result = network.i(ps =  ps16s,
                   N = 1000,
                   r.threshold=0.9,
                   big = T,
                   select_layout = T,
                   method = "pearson",
                   scale = FALSE,
                   layout_net = "model_igraph2",
                   p.threshold=0.05,
                   label = FALSE,
                   path =path ,
                   zipi = F,
                   order = NULL
                   )
#> [1] "Group1"
#> [1] "cor matrix culculating over"
#> [1] "Group2"
#> [1] "cor matrix culculating over"
#> [1] "Group3"
#> [1] "cor matrix culculating over"

p1 = result[[1]]
p1

Hide

dat = result[[2]]
tablename <- paste(path,"/co-occurrence_Grobel_net",".csv",sep = "")
write.csv(dat,tablename)
p = result[[5]]


plotname1 = paste(path,"/network_all.pdf",sep = "")
ggsave(plotname1, p1,width = 10,height = 3,limitsize = FALSE)

plotname1 = paste(path,"/network_all2.pdf",sep = "")
ggsave(plotname1, p,width = 25,height = 8,limitsize = FALSE)

跨域网络流程

细菌，真菌，环境因子三者跨域网络

Hide

#--细菌和环境因子的网络#--------
Envnetplot<- paste("./16s_ITS_Env_network",sep = "")
dir.create(Envnetplot)

ps16s =ps16s%>% ggClusterNet::scale_micro()
psITS = psITS%>% ggClusterNet::scale_micro()


#--细菌和真菌ps对象中的map文件要一样
ps.merge <- ggClusterNet::merge16S_ITS(ps16s = ps16s,
                                       psITS= psITS,
                         NITS = 200,
                         N16s = 200)

map =  phyloseq::sample_data(ps.merge)
#  
map$Group = "one"
phyloseq::sample_data(ps.merge) <- map

#--环境因子导入
data1 = env
envRDA  = data.frame(row.names = env$ID,env[,-1])

envRDA.s = vegan::decostand(envRDA,"hellinger")
data1[,-1] = envRDA.s

Gru = data.frame(ID = colnames(env)[-1],group = "env" )
# head(Gru)

library(phyloseq)
library(sna)
library(ggClusterNet)
library(igraph)
library(tidyverse)

result <- ggClusterNet::corBionetwork(ps = ps.merge,
                        N = 0,
                        r.threshold = 0.4, # 相关阈值
                        p.threshold = 0.05,
                        big = T,
                        group = "Group",
                        env = data1, # 环境指标表格
                        envGroup = Gru,# 环境因子分组文件表格
                        # layout = "fruchtermanreingold",
                        path = Envnetplot,# 结果文件存储路径
                        fill = "Phylum", # 出图点填充颜色用什么值
                        size = "igraph.degree", # 出图点大小用什么数据
                        scale = TRUE, # 是否要进行相对丰度标准化
                        bio = TRUE, # 是否做二分网络
                        zipi = F, # 是否计算ZIPI
                        step = 100, # 随机网络抽样的次数
                        width = 18,
                        label = TRUE,
                        height = 10
)
#> [1] "one"
#>  num [1:34, 1:18] 0.0474 0.0633 0.0188 0.0588 0.069 ...
#>  - attr(*, "dimnames")=List of 2
#>   ..$ : chr [1:34] "pH" "SOC" "TN" "NH4.N" ...
#>   ..$ : chr [1:18] "sample1" "sample10" "sample11" "sample12" ...
#> [1] "1"
#> [1] "2"
#> [1] "3"


p = result[[1]]
p

Hide

# 全部样本网络参数比对
data = result[[2]]
plotname1 = paste(Envnetplot,"/network_all.jpg",sep = "")
ggsave(plotname1, p,width = 15,height = 12,dpi = 72)
# plotname1 = paste(Envnetplot,"/network_all.png",sep = "")
# ggsave(plotname1, p,width = 10,height = 8,dpi = 72)
plotname1 = paste(Envnetplot,"/network_all.pdf",sep = "")
ggsave(plotname1, p,width = 15,height = 12)
tablename <- paste(Envnetplot,"/co-occurrence_Grobel_net",".csv",sep = "")
write.csv(data,tablename)

细菌-环境因子跨域网络

Hide


#--细菌-环境因子网络#-----
Envnetplot<- paste("./16s_Env_network",sep = "")
dir.create(Envnetplot)
# ps16s = ps0 %>% ggClusterNet::scale_micro()
psITS = NULL

#--细菌和真菌ps对象中的map文件要一样
ps.merge <- ggClusterNet::merge16S_ITS(ps16s = ps16s,
                                       psITS = NULL,
                                       N16s = 500)
ps.merge
#> phyloseq-class experiment-level object
#> otu_table()   OTU Table:         [ 500 taxa and 18 samples ]
#> sample_data() Sample Data:       [ 18 samples by 2 sample variables ]
#> tax_table()   Taxonomy Table:    [ 500 taxa by 8 taxonomic ranks ]
map =  phyloseq::sample_data(ps.merge)
#  
map$Group = "one"
phyloseq::sample_data(ps.merge) <- map
data(env)
data1 = env
envRDA  = data.frame(row.names = env$ID,env[,-1])

envRDA.s = vegan::decostand(envRDA,"hellinger")
data1[,-1] = envRDA.s

Gru = data.frame(ID = colnames(env)[-1],group = "env" )
# head(Gru)


# library(sna)
# library(ggClusterNet)
# library(igraph)

result <- ggClusterNet::corBionetwork(ps = ps.merge,
                                      N = 0,
                                      r.threshold = 0.4, # 相关阈值
                                      p.threshold = 0.05,
                                      big = T,
                                      group = "Group",
                                      env = data1, # 环境指标表格
                                      envGroup = Gru,# 环境因子分组文件表格
                                      # layout = "fruchtermanreingold",
                                      path = Envnetplot,# 结果文件存储路径
                                      fill = "Phylum", # 出图点填充颜色用什么值
                                      size = "igraph.degree", # 出图点大小用什么数据
                                      scale = TRUE, # 是否要进行相对丰度标准化
                                      bio = TRUE, # 是否做二分网络
                                      zipi = F, # 是否计算ZIPI
                                      step = 100, # 随机网络抽样的次数
                                      width = 18,
                                      label = TRUE,
                                      height = 10
)
#> [1] "one"
#>  num [1:34, 1:18] 0.0474 0.0633 0.0188 0.0588 0.069 ...
#>  - attr(*, "dimnames")=List of 2
#>   ..$ : chr [1:34] "pH" "SOC" "TN" "NH4.N" ...
#>   ..$ : chr [1:18] "sample1" "sample10" "sample11" "sample12" ...
#> [1] "1"
#> [1] "2"
#> [1] "3"


p = result[[1]]
p

Hide

# 全部样本网络参数比对
data = result[[2]]
plotname1 = paste(Envnetplot,"/network_all.jpg",sep = "")
ggsave(plotname1, p,width = 13,height = 12,dpi = 72)

plotname1 = paste(Envnetplot,"/network_all.pdf",sep = "")
ggsave(plotname1, p,width = 13,height = 12)
tablename <- paste(Envnetplot,"/co-occurrence_Grobel_net",".csv",sep = "")
write.csv(data,tablename)

细菌和真菌跨域网络

Hide

#---细菌和真菌OTU网络-域网络-二分网络#-------
# 仅仅关注细菌和真菌之间的相关，不关注细菌内部和真菌内部相关

Envnetplot<- paste("./16S_ITS_network",sep = "")
dir.create(Envnetplot)


data(psITS)
#--细菌和真菌ps对象中的map文件要一样
ps.merge <- ggClusterNet::merge16S_ITS(ps16s = ps16s,
                                       psITS = psITS,
                                       N16s = 500,
                                       NITS = 500
)

ps.merge
#> phyloseq-class experiment-level object
#> otu_table()   OTU Table:         [ 1000 taxa and 18 samples ]
#> sample_data() Sample Data:       [ 18 samples by 2 sample variables ]
#> tax_table()   Taxonomy Table:    [ 1000 taxa by 8 taxonomic ranks ]

map =  phyloseq::sample_data(ps.merge)

#  
map$Group = "one"
phyloseq::sample_data(ps.merge) <- map


data =  data.frame(ID = c("fun_ASV_205","fun_ASV_316","fun_ASV_118"),
                   c("Coremode","Coremode","Coremode"))

# source("F:/Shared_Folder/Function_local/R_function/my_R_packages/ggClusterNet/R/corBionetwork2.R")


result <- corBionetwork(ps = ps.merge,
                        N = 0,
                        lab = data,
                        r.threshold = 0.8, # 相关阈值
                        p.threshold = 0.05,
                        group = "Group",
                        # env = data1, # 环境指标表格
                        # envGroup = Gru,# 环境因子分组文件表格
                        # layout = "fruchtermanreingold",
                        path = Envnetplot,# 结果文件存储路径
                        fill = "Phylum", # 出图点填充颜色用什么值
                        size = "igraph.degree", # 出图点大小用什么数据
                        scale = TRUE, # 是否要进行相对丰度标准化
                        bio = TRUE, # 是否做二分网络
                        zipi = F, # 是否计算ZIPI
                        step = 100, # 随机网络抽样的次数
                        width = 12,
                        label = TRUE,
                        height = 10,
                        big = TRUE,
                        select_layout = TRUE,
                        layout_net = "model_maptree2",
                        clu_method = "cluster_fast_greedy"
                        
                        
)
#> [1] "one"
#> [1] "1"
#> [1] "2"
#> [1] "3"


tem <- model_maptree(cor =result[[5]],
                     method = "cluster_fast_greedy",
                     seed = 12
)
node_model = tem[[2]]
# head(node_model)

p = result[[1]]
p

Hide

# 全部样本网络参数比对
data = result[[2]]

plotname1 = paste(Envnetplot,"/network_all.pdf",sep = "")
ggsave(plotname1, p,width = 20,height = 19)
tablename <- paste(Envnetplot,"/co-occurrence_Grobel_net",".csv",sep = "")
write.csv(data,tablename)
tablename <- paste(Envnetplot,"/node_model_imformation",".csv",sep = "")
write.csv(node_model,tablename)

tablename <- paste(Envnetplot,"/nodeG_plot",".csv",sep = "")
write.csv(result[[4]],tablename)
tablename <- paste(Envnetplot,"/edge_plot",".csv",sep = "")
write.csv(result[[3]],tablename)
tablename <- paste(Envnetplot,"/cor_matrix",".csv",sep = "")
write.csv(result[[5]],tablename)

细菌真菌任意分类水平跨域网络

Hide


#--细菌真菌任意分类水平跨域网络#----------
library(tidyverse)
# res1path = "result_and_plot/Micro_and_other_index_16s/"
Envnetplot<- paste("./16S_ITS_network_Genus",sep = "")
dir.create(Envnetplot)
#--细菌和真菌ps对象中的map文件要一样
ps.merge <- merge16S_ITS(ps16s = ps16s,
                         psITS = psITS,
                         N16s = 500,
                         NITS = 500
)
ps.merge
#> phyloseq-class experiment-level object
#> otu_table()   OTU Table:         [ 1000 taxa and 18 samples ]
#> sample_data() Sample Data:       [ 18 samples by 2 sample variables ]
#> tax_table()   Taxonomy Table:    [ 1000 taxa by 8 taxonomic ranks ]
map =  phyloseq::sample_data(ps.merge)
#  
map$Group = "one"
phyloseq::sample_data(ps.merge) <- map
tem.0 = ps.merge %>% tax_glom_wt(ranks = "Genus")
tax = tem.0 %>% vegan_tax() %>%
  as.data.frame()
# head(tax)

data =  data.frame(ID = c("Acaulium","Acidocella","Acrobeloides"),
                   c("Acaulium","Acidocella","Acrobeloides"))

result <- corBionetwork(ps = tem.0,
                        N = 0,
                        lab = data,
                        r.threshold = 0.6, # 相关阈值
                        p.threshold = 0.05,
                        group = "Group",
                        # env = data1, # 环境指标表格
                        # envGroup = Gru,# 环境因子分组文件表格
                        # layout = "fruchtermanreingold",
                        path = Envnetplot,# 结果文件存储路径
                        fill = "Phylum", # 出图点填充颜色用什么值
                        size = "igraph.degree", # 出图点大小用什么数据
                        scale = TRUE, # 是否要进行相对丰度标准化
                        bio = TRUE, # 是否做二分网络
                        zipi = F, # 是否计算ZIPI
                        step = 100, # 随机网络抽样的次数
                        width = 12,
                        label = TRUE,
                        height = 10,
                        big = TRUE,
                        select_layout = TRUE,
                        layout_net = "model_maptree2",
                        clu_method = "cluster_fast_greedy"
                        
                        
)
#> [1] "one"
#> [1] "1"
#> [1] "2"
#> [1] "3"

tem <- model_maptree(cor =result[[5]],
                     method = "cluster_fast_greedy",
                     seed = 12
)
node_model = tem[[2]]
head(node_model)

	ID <chr>	group <dbl>	degree <dbl>
1	Microascus	12	13
2	Scopulariopsis	13	12
3	Arcopilus	10	10
4	Cladosporium	13	9
5	Conocybe	13	8
6	Hydrogonium	6	8

6 rows

Hide


library(WGCNA)
otu = tem.0 %>% vegan_otu() %>%
  as.data.frame()
node_model = node_model[match(colnames(otu),node_model$ID),]


MEList = moduleEigengenes(otu, colors = node_model$group)
MEs = MEList$eigengenes

nGenes = ncol(otu)
nSamples = nrow(otu)
moduleTraitCor = cor(MEs, envRDA, use = "p")
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)

#sizeGrWindow(10,6)
pdf(file=paste(Envnetplot,"/","Module-env_relationships.pdf",sep = ""),width=10,height=6)
# Will display correlations and their p-values
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
                   signif(moduleTraitPvalue, 1), ")", sep = "")

dim(textMatrix) = dim(moduleTraitCor)
par(mar = c(6, 8.5, 3, 3))
# Display the correlation values within a heatmap plot
labeledHeatmap(Matrix = moduleTraitCor,
               xLabels = names(envRDA),
               yLabels = names(MEs),
               ySymbols = names(MEs),
               colorLabels = FALSE,
               colors = greenWhiteRed(50),
               textMatrix = textMatrix,
               setStdMargins = FALSE,
               cex.text = 0.5,
               zlim = c(-1,1),
               main = paste("Module-trait relationships"))
dev.off()
#> png 
#>   2


p = result[[1]]
p

Hide

# 全部样本网络参数比对
data = result[[2]]


plotname1 = paste(Envnetplot,"/network_all.pdf",sep = "")
ggsave(plotname1, p,width = 10,height = 10)
tablename <- paste(Envnetplot,"/co-occurrence_Grobel_net",".csv",sep = "")
write.csv(data,tablename)


tablename <- paste(Envnetplot,"/nodeG_plot",".csv",sep = "")
write.csv(result[[4]],tablename)
tablename <- paste(Envnetplot,"/edge_plot",".csv",sep = "")
write.csv(result[[3]],tablename)
tablename <- paste(Envnetplot,"/cor_matrix",".csv",sep = "")
write.csv(result[[5]],tablename)

网络稳定性（抗干扰性）

模块相似度

Hide

#--网络稳定性:模块相似性------
data(ps)

library(tidyfst)

res = module.compare.m(
    ps = ps,
    Top = 500,
    degree = TRUE,
    zipi = FALSE,
    r.threshold= 0.8,
    p.threshold=0.05,
    method = "spearman",
    padj = F,
    n = 3)
#> [1] 80  8
#> [1] 159   8
#> [1] 224   8
#> [1] 0
#> [1] 1
#> [1] 2
#> [1] 3
#> [1] 4
#> [1] 0
#> [1] 1
#> [1] 2
#> [1] 3
#> [1] 4
#> [1] 0
#> [1] 1
#> [1] 2
#> [1] 3
#> [1] 4
#不同分组使用一个圆圈展示，圆圈内一个点代表一个模块，相连接的模块代表了相似的模块。
p = res[[1]]
p

Hide

 #--提取模块的OTU，分组等的对应信息
dat = res[[2]]
head(dat)

	ID <chr>	group <chr>	Group <fct>
1	ASV_256	KOmodel_1	KO
2	ASV_142	KOmodel_1	KO
3	ASV_362	KOmodel_1	KO
4	ASV_89	KOmodel_1	KO
5	ASV_419	KOmodel_1	KO
6	ASV_352	KOmodel_1	KO

6 rows

Hide


#模块相似度结果表格
dat2 = res[[3]]
head(dat2)

	module1 <chr>	module2 <chr>	Both <chr>	P1A2 <chr>	P2A1 <chr>	A1A2 <chr>	p_raw <chr>
KOmodel_41-OEmodel_2	KOmodel_41	OEmodel_2	1	2	6	500	0.0407716775257314
KOmodel_62-OEmodel_1	KOmodel_62	OEmodel_1	1	2	6	500	0.0407716775257314
KOmodel_34-OEmodel_1	KOmodel_34	OEmodel_1	1	2	6	500	0.0407716775257314
KOmodel_22-OEmodel_4	KOmodel_22	OEmodel_4	1	3	5	500	0.0464588019672787
KOmodel_26-OEmodel_8	KOmodel_26	OEmodel_8	1	3	4	501	0.0388304723830171
KOmodel_17-OEmodel_10	KOmodel_17	OEmodel_10	1	4	4	500	0.0483470023594227

6 rows | 1-8 of 9 columns

网络鲁棒性（随机去除节点）

这里通过随机去除部分OTU，计算网络鲁棒性，代表网络抗干扰的能力。按照0.05的步长，每次去除5%的文生物，重新计算鲁棒性，知道最终全部去除。如果有分组列Group，则会按照分组进行鲁棒性计算，并且区分颜色绘制到一个面板中。计算鲁棒性这里使用丰度加成权重和不加权两种方式，左边是不加权，后侧是加权的结果。这里步长不可以修改，因为修改好像也没什么意思。

Hide


#--随即取出任意比例节点-网络鲁棒性#---------
res = Robustness.Random.removal(ps = ps,
                                Top = 500,
                                r.threshold= 0.8,
                                p.threshold=0.05,
                                method = "spearman"
                                
                                )


p = res[[1]]

p

Hide

#提取数据
dat = res[[2]]
head(dat)

	Proportion.removed <dbl>	remain.mean <dbl>	remain.sd <dbl>	remain.se <dbl>	weighted <chr>	Group <chr>
1	0.05	0.5408197	0.01310521	0.001310521	weighted	KO
2	0.10	0.5030601	0.01536428	0.001536428	weighted	KO
3	0.15	0.4674044	0.01755406	0.001755406	weighted	KO
4	0.20	0.4325683	0.01966370	0.001966370	weighted	KO
5	0.25	0.3944809	0.02034929	0.002034929	weighted	KO
6	0.30	0.3542623	0.02243074	0.002243074	weighted	KO

6 rows

网络鲁棒性（去除关键节点）

Hide

#---去除关键节点-网络鲁棒性#------
res= Robustness.Targeted.removal(ps = ps,
                                 Top = 500,
                                 degree = TRUE,
                                 zipi = FALSE,
                                 r.threshold= 0.8,
                                 p.threshold=0.05,
                                 method = "spearman")
p = res[[1]]
p

Hide

#提取数据
dat = res[[2]]

负相关比例

Hide


#--计算负相关的比例#----
res = negative.correlation.ratio(ps = ps,
                          Top = 500,
                          degree = TRUE,
                          zipi = FALSE,
                          r.threshold= 0.8,
                          p.threshold=0.05,
                          method = "spearman")

p = res[[1]]
p

Hide

dat = res[[2]]
#-负相关比例数据
head(dat)

<fct>

ratio

<dbl>

41.15139

34.19023

39.24051

3 rows

Hide

# dir.create("./negative_correlation_ratio/")
# write.csv(dat,
#           paste("./negative_correlation_ratio/","_negative_ratio_network.csv",sep = ""))

组成稳定性

Hide


treat = ps %>% sample_data()
treat$pair = paste( "A",c(rep(1:6,3)),sep = "")
# head(treat)
sample_data(ps) = treat
#一般性的处理，没有时间梯度的，这里设定time为F，意味着每两个群落的组合都进行比较
res = community.stability( ps = ps,time = FALSE)
p = res[[1]]
p

Hide

dat = res[[2]]

# 如果是时间序列，会按照时间的单一流动性方向进行比较，自然就不是两两比对了。
res = community.stability( ps = ps,time = TRUE)
p = res[[1]]
p

网络抗毁性

网络抗毁性:使用了自然连通度的概念，用于反映网络稳定性.具体而言，通过在构建好的网络中移除里面的节点，并计算自然连通度。这样随着取出点的数量逐渐增多，就可以计算得到一连串的网络连通度。这个算法最先来自PNAS的文章：“Reduced microbial stability in the active layer is associated with carbon loss under alpine permafrost degradation”

Hide


library(tidyfst)
library("pulsar")
library(ggClusterNet)
library(phyloseq)
library(tidyverse)

res = natural.con.microp (
    ps = ps,
    Top = 200,
    r.threshold= 0.5,
    p.threshold=0.05,
    method = "spearman",
    norm = F,
    end = 150,# 小于网络包含的节点数量
    start = 0,
    con.method = "pulsar"
    )

p = res[[1]]
p

Hide

dat  = res[[2]]

网络模块化

展示网络模块

这里我们通过指定模块

Hide

library(ggClusterNet)
library(igraph)
library(ggClusterNet)
library(tidyverse)
library(phyloseq)
data(ps)


resu  = module_display.2(
  pst = ps,
  r.threshold= 0.6,
  p.threshold=0.05,
  select.mod = c("model_1","model_2","model_3","model_4"),#选择指定模块可视化
  Top = 500,
  num = 5, # 模块包含OTU数量少于5个的不展示,
  leg.col = 3
)

# 全部模块输出展示
p1 = resu[[1]]
p1

Hide


#--提取率相关矩阵
dat0 = resu[[4]]

#--提取模块信息：这里区分好不同模块的OTU，以及每个模块包含的网络边的数量
dat = resu[[5]]
head(dat)

	ID <chr>	group <chr>	degree <dbl>
1	ASV_66	model_3	46
2	ASV_294	model_3	40
3	ASV_326	model_3	34
4	ASV_60	model_2	31
5	ASV_2	model_3	31
6	ASV_163	model_3	30

6 rows

Hide


#--提取网络的边和节点表格
dat2 = resu[[6]]
names(dat2)
#> [1] "edge" "node"

指定的目标模块展示

Hide

p2 = resu[[2]]
p2

Hide

#按照模块着色，模块之间的边着色为灰色展示整个网络。
p2 = resu[[3]]
p2

Hide

#--注意这个报错信息，为展示的空间不够大，拉一拉展示框即可
# Error in grid.Call(C_convert, x, as.integer(whatfrom), as.integer(whatto),  :
#   Viewport has zero dimension(s)

基于模块的微生物组成分析

对微生物数据分好模块，使用模块分类的数据导入module_composition函数，即可作指定模块的物种组成。

Hide


#--这里我们指定三个模块，绘制这三个模块物种组成图表
select.mod = c("model_1","model_2","model_3")
#--模块信息，共有三列，第一列是OTU，第二列是吗模块的分类，第三列是模块捏的边数量
mod1 = resu$mod.groups %>% filter(group %in% select.mod)
head(mod1)

	ID <chr>	group <chr>	degree <dbl>
1	ASV_66	model_3	46
2	ASV_294	model_3	40
3	ASV_326	model_3	34
4	ASV_60	model_2	31
5	ASV_2	model_3	31
6	ASV_163	model_3	30

6 rows

Hide


#-计算模块物种组成信息
pst = ps %>%
  filter_taxa(function(x) sum(x ) > 0, TRUE) %>%
  scale_micro("rela") %>%
  # subset_samples.wt("Group" ,id3[m]) %>%
  filter_OTU_ps(500)

#-选择模块的微生物组成，通过参数j选择不同的分类水平。
res = module_composition(pst = pst,mod1 = mod1,j = "Family")
p3 = res[[1]]
p3

Hide

#按照属水平进行绘制
res = module_composition(pst = pst,mod1 = mod1,j = "Genus")
p3 = res[[1]]
p3

Hide


#---提取出图数据物种组成数据导出
ps.t = res[[3]]

otu = ps.t %>% vegan_otu() %>% t() %>%
  as.data.frame()
# head(otu)
tax = ps.t %>% vegan_tax() %>%
  as.data.frame()
# head(tax)
tab.res = cbind(otu,tax)
# head(tab.res)

# 没有标准化的出图数据
tab.res.2 = res[[4]]$bundance
# head(tab.res.2)
# 标准化到100%的出图数据
tab.res.3 = res[[4]]$relaabundance
# head(tab.res.3)

模块计算alpha多样性

这里将指定模块和ps对象输入函数module_alpha即可计算三种alpha多样性，这里调用了EasyStat包，使用非参数检验检测了alpha多样性的显著性。

注意：这里的ps对象需要原始count数据，不可用标准化后的数据。

Hide



#--这里我们指定三个模块，绘制这三个模块物种组成图表
select.mod = c("model_1","model_2","model_3")

#--模块信息，两列，第一列是OTU，第二列是模块的分类
mod1 = resu$mod.groups %>% 
  dplyr::select(ID,group) %>%
  dplyr::filter(group %in% select.mod)
head(mod1)

	ID <chr>	group <chr>
1	ASV_66	model_3
2	ASV_294	model_3
3	ASV_326	model_3
4	ASV_60	model_2
5	ASV_2	model_3
6	ASV_163	model_3

6 rows

Hide


#-选择模块的多样性
result = module_alpha (
  ps = ps,
  mod1 = mod1)

p5 = result[[1]]
p5

Hide


#--alpha多样性数据
plotd = result[[4]]$alpha
# alpha多样性显著性表格
sigtab = result[[4]]$sigtab

基于单个样本计算网络属性

这里我们需要输入微生物组数据的ps对象，其次，需要输入这组微生物的相关矩阵；然后netproperties.sample函数会按照单个样本计算网络属性；这里一共有15中网络属性。

这里计算出来后便可以和其他指标进行相关性分析；

Hide


library(ggClusterNet)
library(igraph)
library(ggClusterNet)
library(tidyverse)
library(phyloseq)
data(ps)
#-计算模块物种组成信息
pst = ps %>%
  filter_taxa(function(x) sum(x ) > 0, TRUE) %>%
  scale_micro("rela") %>%
  # subset_samples.wt("Group" ,id3[m]) %>%
  filter_OTU_ps(500)

# 选择和网络属性做相关的指标
select.env = "env1"
#--内置数据无需导入

env = env1 %>% 
  as.data.frame() %>%
  rownames_to_column("id")
dat.f = netproperties.sample(pst = pst,cor = cor)
# head(dat.f)


dat.f$id = row.names(dat.f)
dat.f = dat.f %>% dplyr:: select(id,everything())


subenv = env %>% dplyr::select(id,everything()) %>% dplyr::select(id,select.env )
tab = dat.f %>% left_join(subenv,by = "id")

mtcars2 = reshape2::melt(tab, id.vars=c(select.env,"id"))
lab = mean(mtcars2[,select.env])
preoptab = tab
  
p0_1 = ggplot2::ggplot(mtcars2,aes(x= value,!!sym(select.env), colour=variable)) +
    ggplot2::geom_point() +
    ggpubr::stat_cor(label.y=lab*1.1)+
    ggpubr::stat_regline_equation(label.y=lab*1.1,vjust = 2) +
    facet_wrap(~variable, scales="free_x") +
    geom_smooth(aes(value,!!sym(select.env), colour=variable), method=lm, se=T)+
    theme_classic()

p0_1

计算模块的Zscore值

按照样本计算模块的ZS值，方便施用这一指标和其他相关指标进行相关性分析。

Hide




dat = ZS.net.module(
    pst = ps,
    Top = 500,
    r.threshold= 0.8,
    p.threshold=0.05,
    select.mod =  NULL
    )
# head(dat)


map =sample_data(ps)
map$id = row.names(map)
map = map[,c("id","Group")]
data = map %>%
    as.tibble() %>%
    inner_join(dat,by = "id") %>%
    dplyr::rename(group = Group)

colnames(env)[1] = "id"
subenv = env %>% dplyr::select(id,everything()) %>% dplyr::select(id,select.env )
# head(data)
tab = data %>% left_join(subenv,by = "id")
modenv = tab
mtcars2 = reshape2::melt(tab, id.vars=c(select.env,"group","id"))
# mtcars2$variable
head(mtcars2)

	env1 <dbl>	group <fct>	id <chr>	variable <fct>	value <dbl>
1	7.67	KO	KO1	model_1	7.1549454
2	7.61	KO	KO2	model_1	-0.3726327
3	7.15	KO	KO3	model_1	10.6994573
4	6.89	KO	KO4	model_1	10.3945282
5	7.39	KO	KO5	model_1	35.8009267
6	7.01	KO	KO6	model_1	9.5297211

6 rows

Hide

lab = mean(mtcars2[,select.env])
p1_1 = ggplot2::ggplot(mtcars2,aes(x= value,!!sym(select.env), colour=variable)) +
    ggplot2::geom_point() +
    ggpubr::stat_cor(label.y=lab*1.1)+
    ggpubr::stat_regline_equation(label.y=lab*1.1,vjust = 2) +
    facet_wrap(~variable, scales="free_x") +
    geom_smooth(aes(value,!!sym(select.env), colour=variable), method=lm, se=T)+
    theme_classic()

p1_1

模块ZS值，网络属性，环境因子相关联合到一起

Hide


#网络模块和环境因子相关-网络指标和环境相关
dat.1 = env1 %>% rownames_to_column("id")
res = net.property.module.env (
  pst = pst,
  Top = 500,
  r.threshold= 0.8,
  p.threshold=0.05,
  env = dat.1,
  select.mod = c("model_1","model_2","model_3"),
  select.env = "env1")

#--模块和环境因子之间相关
p9 = res[[1]]
p9

Hide


#-网络属性和环境因子相关
p9 = res[[2]]
p9

时空组:网络相关分析

时空组：多组网络展示-网络性质稳定性等

这部分更新与2023年1月27日完成，可能会存在bug等，如有问题即使告知我解决；

这部分更新主要是为了解决一下几个问题：- 多个网络在一张图上需要调整映射颜色相同 - 分组并不是只有一列，很多时候有时间序列，有空间部位等，优化参数使得时间序列等多个分组的研究展示更加顺畅；

Hide


rm(list=ls())
library(tidyverse)
library(ggClusterNet)
library(phyloseq)
library(igraph)
library(tidyfst)

ps.st = readRDS("./ps_TS.rds")
ps.st
#> phyloseq-class experiment-level object
#> otu_table()   OTU Table:         [ 2432 taxa and 108 samples ]
#> sample_data() Sample Data:       [ 108 samples by 4 sample variables ]
#> tax_table()   Taxonomy Table:    [ 2432 taxa by 7 taxonomic ranks ]
#> phy_tree()    Phylogenetic Tree: [ 2432 tips and 2431 internal nodes ]
#> refseq()      DNAStringSet:      [ 2432 reference sequences ]

#时空组网络-分面网络图-解决填充颜色不一致问题#----

res = Facet.network (
    ps.st= ps.st,# phyloseq对象
    g1 = "Group",# 分组1
    g2 = "space",# 分组2
    g3 = "time",# 分组3
    ord.g1 = c("WT","KO","OE"),# 排序顺序
    ord.g2 = c("B","R") ,# 排序顺序
    ord.g3 = c("T1","T2","T3") ,# 排序顺序
    order = "time", # 出图每行代表的变量
    fill = "Phylum",
    size = "igraph.degree",
    layout_net = "model_maptree2",
    r.threshold=0.8,
    p.threshold=0.01,
    method = "spearman",
    select_layout = TRUE,
    clu_method = "cluster_fast_greedy",
    maxnode = 5
)
#> [1] "B" "R"
#> [1] "T1" "T2" "T3"
#> [1] "WT" "KO" "OE"

p = res[[1]]

p

Hide


# ggsave("cs1.pdf",p,width = 5*5,height = 4*3)

#--时空组网络-网络性质#-------
dat = net.pro.ts(dat = res[[2]][[1]])
# head(dat)

Hide


#--时空组的网络稳定性
#--时空组-稳定性1模块比对#--------
res = module.compare.m.ts(ps.st = ps.st, N = 200,
  degree = TRUE,zipi = FALSE,r.threshold= 0.8,
  p.threshold=0.05,method = "spearman",
  padj = F,n = 3,
  g1 = "Group",# 分组1
  g2 = "space",# 分组2
  g3 = "time",# 分组3
  zoom = 0.3,# 控制小圈大小
  b.x = 1.3,
  b.y = 1.3)


p = res[[1]]
p

Hide

#提取数据
dat = res[[2]]
#-作图数据提取
dat = res[[3]]

#时空组-稳定性2 鲁棒性随即取出#-----
res = Robustness.Random.removal.ts(
  ps.st = ps.st,
  N = 200,
  r.threshold= 0.8,
  p.threshold=0.05,
  method = "spearman",
  order = "time",
  g1 = "Group",
  g2 = "space",
  g3 = "time")

res[[1]][[1]]

Hide

res[[1]][[2]]

Hide



#时空组-稳定性2-2 鲁棒性随即取出#-----
res = Robustness.Targeted.removal.ts(
  ps.st,
  N = 200,
  degree = TRUE,
  zipi = FALSE,
  r.threshold= 0.8,
  p.threshold=0.05,
  method = "spearman",
  order = "time",
  g1 = "Group",# 分组1
  g2 = "space",# 分组2
  g3 = "time")

res[[1]][[1]]

Hide

res[[1]][[2]]

Hide


#时空组-稳定性3负相关比例#-------
res = negative.correlation.ratio.ts(
  ps.st,
  N = 200,
  r.threshold= 0.8,
  p.threshold=0.05,
  method = "spearman",
  order = "time",
  g1 = "Group",# 分组1
  g2 = "space",# 分组2
  g3 = "time",# 分组3
  ord.g1 =NULL,# 排序顺序
  ord.g2 = NULL, # 排序顺序
  ord.g3 = NULL# 排序顺序
)
# 导出图片
res[[1]]

Hide

# 导出数据
dat = res[[2]]
head(dat)

	ID <chr>	ratio <dbl>	group <chr>	time <chr>	space <chr>	label <chr>	Group <fct>
1	R.T1.KO	46.31579	KO	T1	R	T1.R	R.T1.KO
2	R.T1.OE	23.17073	OE	T1	R	T1.R	R.T1.OE
3	R.T1.WT	50.00000	WT	T1	R	T1.R	R.T1.WT
4	R.T2.KO	43.92523	KO	T2	R	T2.R	R.T2.KO
5	R.T2.OE	25.92593	OE	T2	R	T2.R	R.T2.OE
6	R.T2.WT	50.74627	WT	T2	R	T2.R	R.T2.WT

6 rows

Hide



#时空组-稳定性4 网络易损性#------
library("pulsar")
res = natural.con.microp.ts(
  ps.st,
  N = 200,
  r.threshold= 0.8,
  p.threshold=0.05,
  method = "spearman",
  order = "time",
  g1 = "Group",# 分组1
  g2 = "space",# 分组2
  g3 = "time",# 分组3
  ord.g1 =NULL,# 排序顺序
  ord.g2 = NULL, # 排序顺序
  ord.g3 = NULL,# 排序顺序
  norm = F,
  end = N -2,
  start = 0,
  con.method = "pulsar"

)

res[[1]]

Hide


#时空组-稳定性5 组成稳定性#-------
res = community.stability.ts(
  ps.st = ps.st,
  N = 200,
  r.threshold= 0.8,
  p.threshold=0.05,
  method = "spearman",
  order = "space",
  g1 = "Group",# 分组1
  g2 = "space",# 分组2
  g3 = "time",# 分组3
  map.art = NULL, # 人工输入的分组 默认为NULL
  time = F,# 稳定性是否有时间序列
  ord.map = TRUE# map文件是否是已经按照pair要求进行了排序
)

res[[1]]

Hide

# res[[2]]

时空组:多网络单独绘制却填充相同颜色

Hide

#---时空组双网络手动填充相同颜色#-------
library(tidyverse)
library(ggClusterNet)
library(phyloseq)
library(igraph)

ps.st = readRDS("./ps_TS.rds")
ps.st
#> phyloseq-class experiment-level object
#> otu_table()   OTU Table:         [ 2432 taxa and 108 samples ]
#> sample_data() Sample Data:       [ 108 samples by 4 sample variables ]
#> tax_table()   Taxonomy Table:    [ 2432 taxa by 7 taxonomic ranks ]
#> phy_tree()    Phylogenetic Tree: [ 2432 tips and 2431 internal nodes ]
#> refseq()      DNAStringSet:      [ 2432 reference sequences ]



#-----第一组网络绘制

res = Facet.network(
  ps.st= ps.st,# phyloseq对象
  g1 = "Group",# 分组1
  g2 = "space",# 分组2
  g3 = "time",# 分组3
  ord.g1 = c("WT","KO","OE"),# 排序顺序
  ord.g2 = c("B","R") ,# 排序顺序
  ord.g3 = c("T1","T2","T3") ,# 排序顺序
  order = "time", # 出图每行代表的变量
  fill = "Phylum",
  size = "igraph.degree",
  layout_net = "model_maptree2",
  r.threshold=0.8,
  p.threshold=0.01,
  method = "spearman",
  select_layout = TRUE,
  clu_method = "cluster_fast_greedy",
  maxnode = 5
)
#> [1] "B" "R"
#> [1] "T1" "T2" "T3"
#> [1] "WT" "KO" "OE"

p = res[[1]]
p

Hide





#--指定颜色映射，为了多个图使用同一个颜色#-------



#--设定的参数一致
fill = "Phylum"
size = "igraph.degree"
maxnode = 5
row.num = 6

#-从结果中提取网络节点和边数据
node  = res[[2]][[2]]
head(node)

	ID <chr>	X1 <dbl>	X2 <dbl>	elements <chr>	igraph.degree <dbl>	igraph.closeness <dbl>
1	ASV_101	1.8922077	-16.17873	ASV_101	3	1
2	ASV_1050	10.4488008	11.60225	ASV_1050	1	1
3	ASV_113	-0.3111552	18.45045	ASV_113	0	0
4	ASV_1148	11.7082072	-13.80247	ASV_1148	0	0
5	ASV_1149	2.1328120	18.18372	ASV_1149	1	1
6	ASV_115	7.1706596	-16.75757	ASV_115	0	0

6 rows | 1-7 of 23 columns

Hide

edge = res[[2]][[3]]
head(edge)

	X2 <dbl>	Y2 <dbl>	OTU_2 <chr>	OTU_1 <chr>	weight <dbl>	X1 <dbl>	Y1 <dbl>	cor <chr>	group <chr>
1	6.251566	6.418839	ASV_30	ASV_21	-1	5.453796	4.025147	-	R.T3.OE
2	6.251566	6.418839	ASV_30	ASV_29	1	4.842470	7.579802	+	R.T3.OE
3	6.251566	6.418839	ASV_30	ASV_55	1	3.675405	8.282087	+	R.T3.OE
4	6.251566	6.418839	ASV_30	ASV_290	1	3.529926	5.996121	+	R.T3.OE
5	4.842470	7.579802	ASV_29	ASV_21	-1	5.453796	4.025147	-	R.T3.OE
6	4.842470	7.579802	ASV_29	ASV_30	1	6.251566	6.418839	+	R.T3.OE

6 rows | 1-10 of 16 columns

Hide


cb_palette <- c("#ed1299", "#09f9f5", "#246b93", "#cc8e12", "#d561dd", "#c93f00", "#ddd53e",
                "#4aef7b", "#e86502", "#9ed84e", "#39ba30", "#6ad157", "#8249aa", "#99db27", "#e07233", "#ff523f",
                "#ce2523", "#f7aa5d", "#cebb10", "#03827f", "#931635", "#373bbf", "#a1ce4c", "#ef3bb6", "#d66551",
                "#1a918f", "#ff66fc", "#2927c4", "#7149af" ,"#57e559" ,"#8e3af4" ,"#f9a270" ,"#22547f", "#db5e92",
                "#edd05e", "#6f25e8", "#0dbc21", "#280f7a", "#6373ed", "#5b910f" ,"#7b34c1" ,"#0cf29a" ,"#d80fc1",
                "#dd27ce", "#07a301", "#167275", "#391c82", "#2baeb5","#925bea", "#63ff4f")
tem1 = node[,fill] %>% as.matrix() %>% as.vector() %>% unique()
tabf = data.frame(id = tem1,color =cb_palette[1:length(tem1)] )
head(tabf)

	id <chr>	color <chr>
1	Proteobacteria	#ed1299
2	Actinobacteria	#09f9f5
3	Firmicutes	#246b93
4	Bacteroidetes	#cc8e12
5	Unassigned	#d561dd
6	Verrucomicrobia	#c93f00

6 rows

Hide

colnames(tabf)[1] = fill
node1 = node %>% left_join(tabf,by = fill)

p2 <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,
                                  color = cor),
                              data = edge, size = 0.03,alpha = 0.5) +
  geom_point(aes(X1, X2,
                 size = !!sym(size),
                 fill = color ),
             pch = 21, data = node1,color = "gray40") +
  facet_wrap(.~ label,scales="free_y",ncol = row.num ) +
  # geom_text_repel(aes(X1, X2,label = elements),pch = 21, data = nodeG) +
  # geom_text(aes(X1, X2,label = elements),pch = 21, data = nodeG) +
  scale_colour_manual(values = c("#6D98B5","#D48852")) +
  scale_fill_manual(values  = tabf$color,labels = tabf[,fill]) +
  scale_size(range = c(0.8, maxnode)) +
  scale_x_continuous(breaks = NULL) +
  scale_y_continuous(breaks = NULL) +
  theme(panel.background = element_blank(),
        plot.title = element_text(hjust = 0.5)
  ) +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank()
  ) +
  theme(legend.background = element_rect(colour = NA)) +
  theme(panel.background = element_rect(fill = "white",  colour = NA)) +
  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
p2

Hide



# ggsave("cs1.pdf",p2,width = 5*5,height = 4*3)



#-------另外一组网络绘制

Hide


data(ps)

res2 = Facet.network (
  ps.st= ps,# phyloseq对象
  g1 = "Group",# 分组1
  # g2 = "space",# 分组2
  # g3 = "time",# 分组3
  ord.g1 = c("WT","KO","OE"),# 排序顺序
  # ord.g2 = c("B","R") ,# 排序顺序
  # ord.g3 = c("T1","T2","T3") ,# 排序顺序
  order = "time", # 出图每行代表的变量
  fill = "Phylum",
  size = "igraph.degree",
  layout_net = "model_maptree2",
  r.threshold=0.8,
  p.threshold=0.01,

  method = "spearman",
  select_layout = TRUE,
  clu_method = "cluster_fast_greedy",
  maxnode = 5
)
#> [1] "WT" "KO" "OE"

p3 = res2[[1]]
p3

Hide


#-提取节点和边数据
node  = res2[[2]][[2]]
head(node)

	ID <chr>	X1 <dbl>	X2 <dbl>	elements <chr>	igraph.degree <dbl>	igraph.closeness <dbl>
1	ASV_1	6.394069	-0.02808431	ASV_1	1	1
2	ASV_10	3.652617	-21.02187428	ASV_10	0	0
3	ASV_100	7.078794	10.35742308	ASV_100	1	1
4	ASV_101	-3.792132	16.03963922	ASV_101	1	1
5	ASV_102	8.292658	-9.25473746	ASV_102	1	1
6	ASV_103	4.025401	21.04072916	ASV_103	0	0

6 rows | 1-7 of 23 columns

Hide

edge = res2[[2]][[3]]
head(edge)

	X2 <dbl>	Y2 <dbl>	OTU_2 <chr>	OTU_1 <chr>	weight <dbl>	X1 <dbl>	Y1 <dbl>	cor <chr>	group <chr>
1	1.9872552	1.5012669	ASV_22	ASV_32	1	-3.0932354	-2.3169033	+	..OE
2	1.9872552	1.5012669	ASV_22	ASV_77	1	-0.4037329	-0.8876355	+	..OE
3	1.9872552	1.5012669	ASV_22	ASV_95	-1	-0.9859822	2.3370501	-	..OE
4	1.9872552	1.5012669	ASV_22	ASV_214	1	2.3530761	-1.0455685	+	..OE
5	1.9872552	1.5012669	ASV_22	ASV_156	1	-3.4791850	0.5142463	+	..OE
6	-0.4037329	-0.8876355	ASV_77	ASV_32	1	-3.0932354	-2.3169033	+	..OE

6 rows | 1-10 of 16 columns

Hide

node$group %>% unique()
#> [1] "..WT" "..KO" "..OE"

Hide

#----第二组网络用第一组颜色填充#---------


tem1 = node[,fill] %>% as.matrix() %>% as.vector() %>% unique()
tabf2 = data.frame(id = tem1 )
colnames(tabf2)[1] = fill
head(tabf2)

	Phylum <chr>
1	Actinobacteria
2	Proteobacteria
3	Bacteroidetes
4	Unassigned
5	Chloroflexi
6	Firmicutes

6 rows

Hide

# tabf2[1,1] = "wentao"
tabf3 = tabf2 %>% left_join(tabf,by = fill)

num = tabf3$color[is.na(tabf3$color)] %>% length()
if (num == 0) {
  tabf.f = tabf3
}else{
  tem = tabf3 %>% filter(is.na(color))
  tem$color = setdiff(cb_palette,tabf$color)[1:length(tem$color)]
  tabf.f = rbind(tabf3,tem)
}

colnames(tabf.f)[1] = fill
node1 = node %>% left_join(tabf.f,by = fill)

head(node1)

	ID <chr>	X1 <dbl>	X2 <dbl>	elements <chr>	igraph.degree <dbl>	igraph.closeness <dbl>
1	ASV_1	6.394069	-0.02808431	ASV_1	1	1
2	ASV_10	3.652617	-21.02187428	ASV_10	0	0
3	ASV_100	7.078794	10.35742308	ASV_100	1	1
4	ASV_101	-3.792132	16.03963922	ASV_101	1	1
5	ASV_102	8.292658	-9.25473746	ASV_102	1	1
6	ASV_103	4.025401	21.04072916	ASV_103	0	0

6 rows | 1-7 of 24 columns

Hide



p4 <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,
                                  color = cor),
                              data = edge, size = 0.03,alpha = 0.5) +
  geom_point(aes(X1, X2,
                 fill =color,
                 size = !!sym(size)),
             pch = 21, data = node1,color = "gray40") +
  facet_wrap(.~ label,scales="free_y",ncol = row.num ) +
  # geom_text_repel(aes(X1, X2,label = elements),pch = 21, data = nodeG) +
  # geom_text(aes(X1, X2,label = elements),pch = 21, data = nodeG) +
  scale_colour_manual(values = c("#6D98B5","#D48852")) +
  scale_fill_manual(values  = tabf.f$color,labels = tabf.f[,fill]) +
  scale_size(range = c(0.8, maxnode)) +
  scale_x_continuous(breaks = NULL) +
  scale_y_continuous(breaks = NULL) +
  theme(panel.background = element_blank(),
        plot.title = element_text(hjust = 0.5)
  ) +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank()
  ) +
  theme(legend.background = element_rect(colour = NA)) +
  theme(panel.background = element_rect(fill = "white",  colour = NA)) +
  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())
# guides(fill=FALSE)
p4

Hide

# ggsave("cs2.pdf",p4,width = 5*3,height = 4*1)

Hide

library(patchwork)

layout <- "
AAAABB
AAAA##
AAAA##
"
p2+p4 +
  plot_layout(design = layout)

时空组跨域网络：不同分类等级跨域网络

Hide


ps.st = readRDS("./ps_TS.rds")
ps.st
#> phyloseq-class experiment-level object
#> otu_table()   OTU Table:         [ 2432 taxa and 108 samples ]
#> sample_data() Sample Data:       [ 108 samples by 4 sample variables ]
#> tax_table()   Taxonomy Table:    [ 2432 taxa by 7 taxonomic ranks ]
#> phy_tree()    Phylogenetic Tree: [ 2432 tips and 2431 internal nodes ]
#> refseq()      DNAStringSet:      [ 2432 reference sequences ]
#---不同分类等级网络#-------
res = rank.network(
  ps.st= ps.st,# phyloseq对象
  g1 = "Group",# 分组1
  g2 = "space",# 分组2
  g3 = "time",# 分组3
  ord.g1 = NULL, # 排序顺序
  ord.g2 = NULL, # 排序顺序
  ord.g3 = NULL, # 排序顺序
  order = "space", # 出图每行代表的变量
  jj = "Phylum",
  fill = "Phylum",
  method = "spearman",
  clu_method = "cluster_fast_greedy",
  select_layout = TRUE,
  r.threshold=0.8,
  p.threshold=0.01,
  N= 500)

p = res[[1]]
p

Hide

# ggsave("cs1.pdf",p,width = 5*9,height = 4*2)

#--提取相关数据
# res$network.data$cortab$R.T1.WT

时空组跨域网络：细菌、真菌，原生动物，环境因子跨域网络

Hide

#--时空组：跨域网络一体化#---------
#--细菌和真菌ps对象中的map文件要一样
ps.merge <- ggClusterNet::merge16S_ITS(ps16s = ps16s,
                                       psITS = psITS,
                                       N16s = 500,
                                       NITS = 500
)

ps.merge
#> phyloseq-class experiment-level object
#> otu_table()   OTU Table:         [ 1000 taxa and 18 samples ]
#> sample_data() Sample Data:       [ 18 samples by 2 sample variables ]
#> tax_table()   Taxonomy Table:    [ 1000 taxa by 8 taxonomic ranks ]
map =  phyloseq::sample_data(ps.merge)
# ps.merge %>% vegan_tax()


res = corBionetwork.st(
  ps.st= ps.merge,# phyloseq对象
  g1 = "Group",# 分组1
  g2 = NULL,# 分组2
  g3 = NULL,# 分组3
  ord.g1 = NULL, # 排序顺序
  ord.g2 = NULL, # 排序顺序
  ord.g3 = NULL, # 排序顺序
  order = NULL, # 出图每行代表的变量
  fill = "filed",
  size = "igraph.degree",method = "spearman",
  clu_method = "cluster_fast_greedy",
  select_layout = TRUE,layout_net = "model_maptree2",
  r.threshold=0.8,
  p.threshold=0.01,
  maxnode = 5,
  N= 500,scale = TRUE,env = NULL,
  bio = TRUE,minsize = 4,maxsize = 14)

res[[1]]

Hide


# res[[2]]

根际互作生物学研究室 简介

转藏分享

QQ空间 QQ好友新浪微博微信

献花（0） +1

来自：微生信生物 > 《待分类》

举报/认领

0条评论

发表

请遵守用户评论公约

类似文章

微生信生物

关注对话

TA的最新馆藏

高强度人类活动背景下微生物对环境变化的抗性和恢复力
ggCLusterNet更新-模块相似性与模块互联
土壤中微生物对养分获取的相互作用：矿工、清道夫和载体
微生物的抗旱性可能会破坏土壤碳的稳定性
农田土壤细菌群落：分类和功能
微生物组多组学网络分析：统计考虑因素、局限性和机遇（下）

喜欢该文的人也喜欢更多

热门阅读换一换

【原】ggClusterNet包文档更新-增加30%以上网络分析相关内容（2023年2月）

R包安装

导入R包

输入方式

输入方式一

输入方式二

使用ggCLusterNet进行网络分析的过程

corMicro函数用于计算网络相关

制作分组,我们模拟五个分组

PolygonClusterG 根据分组，计算布局位置坐标

nodeadd 节点注释的：用otu表格和分组文件进行注释

计算边

出图

模拟不同的分组–可视化

模拟不同分组效果展示：3个分组

模拟不同的分组查看效果：8个分组

微生物分类分组可视化

微生物分类可视化布局优化1-圆环半径调整PolygonRrClusterG

微生物分类可视化布局优化2-model_filled_circle

微生物分类可视化布局优化3 model_maptree_group

按照网络模块化分组

模块布局算法 model_maptree_group

模块布局算法 model_maptree2

上千OTU相关性计算测试model_igraph2布局

网络属性和节点属性

网络性质计算

节点性质计算

Zipi基于模块对OTU进行分类

扩展-关键OTU挑选

对应随机网络构建和网络参数比对

微生物网络流程

微生物组小网络：model_Gephi.2

微生物大网络：model_maptree2

微生物大网络：model_igraph

微生物大网络：model_igraph2-network.i

跨域网络流程

细菌，真菌，环境因子三者跨域网络

细菌-环境因子跨域网络

细菌和真菌跨域网络

细菌真菌任意分类水平跨域网络

网络稳定性（抗干扰性）

模块相似度

网络鲁棒性（随机去除节点）

网络鲁棒性（去除关键节点）

负相关比例

组成稳定性

网络抗毁性

网络模块化

展示网络模块

基于模块的微生物组成分析

模块计算alpha多样性

基于单个样本计算网络属性

计算模块的Zscore值

模块ZS值，网络属性，环境因子相关联合到一起

时空组:网络相关分析

时空组：多组网络展示-网络性质稳定性等

时空组:多网络单独绘制却填充相同颜色

时空组跨域网络：不同分类等级跨域网络

时空组跨域网络：细菌、真菌，原生动物，环境因子跨域网络

根际互作生物学研究室 简介

根际互作生物学研究室简介