【原】ggClusterNet更新多网络比对教程

写在前面

之前我们使用network.2()函数或者network.i()运行网络分析流程，这个函数其实是基于单个网络运行，然后使用for循环进行网络分析内容的单个保存。其实这达不到网络比对的要求。这里更新的更多的关注不同网络之间的比对。

实战

本次更新函数network.pip

推荐使用sparcc进行网络分析，适用于小样本数据。相关参数和之前函数的都一样，只是输出内容不一样。

rm(list=ls())

#--更新流程，自由度更高的网络分析流程
library(tidyverse)
library(ggClusterNet)
library(phyloseq)
library(igraph)

#--sparcc方法计算相关矩阵#----
tab.r = network.pip(
  ps = ps,
  N = 200,
  # ra = 0.05,
  big = FALSE,
  select_layout = FALSE,
  layout_net = "model_maptree2",
  r.threshold = 0.6,
  p.threshold = 0.05,
  maxnode = 2,
  method = "spearman",
  label = FALSE,
  lab = "elements",
  group = "Group",
  fill = "Phylum",
  size = "igraph.degree",
  zipi = TRUE,
  ram.net = TRUE,
  clu_method = "cluster_fast_greedy",
  step = 100,
  R=10,
  ncpus = 1
  )

这里的输出没有直接输出到文件夹中，全部保存在R环境中，可以进行后续操作，最典型的就是将相关矩阵保存，方便后面各种计算相关矩阵的网络相关分析。

#  建议保存一下输出结果为R对象，方便之后不进行相关矩阵的运算，节约时间
# saveRDS(tab.r,"network.pip.sparcc.rds")
# tab.r = readRDS("./network.pip.sparcc.rds")

#-提取全部图片的存储对象
plot = tab.r[[1]]# 提取网络图可视化结果
p0 = plot[[1]]p0.1 = plot[[2]]
#--随机网络幂率分布
p0.2 = plot[[3]]
#--提取相关矩阵,这是一个list存储的相关矩阵
dat = tab.r[[2]]
cortab = dat$net.cor.matrix$cortab

# 大型相关矩阵跑出来不容易，建议保存，方便各种网络性质的计算
# saveRDS(cortab,"cor.matrix.all.group.rds")
# cor = readRDS("./cor.matrix.all.group.eds")

网络显著性比对

这是使用的是fisher检验

#--网络显著性比较#-----
dat = module.compare.net.pip(
  ps = NULL,
  corg = cortab,
  degree = TRUE,
  zipi = FALSE,
  r.threshold= 0.8,
  p.threshold=0.05,
  method = "spearman",
  padj = F,
  n = 3)
res = dat[[1]]
head(res)

使用network.pip输出结果进行可视化订制

这里提供节点和边文件输出，进行高度自定义网络绘制 例如这里我们看到图例有点长，那何不将网络纵向展示，拉一下，显得图例没那么长 其次正负相关许多人都会觉得不好看，这里修改一下

#--网络自定义可视化#-----

dat = tab.r[[2]]
node = dat$net.cor.matrix$node
edge = dat$net.cor.matrix$edge
head(edge)
head(node)

p <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2
                                  ),
                              data = edge, size = 0.03,alpha = 0.5,
                             color = "yellow") +
  geom_point(aes(X1, X2,
                 fill = Phylum,
                 size = igraph.degree),
             pch = 21, data = node,color = "gray40") +
  facet_wrap(.~ label,scales="free_y",nrow = 3) +
  # 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_hue()+
  scale_size(range = c(0.8, 5)) +
  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())
p

ZiPi可是化订制

#--zipi可视化-定制#-----
dat.z = dat$zipi.data
head(dat.z)
x1<- c(0, 0.62,0,0.62)
x2<- c( 0.62,1,0.62,1)
y1<- c(-Inf,2.5,2.5,-Inf)
y2 <- c(2.5,Inf,Inf,2.5)
lab <- c("peripheral",'Network hubs','Module hubs','Connectors')
roles.colors <- c("#E6E6FA","#DCDCDC","#F5FFFA", "#FAEBD7")
tab = data.frame(x1 = x1,y1 = y1,x2 = x2,y2 = y2,lab = lab)
tem = dat.z$group %>% unique() %>% length()
for ( i in 1:tem) {
  if (i == 1) {
    tab2 = tab
  } else{
    tab2 = rbind(tab2,tab)
  }
}

p <- ggplot() +
  geom_rect(data=tab2,
            mapping=aes(xmin=x1,
                        xmax=x2,
                        ymin=y1,
                        ymax=y2,
                        fill = lab))+
  guides(fill=guide_legend(title="Topological roles")) +
  scale_fill_manual(values = roles.colors)+
  geom_point(data=dat.z,aes(x=p, y=z,color=module)) + theme_bw()+
  guides(color= F) +
  ggrepel::geom_text_repel(data = dat.z,
                           aes(x = p, y = z,
                               color = module,label=label),size=4)+
  # facet_wrap(.~group) +
  facet_grid(.~ group, scale='free') +
  theme(strip.background = element_rect(fill = "white"))+
  xlab("Participation Coefficient")+ylab(" Within-module connectivity z-score")
p

随机网络订制

#--随机网络，幂率分布#-------
dat.r = dat$random.net.data

p3 <- ggplot(dat.r) +
  geom_point(aes(x = ID,y = network,
                 group =group,fill = group),pch = 21,size = 2) +
  geom_smooth(aes(x = ID,y = network,group =group,color = group))+
  facet_grid(.~g,scales = "free") +
  theme_bw() + theme(
    plot.margin=unit(c(0,0,0,0), "cm")
  )
p3

使用network.pip输出的相关矩阵缩短网络稳定性计算时间

# 网络稳定性#--------
##--模块比对#-----
library(tidyfst)

res1 = module.compare.m(
  ps = NULL,
  corg = cortab,
  zipi = FALSE,
  zoom = 0.9,
  padj = F,
  n = 3)
#不同分组使用一个圆圈展示，圆圈内一个点代表一个模块，相连接的模块代表了相似的模块。p1 = res1[[1]]
p1
#--提取模块的OTU，分组等的对应信息
dat1 = res1[[2]]
head(dat1)
#模块相似度结果表格
dat2 = res1[[3]]
dat2$m1 = dat2$module1 %>% strsplit("model") %>%
  sapply(`[`, 1)
dat2$m2 = dat2$module2 %>% strsplit("model") %>%
  sapply(`[`, 1)
dat2$cross = paste(dat2$m1,dat2$m2,sep = "_Vs_")
# head(dat2)
dat2 = dat2 %>% filter(module1 != "none")

p2 = ggplot(dat2) + geom_bar(aes(x = cross,fill = cross)) +
  labs(x = "",
       y = "numbers.of.similar.modules"
  )+ theme_classic()

p2

#---去除关键节点-网络鲁棒性#-----
# 鲁棒性计算需要物种丰富，所以即使计算好了相关矩阵，也需要输入ps对象
library(patchwork)
res2= Robustness.Targeted.removal(ps = ps,
                                 corg = cortab,
                                 degree = TRUE,
                                 zipi = FALSE
)
p3 = res2[[1]]
p3
#提取数据
dat4 = res2[[2]]
#--随即取出任意比例节点-网络鲁棒性#---------
res3 = Robustness.Random.removal(ps = ps,
                                corg = cortab,
                                Top = 0
)
p4 = res3[[1]]
p4
#提取数据
dat5 = res3[[2]]
# head(dat5)
#--计算负相关的比例#----
res4 = negative.correlation.ratio(ps = ps16s,
                                 corg = cortab,
                                 # Top = 500,
                                 degree = TRUE,
                                 zipi = FALSE)

p5 = res4[[1]]
p5
dat6 = res4[[2]]
#-负相关比例数据
# head(dat6)
#----时间群落稳定性-只有pair样本使用#-----
treat = ps %>% sample_data()
treat$pair = paste( "A",c(rep(1:6,3)),sep = "")
# head(treat)
sample_data(ps) = treat
#一般性的处理，没有时间梯度的，这里设定time为F，意味着每两个群落的组合都进行比较
res5 = community.stability( ps = ps,
                           corg = cortab,
                           time = FALSE)
p6 = res5[[1]]
p6
dat7 = res5[[2]]
#--网络抗毁性#------
library("pulsar")
res6 = natural.con.microp (
  ps = ps,
  corg = cortab,
  norm = TRUE,
  end = 150,# 小于网络包含的节点数量
  start = 0
)
p7 = res6[[1]]
p7
dat8  = res6[[2]]

基于输出矩阵运行网络模块化分析

#--这里我们指定三个模块，绘制这三个模块物种组成图表
# select.mod = c("model_1","model_2","model_3")
select.mod ="no"# 可以选择展示那些模块编号从1开始，往后模块大小会逐渐见笑。plots = list() # 保存网络模块化展示结果
plots2 = list()# 保存模块化展示结果2
alldat = list()# 模块数据存储
plots3 = list()# 模块物种组成
plots4 = list()# 模块多样性
alldat2 = list()# 模块物种组成数据
alldat3 = list()# 模块多样性数据
i  = 1
for (i in 1:length(names(cortab))) {
  # 模块网络展示#--------
  resu  = module_display.2(
    pst = ps,
    corg = cortab[[names(cortab)[i]]],
    select.mod = "no",#选择指定模块可视化
    num = 5, # 模块包含OTU数量少于5个的不展示,
    leg.col = 3
  )

  # 全部模块输出展示
  p1 = resu[[1]]+ labs(title = names(cortab)[i]) +
    theme(
          plot.title = element_text(hjust=0.5),#title居中
          plot.margin = unit(c(3,3,4,4),"lines")#绘画区域与边界的距离
    )
  p1
  #--提取模块信息：这里区分好不同模块的OTU，以及每个模块包含的网络边的数量
  dat = resu[[5]]
  head(dat)

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

  #按照模块着色，模块之间的边着色为灰色展示整个网络。p2 = resu[[3]] + labs(title = names(cortab)[i]) +
    theme(
          plot.title = element_text(hjust=0.5)
    )
  p2
  alldat[[names(cortab)[i]]] = dat2
  plots[[names(cortab)[i]]] = p1
  plots2[[names(cortab)[i]]] = p2

  #--模块信息，共有三列，第一列是OTU，第二列是吗模块的分类，第三列是模块捏的边数量
  if (select.mod == "no") {
    mod1 = resu$mod.groups 
  } else{
    mod1 = resu$mod.groups %>%
      filter(group %in% select.mod)
    head(mod1)
  }

  #-计算模块物种组成信息
  pst = ps %>%
    subset_samples.wt("Group",names(cortab)[i]) %>%
    subset_taxa.wt("OTU",colnames(cortab[[names(cortab)[i]]])) %>%
    filter_taxa(function(x) sum(x ) > 0, TRUE) %>%
    scale_micro("rela")

  # #-选择模块的微生物组成#-------
  # res = module_composition(pst = pst,mod1 = mod1,j = "Family")
  # p3 = res[[1]]
  # p3

  #按照属水平进行绘制
  res = module_composition(pst = pst,mod1 = mod1,j = "Genus")
  p3 = res[[1]]
  p3
  plots3[[names(cortab)[i]]] = p3
  #---提取出图数据物种组成数据导出
  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)
  alldat2[[names(cortab)[i]]] = tab.res
  # 没有标准化的出图数据
  tab.res.2 = res[[4]]$bundance
  # head(tab.res.2)
  # 标准化到100%的出图数据
  tab.res.3 = res[[4]]$relaabundance
  # head(tab.res.3)

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

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

  p5 = result[[1]]
  p5
  plots4[[names(cortab)[i]]] = p5
  #--alpha多样性数据
  plotd = result[[4]]$alpha

  alldat3[[names(cortab)[i]]] = plotd

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

}

p1  = ggpubr::ggarrange(plotlist = plots,
                        common.legend = FALSE, legend="right",ncol = 3,nrow = 1)
p1
p2  = ggpubr::ggarrange(plotlist = plots2,
                        common.legend = FALSE, legend="right",ncol = 3,nrow = 1)
# p2
# ggsave("cs.pdf",p2,width = 25,height = 5)

p3  = ggpubr::ggarrange(plotlist = plots3,
                        common.legend = FALSE, legend="right",ncol = 1,nrow =3)
p3
# ggsave("cs.pdf",p3,width = 25,height = 15)

p4  = ggpubr::ggarrange(plotlist = plots4,
                        common.legend = FALSE, legend="right",ncol = 1,nrow = 3)
p4
# ggsave("cs.pdf",p4,width = 35,height = 15)

基于多网络计算基于样本的网络属性

#-基于单个样本进行网络属性计算和合并
for (i in 1:length(names(cortab))) {
  pst = ps %>%
    subset_samples.wt("Group",names(cortab)[i]) %>%
    subset_taxa.wt("OTU",colnames(cortab[[names(cortab)[i]]])) %>%
    filter_taxa(function(x) sum(x ) > 0, TRUE) %>%
    scale_micro("rela")

  dat.f = netproperties.sample(pst = pst,cor = cortab[[names(cortab)[i]]])
  head(dat.f)
  if (i == 1) {
    dat.f2 = dat.f
  }else{
    dat.f2 = rbind( dat.f2, dat.f)
  }

}

计算模块ZScore值用于与其他指标做相关

# 每个模块丰度#------
tab.f = list()
for (i in 1:length(names(cortab))) {
  dat = ZS.net.module(
    pst = ps,
    # Top = 500,
    corg =  cortab[[names(cortab)[i]]],
    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)
  head(data)
  tab.f[[names(cortab)[i]]] = data
}
根际互作生物学研究室 简介