R语言并行计算spearman相关系数

R语言并行计算spearman相关系数,加快共现网络(co-occurrence network)构建速度

   利用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核运行的比较:
自建函数单核与十核运行时间的比较
   我们可以看到计算量越大,多核的性能就越优越。当我们有成千上万个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程序卡死。可通过任务管理器实时监测。
  3. 测试文件可参考付费资源,https://download.csdn.net/download/weixin_43367441/22467741
  4. 更多R语言分析微生物生态学的资源可参考如下链接:
    https://mianbaoduo.com/o/bread/mbd-YpmTlZpr
library(parallel)
detectCores()
#24
#以我的电脑为例是二十四核
  • 25
    点赞
  • 46
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 8
    评论
评论 8
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

道阻且长1994

您的鼓励有助于我创作更多高质量

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值