woodcorpse的个人博客分享 http://blog.sciencenet.cn/u/woodcorpse

博文

周集中团队Nature子刊中网络图布局的R语言可视化复现

已有 3388 次阅读 2021-10-9 23:45 |个人分类:R|系统分类:科研笔记

Network layout in NCC style

2021年初,周集中团队在Nature Climate Change杂志上发表了一篇气候变暖对微生物网络复杂性与稳定性影响的文章见下文

Nature子刊:周集中团队揭示气候变暖增加了微生物网络的复杂性和稳定性!

图1a 通过多个连续年际采样的网络图展示了气候变化对微生物网络的大致影响,网络图的风格有一种简介之美,同时网络之间差异的展示也非常醒目,十分拿人,该文章后续分析的大部分代码都在GitHub进行了公开分享,但其中没有包含上述网络图

图1a 土壤微生物网络随时间的演替

image

从 2009 (Y0) 到 2014 (Y5) 的六年 (Y) 中构建的分子生态网络的可视化。大模块(≥5 个节点)以不同颜色显示,较小的模块以灰色显示。

这篇推送根据上图风格对网络图进行了复刻,内容包括(但不限于):

  • 节点按照模块划分上色
  • 边按照模块划分上色
  • 多个网络图的批量输出

作图的测试数据和R代码如下:

百度云链接:https://pan.baidu.com/s/1zLfLGlEcTUBFeCGGwBJztw
提取码:nty1

因公众号文章不可修改,如以上链接失效,或想获取代码的更新版,请在“宏基因组”公众号后台回复本文关键字“NCCnetwork”获取最新下载地址

单个图示例

##导入数据以及R包
library(igraph)
library(dplyr)
library(Hmisc)
## 读入OTU/ASV表格,列为样本,行为物种
otu_rare <- read.delim('otutab.txt',header = T,row.names = 1,stringsAsFactors = F)
## 定义一些颜色
col_g <- "#C1C1C1"
cols <- c("#DEB99B" ,"#5ECC6D", "#5DAFD9", "#7ED1E4", "#EA9527", "#F16E1D" ,"#6E4821", "#A4B423",
          "#C094DF" ,"#DC95D8" ,"#326530", "#50C0C9", "#67C021" ,"#DC69AF", "#8C384F", "#30455C", "#F96C72","#5ED2BF")

trt_id <- c('KO','OE','WT') ## 定义样本的关键词,然后从样本名抓取处理的样本
split_otu <- lapply(apply(sapply(trt_id,function(x){grep(x,colnames(otu_rare))}),2,FUN = function(x){otu_rare[,x]}),function(x){x[-(which(rowSums(x)==0)),]})

## 有此处主要聚焦展示绘图方法,构建网络时没有对输入数据进行筛选
## 此外,笔者建议在样本量较小的情况下,推荐采用MENA的方法构建网络
g <- lapply(split_otu,function(x){
    occor<-WGCNA::corAndPvalue(t(x)/colSums(x),method = 'pearson')
    # occor<-WGCNA::corAndPvalue(t(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[is.na(occor.r)]=0
    g <-  graph.adjacency(occor.r, weighted = TRUE, mode = 'undirected')
    # 删除自相关
    g <- simplify(g)
    # 删除孤立节点
    g <- delete.vertices(g, which(degree(g)==0) )
    return(g)
})

 save(g,file = 'network.rda')
 load('network.rda')

计算网络模块

## 提取第一个网络演示
g1 <- g[[1]]
# plot(g[[1]])

## 设置网络的weight,为计算模块性做准备
E(g1)$correlation <- E(g1)$weight
E(g1)$weight <- abs(E(g1)$weight)

## 计算网络模块
set.seed(007)
V(g1)$modularity <- membership(cluster_fast_greedy(g1))

添加节点以及边的颜色

按照模块设置节点的颜色

选取包含节点数量前18个模块赋予不同的颜色,剩余模块赋予灰色

V(g1)$label <- V(g1)$name
V(g1)$label <- NA
modu_sort <- V(g1)$modularity %>% table() %>% sort(decreasing = T)
top_num <- 18
modu_name <- names(modu_sort[1:18])
modu_cols <- cols[1:length(modu_name)]
names(modu_cols) <- modu_name
V(g1)$color <- V(g1)$modularity
V(g1)$color[!(V(g1)$color %in% modu_name)] <- col_g
V(g1)$color[(V(g1)$color %in% modu_name)] <- modu_cols[match(V(g1)$color[(V(g1)$color %in% modu_name)],modu_name)]
V(g1)$frame.color <- V(g1)$color

设置边的颜色

边的颜色与模块颜色保持一致

由于边连接了两个节点,如果两个节点同属于一个模块,我们赋予其模块的颜色

如果两个两个节点属于不同模块,我们赋予其灰色

E(g1)$color <- col_g
for ( i in modu_name){
  col_edge <- cols[which(modu_name==i)]
  otu_same_modu <-V(g1)$name[which(V(g1)$modularity==i)]
  E(g1)$color[(data.frame(as_edgelist(g1))$X1 %in% otu_same_modu)&(data.frame(as_edgelist(g1))$X2 %in% otu_same_modu)] <- col_edge
}

计算网络的layout并输出

我们基于layout_with_fr的算法计算layout

当我们的节点数量大于1000时,节点会按照坐标轴排布,出现异常的layout,因此建议设置grid参数为'nogrid'

设置font.main=4,使得标题为斜体加粗

# 计算 layout
sub_net_layout <- layout_with_fr(g1, niter=999,grid = 'nogrid')
## 可视化并输出
par(font.main=4)
plot(g1,layout=sub_net_layout, edge.color = E(g1)$color,vertex.size=2)
title(main = paste0('Nodes=',length(V(g1)$name),', ','Edges=',nrow(data.frame(as_edgelist(g1)))))

# pdf(paste0("Example 1.pdf"), encoding="MacRoman", width=6, height=6)
# par(font.main=4)
# plot(g1,layout=sub_net_layout, edge.color = E(g1)$color,vertex.size=2)
# title(main = paste0('Nodes=',length(V(g1)$name),', ','Edges=',nrow(data.frame(as_edgelist(g1)))))
# dev.off()

image

多图批量产出

pdf(paste0("Example 2.pdf"), encoding="MacRoman", width=15, height=9)
par(mfrow=c(1,3),mar=c(0,0,1,0),font.main=4)
for(i in 1:3){
  g1 <- g[[i]]
  E(g1)$correlation <- E(g1)$weight
  E(g1)$weight <- abs(E(g1)$weight)
  set.seed(007)
  V(g1)$modularity <- membership(cluster_fast_greedy(g1))

  V(g1)$label <- V(g1)$name
  V(g1)$label <- NA
  modu_sort <- V(g1)$modularity %>% table() %>% sort(decreasing = T)

  top_num <- 18
  modu_name <- names(modu_sort[1:18])
  modu_cols <- cols[1:length(modu_name)]
  names(modu_cols) <- modu_name
  V(g1)$color <- V(g1)$modularity
  V(g1)$color[!(V(g1)$color %in% modu_name)] <- col_g
  V(g1)$color[(V(g1)$color %in% modu_name)] <- modu_cols[match(V(g1)$color[(V(g1)$color %in% modu_name)],modu_name)]
  V(g1)$frame.color <- V(g1)$color

  E(g1)$color <- col_g
  for ( i in modu_name){
    col_edge <- cols[which(modu_name==i)]
    otu_same_modu <-V(g1)$name[which(V(g1)$modularity==i)]
    E(g1)$color[(data.frame(as_edgelist(g1))$X1 %in% otu_same_modu)&(data.frame(as_edgelist(g1))$X2 %in% otu_same_modu)] <- col_edge
  }

  sub_net_layout <- layout_with_fr(g1, niter=999,grid = 'nogrid')
  plot(g1,layout=sub_net_layout, edge.color = E(g1)$color,vertex.size=2)
  title(main = paste0('Nodes=',length(V(g1)$name),', ','Edges=',nrow(data.frame(as_edgelist(g1)))))
}
dev.off()

image

作者:flyfly

责编:马腾飞 南京农业大学

审核:刘永鑫 中科院遗传发育所

猜你喜欢

写在后面

为鼓励读者交流、快速解决科研困难,我们建立了“宏基因组”专业讨论群,目前己有国内外5000+ 一线科研人员加入。参与讨论,获得专业解答,欢迎分享此文至朋友圈,并扫码加主编好友带你入群,务必备注“姓名-单位-研究方向-职称/年级”。技术问题寻求帮助,首先阅读《如何优雅的提问》学习解决问题思路,仍末解决群内讨论,问题不私聊,帮助同行。
image

学习扩增子、宏基因组科研思路和分析实战,关注“宏基因组”
image

点击阅读原文,跳转最新文章目录阅读
https://mp.weixin.qq.com/s/5jQspEvH5_4Xmart22gjMA



https://m.sciencenet.cn/blog-3334560-1307370.html

上一篇:ISME:南农张瑞福组揭示芽孢杆菌通过代谢互作刺激常驻根际微生物促进植物生长!
下一篇:RDA_环境因子_群落结构_统计检验_可视化

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-3-28 19:35

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部