R语言并行计算 deviation of null beta diversity（beta多样性零偏差）

  群落构建分析是微生物生态学分析的重要组成部分，成为目前文章发表的热点技术。之前我们介绍了计算beta-NTI（beta nearest taxon index）来进行群落构建分析（https://blog.csdn.net/weixin_43367441/article/details/118515090）。|beta-NTI| >2说明决定性过程主导，其中beta-NTI >2说明OTU的遗传距离发散，为生物交互作用主导，beta-NTI < -2则说明OUT的遗传距离收敛为环境选择主导。利用beta-NTI推测群落构建有一个前提条件，即系统发育树必须要有遗传信号（phylogenetic signal），但某些功能基因（例如nifH）的遗传发育树往往缺乏遗传信号。没有遗传发育信号情况下，我们又该怎样进行群落构建分析呢？beta多样性零偏差就派上用场了。

1. 如何理解beta多样性零偏差
  beta多样性零偏差的计算过程可以概括为以下三点，1.提出零假设，2.构建零模型，3.观察数据与零模型数据进行比较。与我们经常使用的假设检验（t检验、F检验）本质上是相同的，不同的是假设检验的数据分布是已知的，而群落结构的数据我们需要多次构建零模型来模拟。
  首先我们提出零假设，假设群落构建是完全随机的。然后基于此假设构建零模型，由于群落构建过程完全随机，因此各物种在样点间的分布是完全随机、没有偏好的，将每个物种随机分配到各样点，构成一个零模型，重复1000次。最后将观察到的beta多样性数据与零模型的beta多样性数据进行比较，观察数据与零模型的数据越相似，说明构建过程越随机，示意图如下。

2. 如何并行计算beta多样性零偏差

#threads为拟使用的CPU数，dune为otu table，reps为模拟的次数
nulldiv <- function(dune,reps,threads){
  library(reldist)
  library(vegan)
  library(bipartite)
  library(foreach)
  library(doParallel)
  patches <- nrow(dune)
  ##Calculate beta-diversity for metacommunity
  ### Prepare and calculate abundance beta-null deviation metric
  ## Adjusted from Stegen et al 2012 GEB
  bbs.sp.site <- dune
  rand <- reps
  null.alphas <- matrix(NA, ncol(dune), rand)
  null.alpha <- matrix(NA, ncol(dune), rand)
  expected_beta <- matrix(NA, 1, rand)
  null.gamma <- matrix(NA, 1, rand)
  null.alpha.comp <- numeric()
  bucket_bray_res <- NULL
  
  bbs.sp.site = ceiling(bbs.sp.site/max(bbs.sp.site)) 
  mean.alpha = sum(bbs.sp.site)/nrow(bbs.sp.site) #mean.alpha
  gamma <- ncol(bbs.sp.site) #gamma
  obs_beta <- 1-mean.alpha/gamma
  obs_beta_all <- 1-rowSums(bbs.sp.site)/gamma
  
  ##Generate null patches
  registerDoParallel(cores = threads)
  bucket_bray_res <- foreach (randomize = 1:rand,.combine = "cbind") %dopar% {  
    null.dist = dune
    for (species in 1:ncol(null.dist)) {
      tot.abund = sum(null.dist[,species])
      null.dist[,species] = 0
      for (individual in 1:tot.abund) {
        sampled.site = sample(c(1:nrow(bbs.sp.site)), 1)
        null.dist[sampled.site, species] = null.dist[sampled.site, species] + 1
      }
    }
    ##Calculate null deviation for null patches and store
    null.alphas[,randomize] <- apply(null.dist, 2, function(x){sum(ifelse(x > 0, 1, 0))})
    null.gamma[1, randomize] <- sum(ifelse(rowSums(null.dist)>0, 1, 0))
    expected_beta[1, randomize] <- 1 - mean(null.alphas[,randomize]/null.gamma[,randomize])
    null.alpha <- mean(null.alphas[,randomize])
    null.alpha.comp <- c(null.alpha.comp, null.alpha)
    bucket_bray <- as.matrix(vegdist(null.dist, "bray"))
    diag(bucket_bray) <- NA
    bucket_bray_res <- apply(bucket_bray, 2, FUN="mean", na.rm=TRUE)
    
  } ## end randomize loop
  
  ## Calculate beta-diversity for obs metacommunity
  beta_comm_abund <- vegdist(dune, "bray")
  res_beta_comm_abund <- as.matrix(as.dist(beta_comm_abund))
  diag(res_beta_comm_abund) <- NA
  # output beta diversity (Bray)
  beta_div_abund_stoch <- apply(res_beta_comm_abund, 2, FUN="mean", na.rm=TRUE)
  
  # output abundance beta-null deviation
  abund_null_dev <- beta_div_abund_stoch - mean(bucket_bray_res)
  return(abund_null_dev)
}

2. 多核运行可以节省多少时间？
  以前200个OTUs为例，十核运行的时间为28s，单核运行的时间为247秒，快了近十倍，当数据量大时可以帮我们节约更多的时间，有利于我们平时进行数据的探索。

测试数据及其代码见如下链接，两者内容相同：
面包多： https://mianbaoduo.com/o/bread/YpmTlZlp
CSDN：https://download.csdn.net/download/weixin_43367441/72178211
更多R语言分析微生物生态学的资源可参考如下链接：
https://mbd.pub/o/bread/mbd-YpmTlZpr

参考文献
[1] Differentiating between niche and neutral assembly in metacommunities using null models of b-diversity (https://doi.org/10.1111/oik.02803)
[2] 刘红涛,胡天龙,王慧,张燕辉,郭世伟,谢祖彬.典型水稻土细菌亚类群——泛化种、特化种的群落构建及功能潜力[J/OL].土壤学报:1-13[2022-02-07].http://kns.cnki.net/kcms/detail/32.1119.p.20220120.1356.008.html.