【原】ggClusterNet_網絡挖掘與比較-布局算法精簡教程

微生信生物 2022-05-24 發布于北京

展開全文

R包安裝

# 1-包的安裝-------
remotes::install_github("taowenmicro/ggClusterNet",force = TRUE)
library(ggClusterNet)

R包導入

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

查看幫助文檔

# 查看幫助文檔
?network.2()

多組網絡流程

微生物網絡-基于多組pipeline流程

#-2-微生物網絡分析流程#-------

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

result = network.2(ps = ps16s,
                   N = 800,
                   big = TRUE,
                   maxnode = 4,
                   select_layout = TRUE,
                   layout_net = "model_igraph",
                   r.threshold=0.6,
                   p.threshold=0.01,
                   label = FALSE,
                   path = path,
                   zipi = FALSE)

# 多組網絡繪制到一個面板
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)

跨域網絡（細菌-真菌）-pipeline

#2-跨域網絡-16s_ITS#--------

#--創建文件夾用于存儲文件
Envnetplot<- paste("./16S_ITS_network",sep = "")
dir.create(Envnetplot)

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

envRDA = env
row.names(envRDA) = env$ID
envRDA$ID = NULL
head(envRDA)

#--細菌和真菌ps對象中的分組信息要一致
ps.merge <- ggClusterNet::merge16S_ITS(ps16s = ps16s,
                                       psITS = psITS,
                                       N16s = 500,
                                       NITS = 500
)

ps.merge

map =  phyloseq::sample_data(ps.merge)
head(map)

#依據實際情況決定是否修改分組
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"))

result <- corBionetwork(ps = ps.merge,
                        N = 0,
                        lab = data,# 提供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_maptree",
                        clu_method = "cluster_fast_greedy"


)

p = result[[1]]
p

# 全部樣本網絡參數比對
data = result[[2]]

plotname1 = paste(Envnetplot,"/network_all.pdf",sep = "")
ggsave(plotname1, p,width = 10,height = 8)
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)

跨域網絡（細菌-真菌-環境因子）-pipeline

#2-跨域網絡16s-ITS-環境因子的網絡#--------

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

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

library(phyloseq)
#--細菌和真菌ps對象中的map文件要一樣
ps.merge <- ggClusterNet::merge16S_ITS(ps16s = ps16s,
                                       psITS= psITS,
                                       NITS = 200,
                                       N16s = 200)

map =  phyloseq::sample_data(ps.merge)
head(map)
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(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
)

p = result[[1]]
p
# 全部樣本網絡參數比對
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)

ggClusterNet中的網絡布局算法使用

布局算法使用1-PolygonClusterG

#3-布局算法使用1-PolygonClusterG#-------

#-計算相關性
result = corMicro (ps = ps,
                   N = 200,
                   r.threshold=0.8,
                   p.threshold=0.05
)
cor = result[[1]]
dim(cor)
ps_net = result[[3]]
otu_table = as.data.frame(t(vegan_otu(ps_net)))
tax_table = as.data.frame(vegan_tax(ps_net))

#-模擬一個分組
netClu = data.frame(ID = row.names(tax_table),group =rep(1:5,length(row.names(tax_table)))[1:length(row.names(tax_table))] )
netClu$group = as.factor(netClu$group)
set.seed(12)
result2 = PolygonClusterG(cor = cor,nodeGroup =netClu,zoom = 0.8,zoom2 = 0.8 )
node = result2[[1]]
head(node)
# ---node節點注釋
nodes = nodeadd(plotcord =node,otu_table = otu_table,tax_table = tax_table)
#-----計算邊
edge = edgeBuild(cor = cor,node = node)
head(edge)
# colnames(edge)[8] = "cor"
p1 <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = cor),
                              data = edge, size = 0.1,alpha = 0.6) +
  geom_point(aes(X1, X2,fill = Phylum,size = mean),pch = 21, data = nodes) +
  # geom_text_repel(aes(X1, X2,label = elements),pch = 21, data = nodes) +
  # geom_text(aes(X1, X2,label = elements),pch = 21, data = nodes) +
  scale_colour_manual(values = c("#377EB8","#E41A1C")) +
  scale_size(range = c(4, 14)) +
  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())
p1

布局算法使用2-model_igraph

#3-布局算法使用2-model_igraph#-----
p.threshold = 0.05
r.threshold = 0.6
x = ps %>% 
  # scale_micro(method = "TMM") %>%
  vegan_otu() %>% 
  t() %>%
  as.data.frame()
occor<-WGCNA::corAndPvalue(t(x)/colSums(x),method = 'pearson')
mtadj<-multtest::mt.rawp2adjp(unlist(occor$p),proc='BH')
adpcor<-mtadj$adjp[order(mtadj$index),2]
occor.p<-matrix(adpcor,dim(t(x)/colSums(x))[2])
## R value
occor.r<-occor$cor
diag(occor.r) <- 0

# occor.r[occor.p > 0.05|abs(occor.r)<0.4] = 0

occor.r[occor.p > p.threshold | abs(occor.r)< r.threshold] = 0
occor.r[is.na(occor.r)]=0
cor = occor.r

result2 <- model_igraph(cor = cor,
                        method = "cluster_fast_greedy",
                        seed = 12
)
node = result2[[1]]
head(node)

dat = result2[[2]]
head(dat)
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)
#-----計算邊
edge = edgeBuild(cor = cor,node = node)
colnames(edge)[8] = "cor"
head(edge)

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)

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)

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

edge2 = tem4 %>% mutate(color = ifelse(model1 == model2,as.character(model1),"across"),
                        manual = ifelse(model1 == model2,as.character(color1),"#C1C1C1")
)
head(edge2)
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) 

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())

# ggsave("./cs2.pdf",p2,width = 16,height = 14)

根際互作生物學研究室簡介

根際互作生物學研究室是沈其榮院士土壤微生物與有機肥團隊下的一個關注于根際互作的研究小組。本小組由袁軍副教授帶領，主要關注：1.植物和微生物互作在抗病過程中的作用；2 環境微生物大數據整合研究；3 環境代謝組及其與微生物過程研究體系開發和應用。團隊在過去三年中在 isme J， Microbiome， PCE，SBB，Horticulture Research等期刊上發表了多篇文章。歡迎關注微生信生物公眾號對本研究小組進行了解。