igraph 算网络指标_ggClusterNet一条代码完成全内容微生物网络

清疚 2023-01-01 10:51 129阅读 0赞

## 写在前面 ##

今年的毕业生应该是6月底离校，但是这位朋友到了九月份终于完善了手上的工作，对这几年来了一个完善的结束，感谢送给我的一方印章，也祝铠明鸣大展宏图。再有就是ggclusterNet功能也完善了90%，剩下的慢慢补充上来，尤其是代码的优化要抽出些时间多和大佬交流，优化。要知道这部分实在是也很头疼。原计划今年将微生物组分析的四个R包够写好，想来想去似乎目前只有easystat和ggClusterNet可以顺利完成，我也会抓紧时间在今年结束的时候至少完成一个包(phyloBiota)功能。![715cc147ef12850f35d8d1f3525fd625.png][]

## 导入所需要的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
 *  分组网络计算自动保存 先看看出来的结果
    
    ![074b9b7d1403cc9a552db56776b4ede9.png][]
    
    #### 函数测试 ####
    
    这个例子表示的是提取微生物相对丰度在千分之一以上的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 write.csv(data,tablename)

#### 结果展示-网络可视化 ####

![926263c197a6aebc905ef64f6e37a7d2.png][]

#### 结果展示-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
    
    ![6dbec30a05d494120b0e57645ac58721.png][]
    
    #### 结果展示-幂律分布 ####
    
    这里看到不同网络的幂律分布是不同的，和随机网络之间的差距也不是很相同，这表明了不同分组的网络不一定都是无标度网络。
    
    ![0ac5f4856a82ff4148e77985117d02fa.png][]![1d928c02e3368f1fbe9fc19a47d4b538.png][]
    
    ## 二分网络一网打尽 ##
    
    二分网络的介绍在之前我已经讲过了，参见推送：
    
    点击跳转。
    
    或者在专辑中查看。
    
    ### 微生物结合环境因子等其他指标网络 ###
    
    这里是我构造的一组模拟环境因子数据 环境因子的格式和vegan中的一样，也就是样本作为列，指标作为行。

# —导入多个环境变量及其分组 #

data1 data1) row.names(SampleID = NULL data1 data1)) head(Gru

细菌真菌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 ps.merge

同样是一条函数即可计算得到类似于上面的全部内容，由于zipi的计算不够快，所以我这里设置：zipi = FALSE。下面是结果文件的总览，少了zipi的部分：![75aa1534f96893c719194ee9c8f3e4d9.png][]

library(ggpubr)
    path = "./result_bio/"
    dir.create(path)
    
    result                         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 write.csv(data,tablename)

#### 结果展示 ####

这里我就简单展示一下出图，使用的ggclusterNet的函数出图，展示不同模块之间的关系：![e7d14ee5a522a725d11ec509979b60dc.png][]

### 添加主编微信 加入群聊 ###

![102995490f09ebd7fd73094d65e5ec66.png][]

### 关于微生信生物 你想要的都在这里 ###

微生信生物

![76c7a487503926aa5e8603a68b6b282d.png][]

这里的数据我放到了github上：https://github.com/taowenmicro/ggCluster\_decument(包括全部流程和结果和原始数据文件)

[715cc147ef12850f35d8d1f3525fd625.png]: /images/20221120/898536ae92bf479caac7dac5f2905f4a.png
[074b9b7d1403cc9a552db56776b4ede9.png]: /images/20221120/265d7a1478194c8c871870b3530e71d5.png
[926263c197a6aebc905ef64f6e37a7d2.png]: /images/20221120/0daba4eaac8044a5a8a3a9b4abbf7105.png
[6dbec30a05d494120b0e57645ac58721.png]: /images/20221120/5d36882108b6477b99e7652105815128.png
[0ac5f4856a82ff4148e77985117d02fa.png]: https://img-blog.csdnimg.cn/img_convert/0ac5f4856a82ff4148e77985117d02fa.png
[1d928c02e3368f1fbe9fc19a47d4b538.png]: /images/20221120/06dab0383f55442da061ff6e6a2b8423.png
[75aa1534f96893c719194ee9c8f3e4d9.png]: /images/20221120/fa8d4d2f8ff74e7b8d76a66ff46f0b49.png
[e7d14ee5a522a725d11ec509979b60dc.png]: /images/20221120/b1b476dee051431da303d2ab377d6002.png
[102995490f09ebd7fd73094d65e5ec66.png]: /images/20221120/88cd4c440e9d4c6da6af3c82358582bc.png
[76c7a487503926aa5e8603a68b6b282d.png]: /images/20221120/c7286542868c4f5181b62dc518cefaf8.png