【原】ggCLusterNet更新-模塊相似性與模塊互聯

微生信生物 2024-01-22 發布于北京

展開全文

寫在前面

之前好多人咨詢上圖的做法，這個有點麻煩，因為要在不同的facet之間進行繪圖，所以用到了grid，當時我就自己用了一下，沒有放到ggClusterNet里面，昨天到今天我重新寫了一下，從此開始多網絡比對，facet之間繪圖的新篇章；還是要提醒大家一下，使用下面代碼之前先更新ggClusterNwet包;還有就是這個功能只是初步功能，我應該會在接下來的一個月內進行大量的更新；



install.packages("BiocManager")library(BiocManager)install("remotes")install("tidyverse")install("tidyfst")install("igraph")install("sna")install("phyloseq")install("ggalluvial")install("ggraph")install("WGCNA")install("ggnewscale")install("pulsar")install("patchwork")remotes::install_github("taowenmicro/EasyStat")remotes::install_github("taowenmicro/ggClusterNet")

實戰

#--更新流程，自由度更高的網絡分析流程
library(tidyverse)
library(ggClusterNet)
library(phyloseq)
library(igraph)
library("ggalt")
library(ggnewscale)
library(gtable)
library(grid)

network.pip函數計算網絡

#方法計算相關矩陣#----
tab.r = network.pip(
  ps = ps16s,
  N = 500,
  # ra = 0.05,
  big = TRUE,
  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
)

#  從nerwork.pip中提出來的節點和矩陣#-------
dat = tab.r[[2]]
node = dat$net.cor.matrix$node
edge = dat$net.cor.matrix$edge

cortab = dat$net.cor.matrix$cortab
library(tidyfst)
library(grid)

計算模塊及其相似模塊

res1 = module.compare.m(
  ps = NULL,
  corg = cortab,
  zipi = FALSE,
  zoom = 0.2,
  padj = F,
  n = 3)

#不同分組使用一個圓圈展示，圓圈內一個點代表一個模塊，相連接的模塊代表了相似的模塊。p1 = res1[[1]]
# p1
#--提取模塊的OTU，分組等的對應信息
dat1 = res1[[2]]
head(dat1)
colnames(dat1)[2] = "group.module"
#模塊相似度結果表格
dat2 = res1[[3]]
head(dat2)
dat2$module2 = row.names(dat2) %>% strsplit("[-]") %>%
  sapply(`[`, 2)
dat2$module1 = row.names(dat2) %>% strsplit("[-]") %>%
  sapply(`[`, 1)
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")

相似模塊合并

# 相似模塊合并#-------
id = dat2$module1 %>% unique()
for (i in 1:length(id)) {
  tem1 = dat2 %>% filter(module1 == id[i]) %>% filter(p_adj < 0.05) %>% .$module2
  tem2 = c(id[i],tem1)
  tem3 = data.frame(ID = tem2,class = paste0("sameM.",i))

  if (i == 1) {
    tem4 = tem3
  } else{
    tem4 = rbind(tem4,tem3)
  }
}

tem4$ID %>% 
  unique() %>%
  length()

id = tem4$ID %>% table() %>% as.data.frame() %>% 
  arrange(desc(Freq)) %>%
  filter(Freq > 1) %>%
  .$`.` %>% as.character()

for (i in 1:length(id)) {
  tem5 = tem4 %>% filter(ID == id[i]) %>%.$class
  print( tem5)

  sid = tem5[str_detect(tem5,"MNEW")]
  if (length(sid)> 0) {
    for (ii in 1:length(tem5)) {
      tem4$class = gsub(tem5[ii],sid,tem4$class)
    }
  } else{
    for (ii in 1:length(tem5)) {
      tem4$class = gsub(tem5[ii],paste0("MNEW",i),tem4$class)
    }
  }
  tem4 = tem4 %>% distinct(ID,class,.keep_all = TRUE)
}

tem4$class %>% unique()

節點信息和模塊信息合并

node3 = node %>% left_join(dat1,by = c("ID","Group"))
head(node3)

node4 = node3 %>% left_join(tem4,by = c("group.module" = "ID"))
head(node4)
node4  = add.id.facet(node4,"Group")
head(node4)
tt.1 = node4$group.module %>% table() %>% as.data.frame() %>% arrange(desc(Freq)) %>%.$`.` %>%
  head(12) %>% as.character()

node4$label = factor(node4$label,levels =as.character(unique( node4$label)[1:3]))
node4$Group = factor(node4$Group,levels = node4$Group %>% unique())

# p <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = cor
# ),
# data = edge, size = 0.03,alpha = 0.3) +
#   geom_point(aes(X1, X2,
#                  fill = class,
#                  size = igraph.degree),
#              pch = 21, data = node4,color = "gray40") +
#   ggnewscale::new_scale_fill() + 
#   geom_encircle(aes(X1, X2,group = group.module,fill = group.module),
#                 linetype = 2,alpha = 0.1, data = node4 %>% filter(group.module %in%tt.1),
#                 s_shape = 1, expand=0
#                 ) +
#   facet_wrap(.~ Group,scales="free_y",nrow = 1) +
#   # geom_text_repel(aes(X1, X2,label = elements),pch = 21, data = nodeG) +
#   # geom_text(aes(X1, X2,label = ID),pch = 21, data = node4,size = 2) +
#   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 <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = cor
),
data = edge, size = 0.03,alpha = 0.3) +
  geom_point(aes(X1, X2,
                 fill = class,
                 size = igraph.degree),
             pch = 21, data = node4,color = "gray40") +
  ggnewscale::new_scale_fill() + 
  geom_encircle(aes(X1, X2,group = group.module,fill = group.module),
                linetype = 2,alpha = 0.1, data = node4 %>% filter(group.module %in%tt.1),
                expand = 0.03
  ) +
  facet_wrap(.~ Group,scales="free_y",nrow = 1) +
  # geom_text_repel(aes(X1, X2,label = elements),pch = 21, data = nodeG) +
  # geom_text(aes(X1, X2,label = ID),pch = 21, data = node4,size = 2) +
  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())ggsave("cs0.pdf",p,width = 20,height = 6)

g = p
id = node4$class %>% unique()
id = id[!is.na(id)]
for (j in 1:length(id)) {

  node4 %>% filter(class == id[j]) %>% distinct(class,group,.keep_all = TRUE)
  iid = node4 %>% filter(class == id[j]) %>% distinct(class,group,.keep_all = TRUE) %>% .$id.facet



  iid2 = iid %>% strsplit( "[_]") %>% sapply(`[`, 2)  %>% 
    as.numeric()

  iid3 = iid %>% strsplit( "[_]") %>% sapply(`[`, 1) 
  # iid3 = c("Group3","Group1")
  iid3 = levels(node4$Group) [levels(node4$Group) %in% iid3]
  iid4 = (1:length(levels(node4$Group)))[levels(node4$Group) %in% iid3]
  tab = data.frame(Group =  iid3,facet = iid4)
  tab2 = tab$facet[tab$Group%in%c(node4 %>% filter(class == id[j]) %>% distinct(class,group,.keep_all = TRUE) %>% .$group)]


  tem = data.frame(ID = c(0,iid2[1:(length(iid2)-1)]),ID2 = iid2,facet =c(0, tab2[1:(length(tab2)-1)]),tab2) 
  tem2 = tem [-1,]
  tem2


  for (i in 1:dim(tem2)[1]) {
    g <- line.across.facets.network(g, from=tem2[i,3], to=tem2[i,4], from_point_id=tem2[i,1], to_point_id=tem2[i,2])

  }
}

# grid.draw(g)ggsave("cs12.pdf",g,width = 20,height = 6)
根際互作生物學研究室 簡介