R语言并行计算加快土壤微生物生态学分析（1）：spearman相关系数加快共现网络构建速度
共现网络 (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 自建函数与原函数结果比较
温馨提示
library(parallel) detectCores() #24 #以我的电脑为例是二十四核 |
本文来源于微科盟原创作者胡天龙，仅用于学术分享，如有侵权，请联系删除！
人点赞
人收藏
打赏
010-82736610
股票代码: 872612
京公网安备 11010802020745号
