R语言并行计算加快土壤微生物生态学分析(1):spearman相关系数加快共现网络构建速度

2021
09/10

+
分享
评论
微生态
A-
A+

共现网络 (co-occurrence network) 和群落构建 (assembly process) 分析日益成为微生物生态学分析中重要的组成部分,



本文由微科盟胡天龙根据实践经验而整理,希望对大家有帮助。

微科盟原创微文,欢迎转发转载,转载须注明来源《微生态》公众号


   

共现网络 (co-occurrence network)   群落构建 (assembly process) 分析日益成为微生物生态学分析中重要的组成部分,成为目前文章发表的热点技术。这些分析中往往会有大量的循环语句,如果顺序执行的话往往非常耗时。R语言的并行计算可以有效地解决耗时的问题,加快程序运行的速度。本文主要介绍如何利用R语言并行计算spearman相关系数,加快共现网络(co-occurrence network)构建速度。后期的的文章将逐一介绍群落构建分析的并行计算,包括R语言并行计算beta-NTI (beta-nearest taxon index)、β多样性零偏差 (null deviation of beta diversity) 和Raup-Crick距离 (Raup-Crick dissimilarity)
利用spearman相关性分析是构建共现网络的重要方法,但由于OTU table往往有成千上万行,用R自带的corr.test()函数计算较为费时,严重制约我们的分析速度。对spearman相关性分析进行并行化运行可大大节省计算时间,为此我们手写了spearman相关性分析函数来实现并行化运行。为方便讲解,本文以我们  常见的OTU table 数据  为例(联系您所添加的任一  微科盟组学老师  即可  免费  领取  ),对OTU进行两两spearman相关性分析,获得相关系数r和显著性p值。我们将自己手写的函数network_construct()与psych包中的corr.test()函数两者运行时间和计算的结果进行了比较,我们自己的函数network_construct()计算时间远远少于corr.test()函数且结果相同,具体的R代码见下文。





运行步骤


我们经常使用corr.test () 函数计算OTU之间的相关性,但该函数在面对较多的OTU时速度较慢。


library(psych)

system.time(corr.test(otu[,1:500],use="pairwise",method="spearman",adjust="fdr",alpha=.05))

运行时间如下:(elapsed栏为程序运行的时间)


图1 原函数运行时间


为了加快计算速度,我们根据相关系数矩阵的计算特点,进行了并行化重写,代码如下:

#otu_table行名为样点名,列名为OTU名

#threads为使用的CPU核数

#本函数默认使用‘fdr’进行P值矫正,可参考corr.test()函数的R帮助文档

network_construct <- function(otu_table,threads){ 

 library(foreach)  

 library(doParallel)  

 library(abind)

  otu_table2 <- apply(otu_table,2,rank)  

   r <- function(rx,ry){  

    n <- length(rx)  

    lxy <- sum((rx-mean(rx))*(ry-mean(ry)))  

    lxx <- sum((rx-mean(rx))^2)  

    lyy <- sum((ry-mean(ry))^2)   

    r <- lxy/sqrt(lxx*lyy)   

    t <- (r * sqrt(n - 2))/sqrt(1 - r^2)   

    p <- -2 * expm1(pt(abs(t), (n - 2), log.p = TRUE))    

   return(c(r,p))  

  }  

 arraybind <- function(...){    

    abind(...,along = 3,force.array=TRUE)  

  } 

  nc <- ncol(otu_table) 

  registerDoParallel(cores = threads)  

  corr <- foreach (i = 1:nc,.combine = "arraybind") %dopar%{   

    corr1 <- matrix(rep(0,2*nc),nrow = 2,ncol=nc)   

    for(j in 1:nc) {     

       if(j > i) corr1[,j] <- r(otu_table2[,i],otu_table2[,j])  

    }   

      corr <- corr1  

  }  

  rr <- corr[1,,] 

  rr <- rr+t(rr)  

 diag(rr) <- 1 

 pp <- corr[2,,] 

 lp <- lower.tri(pp) 

 pa <- pp[lp]  

 pa <- p.adjust(pa, "fdr")  

 pp[lower.tri(pp, diag = FALSE)] <- pa  

 pp <- pp+t(pp)  

 rownames(pp) <- colnames(otu_table) 

 colnames(pp) <- colnames(otu_table)  

 rownames(rr) <- colnames(otu_table) 

 colnames(rr) <- colnames(otu_table)  

 return(list(r = rr,p = pp))

}


使用我们自己构造的函数后,可利用10核进行计算,运行时间如下:


图2 自建函数十核运行时间


我们可以看到使用corr.test () 函数进行计算需要313秒,而使用R代码进行并行计算仅需0.99秒。我们再比较一下不同线程情况下所需的时间。单核运行代码如下:


system.time(network_construct(otu[,1:500],1))

 

图3 自建函数单核运行时间


即使使用单核进行计算速度,也要比corr.test () 函数快上很多,毕竟该函数里面有很多判断语句, 计算更加费时。当我们加大数据量时,单核和多核的区别就更加明显。下图是单核与10核运行的比较:


 图4 自建函数单核与十核运行时间的比较


我们可以看到计算量越大,多核的性能就越优越。当我们有成千上万个OTU时可以节省很多时间。那么我们运行的结果与corr.test () 函数的结果是否一致呢?


res1<- corr.test(otu[,1:50],use="pairwise",method="spearman",adjust="fdr",alpha=.05)

res2 <- network_construct(otu[,1:50],10)

res1$r[1:5,1:5]

res2$r[1:5,1:5]

 

将我们构建的函数network_construct () 与psych包中的corr.test () 的结果进行比较,结果如下:


图5 自建函数与原函数结果比较


我们计算的结果与原函数的结果完全相同,可以放心使用哈,现在就可以把我们的函数复制过去用于自己的项目啦。
除了共现网络分析,群落构建分析中的零模型构建也需要进行大量的循环,我们也可以将并行计算用于群落构建分析,例如beta-NTI的计算,敬请关注我们下一篇文章。

温馨提示


1. 我们自己构造的函数network_construct () 依赖 “foreach”, “doParallel”, “abind”程序包,请预先安装上述三个包。corr.test () 依赖于“psych”程序包。
2. 运行前可通过如下代码查看计算机CPU核数,核数不是越多越好,要根据内存来选,不然内存不够会让R程序卡死。可通过任务管理器实时监测。


library(parallel)

detectCores()

#24

#以我的电脑为例是二十四核


 

本文来源于微科盟原创作者胡天龙,仅用于学术分享,如有侵权,请联系删除!




本文由“健康号”用户上传、授权发布,以上内容(含文字、图片、视频)不代表健康界立场。“健康号”系信息发布平台,仅提供信息存储服务,如有转载、侵权等任何问题,请联系健康界(jkh@hmkx.cn)处理。
关键词:
table,R语言,生态学,微生物,计算,函数,运行,代码

人点赞

收藏

人收藏

打赏

打赏

我有话说

0条评论

0/500

评论字数超出限制

表情
评论

为你推荐

推荐课程


社群

精彩视频

您的申请提交成功

确定 取消
剩余5
×

打赏金额

认可我就打赏我~

1元 5元 10元 20元 50元 其它

打赏

打赏作者

认可我就打赏我~

×

扫描二维码

立即打赏给Ta吧!

温馨提示:仅支持微信支付!