Sloan中性群落模型(NCM)推断群落构建原理及其R实现

许多SCI文章都使用了一种叫NCM的方法来对微生物群落构建（中性理论 or 生态位理论？）进行预测，本文将对NCM的实现原理及其R实现做介绍。在介绍其原理的过程中会有一些数学公式，不要对这些数学公式产生抵触情绪，我会解释其每一项的生态学含义。与原理相比，其R实现则相对简单，阅读完后可将原理与R实现进行对比学习。

道阻且长1994

12680人浏览 · 2022-03-23 15:48:16

道阻且长1994 · 2022-03-23 15:48:16 发布

1.中性理论原理介绍

In order for an assemblage to change, an individual must die or leave the system. It is then immediately replaced by an individual either by a birth from within the community or by an immigrant from a source community

只有当物种死亡或离开这个系统时，群落结构才会发生改变。此时，离开个体的生态位就会空余出来，其它个体会通过来自群落外的迁移或群落内部的繁殖来填补空出的生态位。因此可以把群落的动态描述为死亡——繁殖/扩散——死亡这样的循环。
根据上述的原理，如果我知道 t0 时刻：

物种 i 的在本地群落（local community）内的个体数为 $N_i$ ，即OTU table（行名：sample ID，列名：species ID）中的每个单元格的数字。
物种 i 在宏群落（meta-community）内的相对丰度为 $p_i$ ，即OTU table中列的和/OTU表的和
该本地群落（local community）总的个体数为 $N_T$ ，即OTU table中行的和
个体死亡后被宏群落中的物种替换的概率为 $m$
物种 i 的繁殖率与所有物种的平均繁殖率的比值 $\alpha_i$

当群落中1个个体（individual）死亡时，我们记此时的时间为 t1。此时会有以下几种可能：

可能1： 死亡的是物种 i 。其概率为: $N_i/N_T$
- 可能1.1： 物种 i 的个体死亡后空余的生态位被群落内的物种通过繁殖的方式取代。其概率为: $1 - m$
  - 可能1.1.1： 被群落内相同的物种取代。其概率为: $(1+\alpha_i)(N_i-1)/(N_T - 1)$
  - 可能1.1.2： 被群落内不相同的物种取代。其概率为: $(1-\alpha_i)(N_T - N_i)/(N_T - 1)$
- 可能1.2： 物种 i 的个体死亡后空余的生态位被群落外的物种通过扩散的方式取代。其概率为: $m$
  - 可能1.2.1： 被群落外相同的物种取代。其概率为: $p_i$
  - 可能1.2.2： 被群落外不相同的物种取代。其概率为: $1 - p_i$
可能2： 死亡的不是物种 i。其概率为: $N_T - N_i)/N_T$
- 可能2.1： 非物种 i 的个体死亡后空余的生态位被群落内的物种通过繁殖的方式取代。其概率为: $1 - m$
  - 可能2.1.1： 被群落内的物种 i 取代。 $(1+\alpha_i)N_i/(N_T - 1)$
  - 可能2.1.2： 被群落内的非物种 i 取代。 $(1-\alpha_i)(N_T - N_i - 1)/(N_T - 1)$
- 可能2.2： 非物种 i 的个体死亡后空余的生态位被群落外的物种通过扩散的方式取代。其概率为: $m$
  - 可能2.2.1： 被群落外的物种 i 取代。 $p_i$
  - 可能2.2.2： 被群落外的非物种 i 取代。 $1 - p_i$

那么，与 t0 时刻相比 t1 时刻物种 i 的数量可能：减少 1个、不变和增加一个，即 t1 时刻物种 i 的数量可能为： $N_i-1$ 、 $N_i$ 和 $N_i+1$ ，其可能的概率如下：

物种 i 增加一个的概率 = 可能2 * （可能2.1 * 可能2.1.1 + 可能2.2 * 可能2.2.1）

物种 i 不变的概率 = 可能1 * （可能1.1 * 可能1.1.1 + 可能1.2 * 可能1.2.1）+ 可能2 * （可能2.1 * 可能2.1.2 + 可能2.2 * 可能2.2.2）

物种 i 减少一个的概率 = 可能1 * （可能1.1 * 可能1.1.2 + 可能1.2 * 可能1.2.2）

写成数学表达式就是文献中的公式了：
在这里插入图片描述
但我认为文献中的公式有点小问题，我认为应该是下面的公式，与文献的公式有些出入,我用红色标出来了：

$Pr(\frac{N_i+1}{N_i})=(\frac{N_T-N_i}{N_T})[mp_i+(1+\alpha_i)(1-m)(\frac{N_i}{N_T-1})$
$Pr(\frac{N_i}{N_i})=\frac{N_i}{N_T}[mp_i+\textcolor{red}{(1+\alpha_i)}(1-m)(\frac{N_i-1}{N_T-1})]+(\frac{N_T-N_i}{N_T})[m(1-p_i)+\textcolor{red}{(1-\alpha_i)}(1-m)(\frac{N_T-N_i-1}{N_T-1})]$
$Pr(\frac{N_i-1}{N_i})=\frac{N_i}{N_T}[m(1-p_i)+(1-\alpha_i)(1-m)(\frac{N_T-N_i}{N_T-1})]$

这个公式的意义就是当我们知道 t0 时刻的物种 i 的丰度时，我们就可以推断下一时刻 t1 时物种 i 的可能的数量。此时，令 $x_i=N_i/N_T$ ，及 $x_i$ 表示物种 i 的相对丰度。根据forward Kolmogorov equation， $x_i$ 在 t 时刻的概率密度函数为： $x_i$ 的概率密度 $\phi_i(x_i,t)$ 随时间 t 的变化为：

在这里插入图片描述
$M_{\delta x}$ 和 $V_{\delta x}$ 中 $M$ 表示期望， $V$ 表示方差， $\delta x$ 表示x的变化量。因此 $M_{\delta x}$ 表示物种 i 相对丰度变化量的期望， $V_{\delta x}$ 表示物种 i 相对丰度变化量的方差。因此：

$M_{\delta x}=1 * Pr(\frac{N_i+1}{N_i})+(-1)*Pr(\frac{N_i-1}{N_i})$

$V_{\delta x}=1 * Pr(\frac{N_i+1}{N_i})+1*Pr(\frac{N_i-1}{N_i})$

可得：
在这里插入图片描述
这里 $V_{\delta x}$ 与 $M_{\delta x}$ 相比可以忽略不计，因为 $N_T$ 在OTU表中每个样点往往有成千上万条序列，即 $N_T>1000$ ，而 $0<x_i<1$ 以 $0<p_i-x_i<1$ 。

且 $\frac{ \partial \phi }{\partial x} = 0$ ，获得 $\phi_i$ 的结果为

在这里插入图片描述
该分布是一个beta分布的概率密度函数。物种 i 占居的频率（occurrence frequency: row sums of binary OTU table/number of sites）为其概率密度函数的积分。

在这里插入图片描述

此时该分布是一个beta分布，我们就可以在R语言中利用beta分布对其进行拟合，获得参数m的评估值。

2. R语言实现

##Fits the neutral model from Sloan et al. 2006 to an OTU table and returns several fitting statistics. Alternatively, will return predicted occurrence frequencies for each OTU based on their abundance in the metacommunity
#Install the following packages if they haven't been availabled in your computer yet 
library(Hmisc)
library(minpack.lm)
library(stats4)
#using Non-linear least squares (NLS) to calculate R2:
#spp: A community table with taxa as rows and samples as columns
spp<-read.csv('spp.txt',head=T,stringsAsFactors=F,row.names=1,sep = "\t")
spp<-t(spp)
#get the mean of abundance of each sample
N <- mean(apply(spp, 1, sum))
#get the mean of species relative abundance in metacommmunity
p.m <- apply(spp, 2, mean)
p.m <- p.m[p.m != 0]
#get the percentage of each species in metacommmunity
p <- p.m/N
#get the binary data of community abundance matrix
spp.bi <- 1*(spp>0)
#get the frequncy of species occurrence in metacommunity
freq <- apply(spp.bi, 2, mean)
freq <- freq[freq != 0]
#get a table record species percentage and occurrence frequency in metacommunity
C <- merge(p, freq, by=0)
#sort the table according to occurence frquency of each species
C <- C[order(C[,2]),]
C <- as.data.frame(C)
#delete rows containning zero
C.0 <- C[!(apply(C, 1, function(y) any(y == 0))),]
p <- C.0[,2]
freq <- C.0[,3]
names(p) <- C.0[,1]
names(freq) <- C.0[,1]
d = 1/N
##Fit model parameter m (or Nm) using Non-linear least squares (NLS)
m.fit <- nlsLM(freq ~ pbeta(d, N*m*p, N*m*(1 -p), lower.tail=FALSE),start=list(m=0.1))
m.fit #get the m value
m.ci <- confint(m.fit, 'm', level=0.95)
freq.pred <- pbeta(d, N*coef(m.fit)*p, N*coef(m.fit)*(1 -p), lower.tail=FALSE)
pred.ci <- binconf(freq.pred*nrow(spp), nrow(spp), alpha=0.05, method="wilson", return.df=TRUE)
# get the R2 value
Rsqr <- 1 - (sum((freq - freq.pred)^2))/(sum((freq - mean(freq))^2))
Rsqr

#data visulization
#Drawing the figure using grid package:
#p is the mean relative abundance
#freq is occurrence frequency
#freq.pred is predicted occurrence frequency
bacnlsALL <-data.frame(p,freq,freq.pred,pred.ci[,2:3])
inter.col<-rep('black',nrow(bacnlsALL))
inter.col[bacnlsALL$freq <= bacnlsALL$Lower]<-'#A52A2A'#define the color of below points
inter.col[bacnlsALL$freq >= bacnlsALL$Upper]<-'#29A6A6'#define the color of up points
library(grid)
grid.newpage()
pushViewport(viewport(h=0.6,w=0.6))
pushViewport(dataViewport(xData=range(log10(bacnlsALL$p)), yData=c(0,1.02),extension=c(0.02,0)))
grid.rect()
grid.points(log10(bacnlsALL$p), bacnlsALL$freq,pch=20,gp=gpar(col=inter.col,cex=0.7))
grid.yaxis()
grid.xaxis()
grid.lines(log10(bacnlsALL$p),bacnlsALL$freq.pred,gp=gpar(col='blue',lwd=2),default='native')

grid.lines(log10(bacnlsALL$p),bacnlsALL$Lower ,gp=gpar(col='blue',lwd=2,lty=2),default='native') 
grid.lines(log10(bacnlsALL$p),bacnlsALL$Upper,gp=gpar(col='blue',lwd=2,lty=2),default='native')  
grid.text(y=unit(0,'npc')-unit(2.5,'lines'),label='Mean Relative Abundance (log10)', gp=gpar(fontface=2)) 
grid.text(x=unit(0,'npc')-unit(3,'lines'),label='Frequency of Occurance',gp=gpar(fontface=2),rot=90) 
#grid.text(x=unit(0,'npc')-unit(-1,'lines'), y=unit(0,'npc')-unit(-15,'lines'),label='Mean Relative Abundance (log)', gp=gpar(fontface=2)) 
#grid.text(round(coef(m.fit)*N),x=unit(0,'npc')-unit(-5,'lines'), y=unit(0,'npc')-unit(-15,'lines'),gp=gpar(fontface=2)) 
#grid.text(label = "Nm=",x=unit(0,'npc')-unit(-3,'lines'), y=unit(0,'npc')-unit(-15,'lines'),gp=gpar(fontface=2))
#grid.text(round(Rsqr,2),x=unit(0,'npc')-unit(-5,'lines'), y=unit(0,'npc')-unit(-16,'lines'),gp=gpar(fontface=2))
#grid.text(label = "Rsqr=",x=unit(0,'npc')-unit(-3,'lines'), y=unit(0,'npc')-unit(-16,'lines'),gp=gpar(fontface=2))
draw.text <- function(just, i, j) {
  grid.text(paste("Rsqr=",round(Rsqr,3),"\n","Nm=",round(coef(m.fit)*N)), x=x[j], y=y[i], just=just)
  #grid.text(deparse(substitute(just)), x=x[j], y=y[i] + unit(2, "lines"),
  #          gp=gpar(col="grey", fontsize=8))
}
x <- unit(1:4/5, "npc")
y <- unit(1:4/5, "npc")
draw.text(c("centre", "bottom"), 4, 1)