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)
)
ps
rank.names(ps)
#或者直接从网站读取,但是由于github在国外,所以容易失败
metadata = read.delim("https://raw.githubusercontent.com/taowenmicro/R-_function/main/metadata.tsv",row.names = 1)
otutab = read.delim("https://raw.githubusercontent.com/taowenmicro/R-_function/main/otutab.txt", row.names=1)
taxonomy = read.table("https://raw.githubusercontent.com/taowenmicro/R-_function/main/taxonomy.txt", row.names=1)
# tree = read_tree("https://raw.githubusercontent.com/taowenmicro/R-_function/main/otus.tree")
# rep = readDNAStringSet("https://raw.githubusercontent.com/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)
)
按照丰度过滤微生物表格,并却计算相关矩阵,按照指定的阈值挑选矩阵中展示的数值。调用了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
不同的模块按照分组聚集成不同的圆,并且圆形的大小一样。如果一个分组只有一个点,则这个点坐落在圆心位置。
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函数只是提供了简单的用注释函数,用户可以自己在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)
这是网络布局的基础,无论是什么聚类布局,都需要制作一个分组文件,这个文件有两列,一列是节点,一列是分组信息,这个分组信息名称必须设定为: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)
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)
结果发现这些高丰度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)
用实心点作为每个模块的布局方式
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)
用实心点作为每个模块布局方式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)
按照物种组成分类完成网络分析其实并不常用,更多的是按照模块分组,进行网络可视化。
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:不在可以将每个模块独立区分,而是将模块聚拢,并在整体布局上将离散的点同这些模块一同绘制到同心圆内。控制了图形的整体布局为圆形。
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)
这个布局最近几年文章上使用非常多。
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
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)
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)
使用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)
大网络运算时间会比较长,这里我没有计算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)
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在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)
ID <fct> | ratio <dbl> | |
---|---|---|
1 | KO | 41.15139 |
2 | OE | 34.19023 |
3 | WT | 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)
这里将指定模块和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
按照样本计算模块的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
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]]
根际互作生物学研究室 简介
联系客服