微生物联合挖掘教程连载-PC-10-网络稳定性分析、分组之间距离比较和可视化、网络模块化分析

library(ggClusterNet)library(phyloseq)library(tidyverse)library(igraph)library(tidyfst)
netpath = paste(otupath,"/network_stab/",sep = "")dir.create(netpath)ps = ps0res = module.compare.m(  ps = ps,  Top = 200,  degree = TRUE,  zipi = FALSE,  r.threshold= 0.8,  p.threshold=0.05,  method = "spearman",  padj = F,  n = 3)
p = res[[1]]p
dat = res[[2]]head(dat)dat2 = res[[3]]head(dat2)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)
p2 = ggplot(dat2) + geom_bar(aes(x = cross,fill = cross)) +  labs(x = "",       y = "numbers.of.similar.modules"  )+ theme_classic()
p2getwd()library(ggplot2)#--发现分组1和分组3网络更相似一些FileName <- paste(netpath,"module.compare.groups.pdf", sep = "")ggsave(FileName, p, width = 10, height = 10)
FileName <- paste(netpath,"numbers.of.similar.modules.pdf", sep = "")ggsave(FileName, p2, width = 8, height = 8)
FileName <- paste(netpath,"module.otu.csv", sep = "")write.csv(dat,FileName, quote = F)
FileName <- paste(netpath,"module.compare.groups.csv", sep = "")write.csv(dat2,FileName, quote = F)
```### 鲁棒性```{r}

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

# dir.create("./Robustness_Random_removal/")path = paste(netpath,"/Robustness_Random_removal/",sep = "")fs::dir_create(path)write.csv(dat,          paste(path,"random_removal_network.csv",sep = ""))ggsave(paste(path,"random_removal_network.pdf",sep = ""),  p,width = 8,height = 4)

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]]pdat = res[[2]]
path = paste(netpath,"/Robustness_Targeted_removal/",sep = "")fs::dir_create(path)write.csv(dat,          paste(path,"Robustness_Targeted_removal_network.csv",sep = ""))ggsave(paste(path,"Robustness_Targeted_removal_network.pdf",sep = ""),  p,width = 8,height = 4)

group = "Group"  map = as.data.frame(sample_data(ps))  alppath = paste(otupath,"/distance/",sep = "")  dir.create(alppath)  gro = map[,group] %>% unique()  colnames(gro) = "group"  conbgroup = combn(gro$group,2)  # 计算包括终点均值的所有样品bray距离  bray_curtis = vegan::vegdist(vegan_otu(ps), method = "bray")  bray_curtis = as.matrix(bray_curtis)    for (i in 1:dim(conbgroup)[2]) {    a = conbgroup[,i]    map = as.data.frame(sample_data(ps))    # head(map)        chose1 = map[as.matrix(map[,group]) %>% as.vector() == a[1],] %>% row.names()    chose2 = map[as.matrix(map[,group]) %>% as.vector() == a[2],] %>% row.names()        dat = data.frame(group = paste(a[1],a[2],sep = "_VS_"), Distance =bray_curtis[chose1,chose2] %>% as.dist() %>% as.vector() )    head(dat)        if (i == 1) {      table = dat    }        if (i != 1) {      table = rbind(table,dat)    }  }    head(table)  table$id = 1:dim(table)[1]  data <- table %>% dplyr::select(id,everything())      result = EasyStat::MuiKwWlx(data = data,num = c(3))  FileName <- paste(alppath,"/distance_label.csv", sep = "")  write.csv(result,FileName,sep = "")  FileName <- paste(alppath,"/distance_index.csv", sep = "")  write.csv(data,FileName,sep = "")    result1 = EasyStat::FacetMuiPlotresultBox(data = data,num = c(3),result = result,sig_show ="abc",ncol = 1 )  p1_1 = result1[[1]] +     # scale_x_discrete(limits = axis_order) +     mytheme2 +    guides(fill = guide_legend(title = NULL)) +    scale_fill_manual(values = colset4)  p1_1  res = EasyStat::FacetMuiPlotresultBar(data = data,num = c(3),result = result,sig_show ="abc",ncol = 1)  p1_2 = res[[1]]+     # scale_x_discrete(limits = axis_order) +     guides(color = FALSE) +    mytheme2+     guides(fill = guide_legend(title = NULL))+    scale_fill_manual(values = colset4)  p1_2  res = EasyStat::FacetMuiPlotReBoxBar(data = data,num = c(3),result = result,sig_show ="abc",ncol = 1)  p1_3 = res[[1]]+     # scale_x_discrete(limits = axis_order) +     mytheme2 +     guides(fill = guide_legend(title = NULL))+    scale_fill_manual(values = colset4)  p1_3  FileName <- paste(alppath,"distance_box", ".pdf", sep = "")  ggsave(FileName, p1_1, width = ((4+ gnum) ), height =8,limitsize = FALSE)    FileName <- paste(alppath,"distance_bar", ".pdf", sep = "")  ggsave(FileName, p1_2, width = ((4+gnum) ), height = 8,limitsize = FALSE)    FileName <- paste(alppath,"distance_boxbar", ".pdf", sep = "")  ggsave(FileName, p1_3, width = ((4+gnum) ), height = 8,limitsize = FALSE)    FileName <- paste(alppath,"distance_box", ".jpg", sep = "")  ggsave(FileName, p1_1, width = (( 4+gnum) ), height =8,limitsize = FALSE)    FileName <- paste(alppath,"distance_bar", ".jpg", sep = "")  ggsave(FileName, p1_2, width = ((4+ gnum) ), height = 8,limitsize = FALSE)    FileName <- paste(alppath,"distance_boxbar", ".jpg", sep = "")  ggsave(FileName, p1_3, width = ((4+ gnum) ), height = 8,limitsize = FALSE)