R语⾔并⾏计算spearman相关系数
利⽤spearman相关性分析是构建共现⽹络的重要⽅法，但由于OTU table往往有成千上万⾏，⽤R⾃带的st()函数计算较为费时，严重制约我们的分析速度。对spearman相关性分析进⾏并⾏化运⾏可⼤⼤节省计算时间，为此我们⼿写了spearman相关性分析函数来实现并⾏化运⾏。为⽅便讲解，本⽂以我们常见的OTU table 数据为例，对OTU进⾏两两spearman相关性分析，获得相关系数r和显著性p值。我们将⾃⼰⼿写的函数network_construct()与psych包中的st()函数两者运⾏时间和计算的结果进⾏了⽐较，我们⾃⼰的函数network_construct()计算时间远远少于st()函数且结果相同，具体的R代码见下⽂。
我们经常使⽤st()函数计算OTU之间的相关性，但该函数在⾯对较多的OTU时速度较慢。`
library(psych)
进修自我鉴定system.st(otu[,1:500],u="pairwi",method="spearman",adjust="fdr",alpha=.05))
运⾏时间如下：（elapd栏为程序运⾏的时间）
为了加快计算速度，我们根据相关系数矩阵的计算特点，进⾏了并⾏化重写，代码如下：
#otu_table ⾏名为样点名，列名为OTU 名
#threads 为使⽤的CPU 核数
#本函数默认使⽤‘fdr’进⾏P 值矫正，可参考st()函数的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 <- i (pp )
pa <- pp [lp ]
pa <- p.adjust (pa , "fdr")
pp [i (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核进⾏计算，运⾏时间如下：
我们可以看到使⽤st()函数进⾏计算需要313秒，⽽使⽤R代码进⾏并⾏计算仅需0.99秒。 我们再⽐较⼀下不同线程情况下所需的时间。单核运⾏代码如下:
女生动漫头像黑白system.time (network_construct (otu [,1:500],1))
............
即使使⽤单核进⾏计算速度，也要⽐st()函数快上很多，毕竟该函数⾥⾯有很多判断语句， 计算更加费时。当我们加⼤数据量时，单核和多核的区别就更加明显。下图是单核与10核运⾏的⽐较：
我们可以看到计算量越⼤，多核的性能就越优越。当我们有成千上万个OTU时可以节省很多时间。 那么我们运⾏的结果与st()函数的结果是否⼀致呢？
res1 <- st(otu[,1:50],u="pairwi",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包中的st()的结果进⾏⽐较，结果如下：
我们计算的结果与原函数的结果完全相同，可以放⼼使⽤哈，现在就可以把我们的函数复制过去⽤于⾃⼰的项⽬啦。
除了共现⽹络分析，群落构建分析中的零模型构建也需要进⾏⼤量的循环，我们也可以将并⾏计算⽤于群落构建分析，例如beta-NTI的计算，敬请关注我们下⼀篇⽂章。
温馨提⽰：
1. 我们⾃⼰构造的函数network_construct()依赖“foreach”,“doParallel”,“abind”程序包，请预先安装上述三个包。st()
依赖于“psych”程序包。
2. 运⾏前可通过如下代码查看计算机CPU核数，核数不是越多越好，要根据内存来选，不然内存不够会让R程序卡死。可通过任务管理
器实时监测。
重启电脑3. 测试⽂件可参考付费资源，
4. 更多R语⾔并⾏计算分析微⽣物⽣态学的资源可参考如下链接：
library(parallel)
主席台座次安排图detectCores()
#24
#以我的电脑为例是⼆⼗四核
