ggClusterNet---一条代码完成全内容微生物网络
The following article is from 微生信生物 Author 文涛聊科研
写在前面
今年的毕业生应该是6月底离校，但是这位朋友到了九月份终于完善了手上的工作，对这几年来了一个完善的结束，感谢送给我的一方印章，也祝铠鸣大展宏图。再有就是ggclusterNet功能也完善了90%，剩下的慢慢补充上来，尤其是代码的优化要抽出些时间多和大佬交流，优化。要知道这部分实在是也很头疼。原计划今年将微生物组分析的四个R包够写好，想来想去似乎目前只有easystat和ggClusterNet可以顺利完成，我也会抓紧时间在今年结束的时候至少完成一个包（phyloBiota）的功能。
导入所需要的R包
#--导入所需R包#-------
library(phyloseq)
library(igraph)
library(network)
library(sna)
library(ggplot2)
library(ggClusterNet)数据导入
我在ggClusterNet包中加入了一个phyloseq封装的微生物组数据，这份数据来自刘永鑫老师的示例文件。这里直接使用data()提取即可
data(ps)网络分析一网打尽
为什么说是一网打尽呢，我使用了network函数打包了ggclusterNet封装的各种功能，一般情况下网络所需要的内容都基本上囊括了，一条代码即可分析完成网络，并且可以按照分组同时分析完成多组网络。完成多组网络的共同比对和输出。
network 函数
计算相关矩阵生成边和节点文件保存为Gephi格式；
使用sne中十多种layout方法并用ggplotv出图，组图
计算节点性质
计算网络性质
计算随机网络和样本网络比对，检测
计算zipi
分组网络计算自动保存 先看看出来的结果
函数测试
这个例子表示的是提取微生物相对丰度在千分之一以上的OTU，，然后对其做 皮尔森相关（推荐使用斯皮尔曼相关），相关阈值选择大于0.6的，并且显著性阈值低于0.05的相关关系，并且使用ggrapg可视化网络，通过igraph计算网络节点和整体属性，后计算zipi模块化信息，将全部结果保存在path中。
```{R} path = “./result” dir.create(path) result = network(ps = ps,N = 0.001,r.threshold=0.6,p.threshold=0.05,label = FALSE,path = path ,zipi = TRUE)
除了保存的结果，这里还有全部分组的网络可视化输出，共比较，和全部网络全局参数一起。
```{R}
# 全部样本的网络比对
p = result[[1]]
p
# 全部样本网络参数比对
data = result[[2]]
plotname1 = paste(path,"/network_all.jpg",sep = "")
ggsave(plotname1, p,width = 22,height = 22)
plotname1 = paste(path,"/network_all.pdf",sep = "")
ggsave(plotname1, p,width = 22,height = 22)
tablename <- paste(path,"/co-occurrence_Grobel_net",".csv",sep = "")
write.csv(data,tablename)结果展示-网络可视化
结果展示-Gephi输入的两个文件
上面的可视化如果你不喜欢，就可以在gephi中出图，使用下面的这两个文件。
KO_Gephi_edge.csv
KO_Gephi_node.csv
···
结果展示-网络性质同随机网络比对
KO_net_VS_erdos_properties.csv
WT_net_VS_erdos_properties.csv
···
结果展示-节点属性
KO_node_properties.csv
WT_node_properties.csv
···
结果展示-zipi计算结果和图形
KOZiPi.csv
KO_ZiPi.pdf
结果展示-幂律分布
这里看到不同网络的幂律分布是不同的，和随机网络之间的差距也不是很相同，这表明了不同分组的网络不一定都是无标度网络。
二分网络一网打尽
二分网络的介绍在之前我已经讲过了，参见推送：
点击跳转。
或者在专辑中查看。
微生物结合环境因子等其他指标网络
这里是我构造的一组模拟环境因子数据 环境因子的格式和vegan中的一样，也就是样本作为列，指标作为行。
—导入多个环境变量及其分组
data1 <- read.delim(“..//ori_data/bios_network_normal//env.txt”) head(data1) row.names(data1) <- data1SampleID = NULL data1 <- as.data.frame(t(data1)) head(data1) data1$ID = row.names(data1) data1 <- data1 %>% select(ID,everything())
Gru <- read.delim(“..//ori_data/bios_network_normal//env_group.txt”)细菌真菌phyloseq对象我们已经准备好了，这里可以参照学习， 如果不会构建phyloseq对象，可以参见之前的推文。```{R}# 导入细菌ps，通过相对丰度value1来过滤otu表格ps16 = readRDS("..//ori_data/bios-network_micro/16S/ps_16S.rds")#导入真菌ps，通过相对丰度value1来过滤otu表格psIT = readRDS("..//ori_data/bios-network_micro/ITS/ps_ITS.rds")ps.merge <- merge16S_ITS (ps16,psIT,N16s = 0.005,NITS = 0.005)ps.merge
同样是一条函数即可计算得到类似于上面的全部内容，由于zipi的计算不够快，所以我这里设置：zipi = FALSE。下面是结果文件的总览，少了zipi的部分：
library(ggpubr)
path = "./result_bio/"
dir.create(path)
result <- corBionetwork(ps = ps.merge,
N = 0.001,
r.threshold = 0.9, # 相关阈值
p.threshold = 0.05,
label = FALSE,
group = "Group",
env = data1, # 环境指标表格
envGroup = Gru,# 环境因子分组文件表格
# layout = "fruchtermanreingold",
path = path,# 结果文件存储路径
fill = "group", # 出图点填充颜色用什么值
size = "igraph.degree", # 出图点大小用什么数据
scale = TRUE, # 是否要进行相对丰度标准化
bio = TRUE, # 是否做二分网络
zipi = FALSE, # 是否计算ZIPI
step = 100, # 随机网络抽样的次数
width = 12,
height = 10
)p = result[[1]]
p
# 全部样本网络参数比对
data = result[[2]]
plotname1 = paste(path,"/network_all.jpg",sep = "")
ggsave(plotname1, p,width = 22,height = 22)
plotname1 = paste(path,"/network_all.pdf",sep = "")
ggsave(plotname1, p,width = 22,height = 22)
tablename <- paste(path,"/co-occurrence_Grobel_net",".csv",sep = "")
write.csv(data,tablename)结果展示
这里我就简单展示一下出图，使用的ggclusterNet的函数出图，展示不同模块之间的关系：
添加主编微信 加入群聊
关于微生信生物 你想要的都在这里
这里的数据我放到了github上：https://github.com/taowenmicro/ggCluster_decument（包括全部流程和结果和原始数据文件）
猜你喜欢
10000+：菌群分析 宝宝与猫狗 梅毒狂想曲 提DNA发Nature Cell专刊 肠道指挥大脑
文献阅读 热心肠 SemanticScholar Geenmedical
16S功能预测 PICRUSt FAPROTAX Bugbase Tax4Fun
生物科普: 肠道细菌 人体上的生命 生命大跃进 细胞暗战 人体奥秘
写在后面
为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外5000+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。PI请明示身份，另有海内外微生物相关PI群供大佬合作交流。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍未解决群内讨论，问题不私聊，帮助同行。
学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”