mbio: 随机采样过程会高估微生物群落的beta多样性

2013年发表于mbio。是微生物生态领域非常经典和重要的一篇文章！

β多样性，即两地点之间物种组成的差异，对于理解物种多样性的时空格局以及控制群落组成和结构的机制至关重要。

然而测序技术中的测序错误、偏差、重现性和定量能力低导致定量β多样性十分困难，且会极大地高估群落的多样性。

本文针对群落大小已知和未知两种情况，基于一般采样理论(general sampling theory)分别提出了数学框架，模拟随机采样对β多样性的影响。

此数学框架能够准确预测技术重复之间低重现性的问题，表明随机采样过程是导致技术重复之间变化的主要因素。

另外，根据泊松随机抽样(Poisson random sampling)预测的值与观测到的OTU overlap相一致，进一步证明了技术重复之间重现性低是由于随机采样过程中的人为因素导致。

最后，建立了一个数学框架来预测一定的采样量可以达到的样本之间重叠的期望值。

利用本文的模型进行预测，为了达到技术重复之间高度的重现性，目前的测序量还需要再提升好几个数量级。

由于大多数的生态学研究都包含随机采样过程，本研究的结论也可以推广到其他的生态研究中。

基于扩增子测序的方法重现性低，尤其是对低丰度的物种。技术重复之间变异大的原因可能是由于随机采样过程导致。

由于微生物群落中大部分物种丰度很低，因此对他们的检测就十分困难。在此基础上，resample过程想要得到低丰度物种就更困难了。

采样带来的偏差取决于群落的复杂度和样本量。群落越复杂偏差会越严重；而通过增加样本量来减少偏差还未有理论支持。

假设：随机采样过程是导致技术重复之间高度变化的主要原因。

证明：先构建了一个理论框架模拟随机采样过程，并预测为了达到理想的重现性所需要的采样量（sampling efforts）；其次利用一个抽球的例子解释随机采样过程带来的偏差，同时检验了这个框架能否应用于预测不同技术重复之间OTU overlap低的情况。

Mathematical framework

(i) Sampling individuals from a large regional community.

大尺度上的物种丰度分布与小样本观测到的物种丰度之间的关系至关重要。本文采用一般采样理论模拟随机采样过程的偏差。

假设一个物种在样本中的个体数取决于物种在大群落中的丰度、采样量、以及个体的空间分布。

先假设个体都是随机分布。N表示总个体数（如16S rRNA gene sequences）, n为物种数量（如OTUs）,丰度分别为x1, x2, . . ., xn。

随机采样得到一个个体，属于第i个物种ith的概率为xi/N。

从群落中随机采样得到m个个体，则属于第i个物种ith的个体数为mxi/N。

那么，在样本m中，得到k个个体中属于第i个物种ith的概率符合二项分布。

若k=0，

令xi/N = axi/m,其中 a = m/N, 表示采样比例,（2）可以近似为Poisson分布的指数形式：

这里先复习一下高等数学：

其中

显然/易证/易得/易知，(3)式成立。

因此，至少存在一个个体属于i个物种ith的概率为

Poisson分布是最简单的采样模型。根据一般采样理论，样本中观测到的丰度分布可以表示为

φa(m)是采样量为m时观测到的物种丰度分布；

φ(x)表示群落中丰度为x的物种丰度分布；

θ是参数向量。

(ii) Expected species overlap among samples with the size of the large community known.

群落大小N已知，两个和三个样本overlap的物种数分别表示为

其中a1,a2,a3分别为三个样本占总体的比例，即a1 = m1/N，a2,a3类似。

//////////

Jaccard和Bray-Curtis用来计算不同样本之间的overlap。

基于Jaccard相似性指数，两个样本共有的OTU（OJ2a1,a2,θ）计算如下：

其实就是P = P(AB)/(P(A∪B)-P(A∩B))

三个样本类似，公式略过。

(iii) Expected species overlap among samples with the size of the large community unknown.

群落大小N未知，计算期望的物种overlap。

大多数情况下群落个体数未知。即N未知。

泊松采样理论条件下，物种的丰度分布遵循尺度不变性，因此样本的丰度分布（y）可以通过群落丰度分布(x)进行估计。y = px，p为群落中被采样的比例。因此公式8可以被写为：

a1*, a2*, and θ*是N未知时的采样比例和参数向量。

随机采样两次，令p = a1 + a2，a1* = a1/(a1+a2), a2* =a2/(a1+ a2)。

m1和m2是两个样本中观测到的总个体数。所以a1=m1/N, a2=m2/N。a1*=m1/(m1 + m2), a2* =m2/(m1+ m2)。

因此物种overlap可以根据y而不是x得到，这样就不需要知道N。

(iv) Predicting sampling efforts for achieving a desired overlap among replicate samples.

假设两个样本分别需要采集m1' 和m2'个个体，保证所需要的overlap。为了简化计算，令m1'= m2'= m'。采样比例A= m'/N。

N已知时，预测的两样本overlap是：

N未知时，y= px = (a1 + a2)x。

N已知或未知时，共有物种数的公式可以统一为：

令A* = m=/(m1 + m2)，预测的overlap模型是：

结果

采用下图的方法模拟随机采样过程。三个罐子完全相同。N个球n种颜色。不同颜色的丰度分布完全相同。一个球就代表一个16S rRNA序列，一种颜色代表一个OTU。随机抽m个球出来，若m=N,三个样本的overlap为100%。但是实际上overlap取决于采样量、球的丰度分布、群落的复杂性。由于球的丰度分布相同，只有随机采样过程会给结果带来影响。

通过模拟实验与理论公式得到的结果进行比较，发现两者没有明显的差异，证明了模型的准确性。且N已知未知也没有显著差异。

模拟实验与理论公式的卡方检验结果。运用五种不同的物种丰度分布，在不同样本个数（2，3）和N已知与未知条件下检验拟合度。

蓝线为N已知条件的公式计算；红线为N未知的公式计算。点为实际观测值。A为2样本模拟，B为3样本模拟。理论和实际符合的非常好。N已知未知无显著差异。

样本量低时技术重复之间差异很大，但是随着样本量增加很快能达到平台期。

要达到一定的overlap所需要的序列数

A为2样本，B为3样本。5万条序列期望达到的overlap都高于80%。

对于2个样本，要达到90%的overlap理论需要71400条序列；

对于3个样本，要达到90%的overlap理论需要63770条序列。

目前绝大多数的微生物研究，尤其是土壤群落相关的研究，采样量都远远上面的理论值。

增加样本量和增加（生物学和技术）重复是最有效的改善重复之间重现性低的方法。