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

library(psych)

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

#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))

}

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

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]

library(parallel)

detectCores()

#24

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