写在前面
邓晔老师GCB出现了不同分类等级的网络：
本次我们使用ggCLusterNet配合写一点代码实现这个过程；
计算相关
注意这里的相关并不是在不同分类等级合并后再做相关，而是在OTU水平做相关，然后合并全部相关结果。
result = cor_Big_micro(ps = pst,
N = 500,
r.threshold= 0.8,
p.threshold=0.05,
method = "spearman")
cor = result[[1]]
dim(cor)
otu_table = as.data.frame(t(vegan_otu(pst)))
tax_table = as.data.frame(vegan_tax(pst))
netClu = data.frame(ID = row.names(tax_table),group = tax_table[[j]] )
netClu$group = as.factor(netClu$group)
result2 = model_maptree_group(
cor = cor,
nodeGroup =netClu,
seed = 12)
node = result2[[1]]# branch=result2[[2]]
# head(branch)
# # ---node节点注释
# nodes = nodeadd(plotcord =node,otu_table = otu_table,tax_table = tax_table)
# head(nodes)
#
# #-----计算边
# edge = edgeBuild(cor = cor,node = node)
# head(edge)
# head(nodes)
# head(branch)
# library("ggalt")
# # BiocManager::install("ggalt")
# p1 <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2,color = cor),
# data = edge, size = 0.1,alpha = 0.01) +
# geom_point(aes(X1, X2,fill = Class,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")) +
# geom_text(aes(x = x, y = y,label = elements), data = branch,pch = 21,size = 5) +
# # geom_encircle(aes(X1, X2,group = Phylum),
# # linetype = 1,alpha = 1, data = nodes,fill = NA,color = "black") +
# scale_size(range = c(8,24))+
# 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#-----计算边
edge = edgeBuild(cor = cor,node = node)
netClu$group %>% unique()
head(edge)
#-去除了分类等级内部相关
#-内部OTU之间的相关都去除
for (i in 1: length(levels(netClu$group))) {
tem1 = levels(netClu$group)[i]
tem2 = netClu$ID[netClu$group == tem1]
if (i == 1) {
tem3 = edge %>%
filter(!(OTU_2 %in% tem2 & OTU_1 %in% tem2))
dim(tem3)
} else if(i > 1) {
tem3 = tem3 %>%
filter(!(OTU_2 %in% tem2 & OTU_1 %in% tem2))
dim(tem3)
}
}
head(tem3)
head(netClu)
tem4 = tem3 %>% left_join(netClu,by = c("OTU_2" = "ID")) %>%
rename(OTU_2g = group) %>%
left_join(netClu,by = c("OTU_1" = "ID")) %>%
rename(OTU_1g = group)
head(tem4)
tem4$tem = paste(tem4$OTU_1g,tem4$OTU_2g,sep = "@")
tem4$weight = abs(tem4$weight)
tem5 = tem4 %>% group_by(tem) %>%
summarise(sum(weight)) %>%
rename( weight = `sum(weight)`) %>%
as.data.frame()
tem5$from = sapply(strsplit(tem5$tem, "[@]"), `[`, 1)
tem5$to = sapply(strsplit(tem5$tem, "[@]"), `[`, 2)
tem5 = tem5 %>% select(from,to,weight)
# 这里将weight的进行标准化，这个也已经不是传统意义
tem5$weight = tem5$weight/sum(tem5$weight)
head(tem5)#--列表边矩阵
mat = tidyfst::df_mat(tem5,from,to,weight)
mat[is.na(mat)] = 0
#-重新做好了 cor矩阵
cor = mat
ps_net = pst %>% tax_glom_wt(j)
ps_net = ps_net %>%
subset_taxa(row.names(tax_table(ps_net))%in%colnames(cor))
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 = "A")
netClu$group = as.factor(netClu$group)
set.seed(12)
library(sna)
# --随机模块化布局-模块多的话，这个布局浪费是非常多的，因为要将每一个模块分开，迭代很多次
result2 = PolygonClusterG(cor = cor,nodeGroup = netClu )
node = result2[[1]]
head(node)
# ---node节点注释
nodes = nodeadd(plotcord =node,otu_table = otu_table,tax_table = tax_table)
head(nodes)
#-----计算边
edge = edgeBuild(cor = cor,node = node)
head(edge)
head(nodes)
# node2$size[is.na(node2$size)] = 1
# colnames(edge)[8] = "cor"
library(ggnewscale)
head(edge)
library(ggrepel)p1 <- ggplot() + geom_segment(aes(x = X1, y = Y1, xend = X2, yend = Y2, color = weight, size = weight ), data = edge,alpha = 0.3) + scale_color_gradientn(colours =c("white","grey80"))+ ggnewscale::new_scale_color() + geom_point(aes(X1, X2,fill = Phylum,size = mean),pch = 21, data = nodes) + ggrepel::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(2, 5)) + 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
ggsave("./cs.pdf",p1,width = 8,height = 6)p1 <- ggplot() +geom_segment(aes(x = X1,y = Y1,xend = X2,yend = Y2,color = weight,size = weight),data = edge,alpha = 1) +scale_color_gradientn(colours =c("white","grey80"))+ggnewscale::new_scale_color() +geom_point(aes(X1, X2,fill = Phylum),pch = 21, data = nodes,size = 4) +ggrepel::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(2, 5)) +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根际互作生物学研究室 简介
联系客服
