• 软件 | 快速计算网络自然连通度评估群落稳定性

    Fastn-快速计算网络自然连接度

    - 软件介绍 -

    从网络图的邻阶矩阵（0、1矩阵，0代表指标间没有边，1代表指标间有边）中删除节点并求剩余矩阵“特征值”被用来表征网络抗毁性，该方法被用在比较不同组的微生物共现网络图的稳定性中。用法为：删除网络图中的1、2、3、4、5、6、7、......N个节点并求取每次删后矩阵的特征值并做转化，将最终结果称为自然连通度。在实际中，一个网络图包含多个节点，那么如何使得某次删除的节点具有代表性呢？一种方法是根据节点的Degree大小来排序，Degree从大到小依次删除节点并计算剩余节点的邻接矩阵的自然连通度。然而，当网络图中的大部分节点的Degree比较接近甚至相等时，按照Degree排序删除节点的做法将不再有效。实施上，我们可以近似认为用删除节点的方式是模拟环境变化使得群落中某些物种消失（或某些物种之间的关系消失，即删除边）之后评估群落的稳定性如何。那么，我们事先不知道哪些物种会消失，例如，当群落中的N个物种消失时，我们很难知道这些N个物种的组合里具体是哪些物种，因此我们使用随机抽样的方法对节点进行随机抓取并删除的方法来模拟物种的随机消失，即计算随机抓取N个节点并删除之后的自然连通度。为了提高结果的可重复性，我们设定删除N个节点时可以迭代1000次（甚至更多次）并分别计算每次删除之后的自然连通度，最后取均值。当网络图节点数目大时，邻接矩阵特征值的计算量非常大，我们在实际计算过程中发现，当删除N个节点时，其自然连通度与（N-1）和（N+1）的结果很接近，因此，在尽量不改变原始结果的情况下，我们设定了步长（step）。当步长大于1，假设为3，那么删除的节点数将为1、4、7、10......（N+step），在这种设置下，计算所需时间将比原始数据（步长为1）至少减少3倍。
      基于上述内容，我们基于C++语言和调用Eigen库设置了删除节点并计算自然连通度的软件--fastnc（Fast calculation of natural connectivity），同时使用OpenMP并行和mkl库加速。经测评，fastnc计算速度快，并且，尽管我们使用了随机抽样的方法，在1000次迭代后，计算的结果仍然可重复 （具体的自然连通度数值会稍微有差别）。我们在将fastnc代码封装为命令行结构时，参考了fastspar软件（用C++实现sparcc算法的快速软件）。fastnc可通过conda安装：

    conda create -n fastnc -c wqssf102 -c conda-forge -c intel fastnc -y

    - 软件参数 -

    1. fastnc --h
    2. Program: FastNC (use c++ to calculate the natural connectivity)
    3. Contact: Qiusheng WU (565715597@qq.com)
    4.  
    5.  
    6. Usage:
    7.   fastnc [options] --adj_table  --outfile 
    8.   -c , --adj_table 
    9.                 The result from get.adjacency function of igraph package
    10.   -o , -outfile 
    11.                 Result output table
    12. Options:
    13.   -t <float>, -threshold <float>
    14.                 The threshold for deletion of node (default: 0.8)
    15.   -s , -step 
    16.                 The step for deletion of node (default: 1)
    17.   -g <float>, -edge <float>
    18.                 The threshold for deletion of edge (default: 0)
    19.   -n , -number 
    20.                 Number of iterations (default: 1000)
    21.   -j , -job 
    22.                 Number of jobs (default: 4)
    23. Other:
    24.   -h        --help
    25.                 Display this help and exit
    1. -c 为邻接矩阵，必须给出
    2.  
    3.  
    4. -o 为输出文件，必须给出
    5.  
    6.  
    7. -t 为删除物种（节点的阈值），默认为0.8，即从1%删到80%
    8.  
    9.  
    10. -s 为步数，默认为1，假如step为3，那么就是依次删除1个节点，4个节点，7个节点，......,即第N次是删除（N+step）个节点（N>1）
    11. 因此，当step大于1时，计算时间将大大缩短，适用于网络图节点多的情况
    12.  
    13.  
    14. -g 删除边的比例（默认不在删除节点及其对应的边之外再删除边），即，先删掉节点，在剩余矩阵中随机删除边，删除边的情况计算量将特别大
    15.  
    16.  
    17. -n 为迭代次数，默认1000，因为我们不知道环境变化会使哪些物种消失，因此这里采用了随机抽样，建议至少迭代1000次
    18.  
    19.  
    20. -j 线程数，默认为4

    邻接矩阵获得的方法之一是通过R软件的igraph包导出，如：

    1. setwd("D:\\fastnc")
    2. ##读取物种数据
    3. sp <- read.csv("sp.csv",header = T,row.names = 1)
    4. ##sp的数据结构
    5. head(sp,c(2,2))
    6.     g__Candidatus.Udaeobacter g__uncultured__f__uncultured__o__Acidobacteriales
    7. site1               0.001288199                                       0.001502898
    8. site2               0.002008320                                       0.000788983
    9. ##读取分组文件
    10. grp <- read.csv("grp.csv",header = T,row.names = 1)
    11. ##样本分组文件
    12. head(grp)
    13.       sampleid group
    14. site1    site1    AA
    15. site2    site2    AA
    16. site3    site3    AA
    17. site4    site4    AA
    18. site5    site5    AA
    19. site6    site6    AA
    20. ##
    21. sp <- sp[grp$sampleid,]
    22. #######按照分组计算物种之间的关系
    23. library(dplyr)
    24. library(igraph)
    25. ##calcor.R为已经整理好的代码，方便计算相关性
    26. source("calcor.R",encoding = "utf-8")
    27. corres <- calcor(codt = sp,group = grp,r_th = 0.6,p_th = 0.5,type = "cor")
    28. ##提取网络图的邻接矩阵
    29. for (i in unique(grp$group)) {
    30.   adjtab <- as.matrix(get.adjacency(corres[[i]]$gg))
    31.   ####导出01矩阵，然后使用服务器计算
    32.   write.table(adjtab,paste("NC/",i,".txt",sep = ""),quote=F,row.names=F,col.names=F,sep="\t")
    33. }

    Linux环境下计算：

    1. ##激活环境
    2. conda activate fastnc 
    3. ##开始计算，当网络总结点数除以step（即-s参数）小于100时候。不建议使用step，默认即可。
    4. ##因为当step大于1时，软件会按照步长删除节点，步长太大，最后剩余的节点过少，就会影响最终绘图时拟合的规律
    5. time fastnc -c AA.txt -j 28 -n 1000 -s 3 -o AA_res.txt 
    6. #real   2m25.998s
    7. #user   64m28.868s
    8. #sys    0m7.128s
    9. ##
    10. time fastnc -c BB.txt -j 28 -n 1000 -s 3 -o BB_res.txt 
    11. #real   2m30.020s
    12. #user   65m51.200s
    13. #sys    0m8.165s
    14. ##关闭conda环境
    15. conda deactivate
    16. ##在实际跑数据中，step是否使用取决于服务器的配置，其次是节点的数量。

    最终结果如下：

    1. nc_index    del_number
    2. 0.376387425101  0.002493765586
    3. 0.368335959568  0.009975062344
    4. 0.361530920709  0.017456359102
    5. 0.356399175783  0.024937655860
    6. 0.352021378405  0.032418952618

    nc_index为删除M个节点后剩余矩阵的“特征值”，del_number为删除节点所占总节点的百分比。

      最后用R绘图：

    1. library(ggpmisc)
    2. grp1 <- read.csv("NC/AA_res.txt",header = T,sep="\t")
    3. grp1$grp <- "AA"
    4. grp2 <- read.csv("NC/BB_res.txt",header = T,sep="\t")
    5. grp2$grp <- "BB"
    6. ####将每组合并
    7. grpnc <- rbind(grp1,grp2)
    8. ##指定图例的顺序
    9. grpnc$grp <- factor(grpnc$grp,levels = c("AA","BB"))
    10. ggplot(grpnc, aes(del_number, nc_index,color=grp)) +##grp为将多个网络图合并时的分组
    11.   geom_point() +
    12.   geom_smooth(formula = y~x,se = FALSE,method = "lm",show.legend = F)+
    13.   stat_poly_eq(formula = y~x,size=5,family="serif",method = "lm",
    14.                output.type = "numeric",
    15.                parse = T,label.x = 0.3,hjust=0,
    16.                mapping =aes(label = paste("italic(R)^2~`=`","~",sprintf("\"%#.*f\"",3,after_stat(`r.squared`)),"*\"\"*",
    17.                 ifelse(after_stat(`p.value`)<=0.001,"'***'",
    18.                 ifelse(after_stat(`p.value`)<=0.01,"'**'",
    19.                  ifelse(after_stat(`p.value`)<=0.05,"'*'",NA))),"*\", \"*",
    20.                   "italic(Slope)~`=`","~",sprintf("\"%#.*f\"",3,after_stat(b_1)),sep = "")))+
    21.   labs(x="Proportion of removed nodes",
    22.        y="Natural connectivity",
    23.        color="Group")+
    24.   scale_color_manual(values = c("AA"="#7570B3","BB"="#1B9E77"))+
    25.   theme_bw()+theme(panel.grid = element_blank())+
    26.   theme(axis.title = element_text(size = 20,colour = "black"))+
    27.   theme(axis.text = element_text(colour = "black",size = 20,angle = 0,vjust = 0.5))+
    28.   theme(legend.title =  element_text(size = 20,colour = "black"), 
    29.         legend.text = element_text(size = 20,colour = "black"),
    30.         strip.text =  element_text(size = 20,colour = "black"),
    31.         text = element_text(family = "serif"))
    32. ##将结果导出
    33. ggsave("nc.pdf",width = 8,height = 6)
    34. ggsave("nc.tiff",width = 8,height = 6)

    27005f68bedbef24a528910ad06e0c53.png

    从上图可知，AA组自然连通度下降趋势快于BB组。

    - 计算速度与时间 -

    1、我们建议将网络图的节点控制在500个以下，网络节点太多导致计算量太大，服务器算不过来。  

     2、若用户的网络包含的节点比较多（如几千个节点），假如用户想要迭代1000次，当服务器的运行内存不大时，我们建议分为4个任务提交到服务器中，每个任务执行250次的迭代，最后取4个任务计算结果的均值即可，这样会大大缩减用户等待的时间。特别地，可以通过设置step参数（如step设置为3）来扩大删除节点的步长，这样耗时将显著减少。总之，step是一个很有用的参数。我们的软件还在开发中，若你对软件的功能有需求或发现软件问题，请联系作者：565715597@qq.com
      fastnc软件目前只在Linux系统下可以正常运行。我们已经将软件源码上传到github上：

    https://github.com/wqssf102/fastnc

    推荐使用conda安装

    src目录为软件的源码，源码需要编译才能使用，用户若想自己编译，那么需要配置依赖的库和软件,编译方式如下所示：

    1. MKLROOT=/opt/intel/oneapi/mkl/2022.1.0
    2. ##
    3. g++ -std=c++11 -O3 -fopenmp -march=native -mavx -mfma -o fastnc fastnc.cpp fastnc_opts.cpp common.cpp -DMKL_ILP64 -m64\
    4.  ${MKLROOT}/lib/intel64/libmkl_scalapack_ilp64.a \
    5.  -Wl,--start-group ${MKLROOT}/lib/intel64/libmkl_cdft_core.a \
    6.  ${MKLROOT}/lib/intel64/libmkl_intel_ilp64.a \
    7.  ${MKLROOT}/lib/intel64/libmkl_gnu_thread.a \
    8.  ${MKLROOT}/lib/intel64/libmkl_core.a\
    9.  ${MKLROOT}/lib/intel64/libmkl_blacs_openmpi_ilp64.a \
    10.  -Wl,--end-group -lgomp -lpthread -lm -ldl

    参考文献

    [1] Wu MH, Chen SY, Chen JW, Xue K, Chen SL, Wang XM, Chen T, Kang SC, Rui JP, Thies JE, Bardgett RD, Wang YF. Reduced microbial stability in the active layer is associated with carbon loss under alpine permafrost degradation. Proc Natl Acad Sci U S A. 2021 Jun 22;118(25):e2025321118. doi: 10.1073/pnas.2025321118.

    [2] G.-s. Peng, J. Wu, Optimal network topology for structural robustness based on natural connectivity. Physica A 443, 212–220 (2016). doi: 10.1016/j.physa.2015.09.023

    - 作者简介 -

    吴求生

    中科院地化所在读博士生

    研究方向为环境生物地球化学，具备一定的扩增子、宏基因组上下游数据分析经验，导师为袁权研究员。

    猜你喜欢

    iMeta简介 高引文章 高颜值绘图imageGP 网络分析iNAP
    iMeta网页工具 代谢组MetOrigin 美吉云乳酸化预测DeepKla
    iMeta综述 肠菌菌群 植物菌群 口腔菌群 蛋白质结构预测

    10000+：菌群分析 宝宝与猫狗 梅毒狂想曲 提DNA发Nature

    系列教程：微生物组入门 Biostar 微生物组  宏基因组

    专业技能：学术图表 高分文章 生信宝典 不可或缺的人

    一文读懂：宏基因组 寄生虫益处 进化树 必备技能：提问 搜索  Endnote

    扩增子分析：图表解读 分析流程 统计绘图

    16S功能预测   PICRUSt  FAPROTAX  Bugbase Tax4Fun

    生物科普:  肠道细菌 人体上的生命 生命大跃进  细胞暗战 人体奥秘  

    写在后面

    为鼓励读者交流快速解决科研困难，我们建立了“宏基因组”讨论群，己有国内外6000+ 科研人员加入。请添加主编微信meta-genomics带你入群，务必备注“姓名-单位-研究方向-职称/年级”。高级职称请注明身份，另有海内外微生物PI群供大佬合作交流。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍未解决群内讨论，问题不私聊，帮助同行。

    a5f80204f850f22f0bc79594b31995f5.jpeg

    点击阅读原文，跳转最新文章目录阅读

  • 相关阅读:
    稀土掺杂氟化物纳米荧光微球/稀土荧光磁性纳米微球Fe3O4@PHEMA-RE的制备方法
    Moment.js 如何对时间进行比较获得不同的天数
    逗号分隔String字符串 - 数组 - 集合，相互转换
    SRC中逻辑漏洞检查总结
    深入理解JavaScript堆栈、事件循环、执行上下文和作用域以及闭包
    SQL 文本函数
    计算机基础--Git
    【LeetCode】No.88. Merge Sorted Array -- Java Version
    前端裸辞躺平三个月自己的一点想法哈哈哈
    mongdb迁移方案及比对方案
  • 原文地址：https://blog.csdn.net/woodcorpse/article/details/127699024