编译：小白同学，编辑：小菌菌、江舜尧。

原创微文，欢迎转发转载。

1 引言

测序技术和分子生物学方法的进步推动了微生物组的研究，从采样设计到生物学解释都有所进步。然而，由于测序数据可能存在测定误差且需要环境因子等辅助数据，用不恰当的方法分析土壤微生物可能得到错误的推论。此外，在整个过程中（DNA提取、扩增、测序和生物信息学分析），随着系统误差和随机误差的累积也会造成结果的偏差。观察到的序列丰度间差异并不仅仅是样本中DNA真实差异，所以不能从原始数据的比较中得出结论。另一方面，土壤微生物间存在复杂的相互作用。一旦某一物种相对丰度升高，必然导致其他物种丰度降低，在统计学上表现为显著的负相关却未必存在相关关系，所以不能仅用相关性来表示这些相互作用。尤其在群落中有优势类群时，这种假阳性相关关系的问题尤为突出。考虑到这些方法和挑战，本文回顾了当前用于微生物群落分析的主要方法，旨在整合基于模型的分析方法，主要目的在于强调正确理解各个方法前提假设的重要性，尤其是：

（i）应用统计学方法分析土壤微生物群落，

（ii）这些方法切实可行的应用场景，

（iii）以及使用各方法所需的假设检验。

此外，通过缺陷检出的方法或可处理测序偏好性，有待进一步开发以应用于群落分析。通过分析土壤微生物组，不仅应得到一些推断性的信息，也应该找到误差来源和处理方法。为进行解释说明，这里用gjam R包中的树苗内生微生物组数据和GOMMS包中的模拟数据进行举例。

2 扩增子测序数据中能有何发现？

微生物生态学主要研究微生物与环境之间相互作用关系，应用分子生物学方法主要围绕以下五个主题进行研究：

（i）评估土壤微生物多样性（即系统发育多样性）；

（ii）识别优势物种和稀有微生物；

（iii）探讨实验处理或地理距离对土壤微生物群落的影响；

（iv）识别对土壤因子敏感的微生物类群；

（v）构建微生物类群之间及其内部的关系以了解群落结构。

高通量测序技术的发展对这些研究大有裨益，另一方面，海量信息也需要更合适的数据分析方法，以整合并解读这千百种微生物的测序结果。起初，数据处理方法源于经典的群落生态学：

（i）评估土壤微生物多样性的alpha和beta多样性指数等；

（ii）用于识别优势物种和稀有物种的相对丰度和对数比等度量方法；

（iii）比较土壤微生物群落

（iv）通过多变量分析（聚类分析、主成分分析等）方法确定敏感微生物类群。

主要思路是通过基于距离的排序进行降维，将群落数据以散点图等方式进行简洁明了的可视化，随后可以进行基于置换的分析（PERMANOVA、ANOSIM、Manteltest等）。排序后也能够识别相关性，但开始主要采用更直接的相关性分析（如Spearman或Pearson）探讨微生物类群间关系并进行网络分析。不同指标在衡量类群间相互作用时具有相似性，明确用什么指标量化微生物群落很有必要，但还未有确切的方案对其进行指导。度量标准的选择通常有赖于方法所提供的潜在信息，如Bray-Curtis距离反映了两样本中属于不同类群生物的百分比。基于所提供的潜在信息，可能有人在不知道该指标是否适用于物种组成数据的情况下就用了这种方法。为了能选择合适的方法，首先需要了解测序数据的特征，及现有数据分析方法的局限性和可能造成的偏差。

3 传统分析方法受限于测序数据的特性

测序数据不同于宏观生物在传统生态学中的数据：首先，土壤样本中的基因序列与计数数据不同，基于距离的方法如果考虑比例则容易产生误差；此外，宏基因组工作流程的每一步（DNA提取、PCR和测序）都会积累特异性误差。为了控制测序偏好性的影响，微生物组数据通常将序列条数除以总数转换成比例，或者其他度量。另一种常用方法通过‘抽平’使样品序列条数一致，代价是损失了大量的信息。这些方法的目的是为了满足经典的多元分析和统计检验，使用随机置换进行推理分析。因此，这些方法侧重于调整数据以适应特定类型的分析，而不是使用合适的分析方法来适应数据的类型。

3.1 微生物群落数据的关键特性

测序之前，各个实验环节中都存在潜在误差来源，可能会对所得测序结果有一定影响。不同土壤类型的DNA提取效率有差异；16S扩增子测序受引物选择所限，能识别的物种数量远少于鸟枪法测序，且会放大样本间的多样性；序列聚类方法也颇受质疑。鉴于测序平台和分子方法的局限性，OTU数据不能反映微生物群落丰度，而只能鉴定微生物门类。从本质上说，测序数据过度分散，且每个样本的测序深度都不相同，这些特性源于高通量测序带来的非独立性。每个测序平台都有内存大小的限制，无论DNA样本中实际数据量有多大，输出的数据量将取决于测序平台的内存限制。因此，就土壤微生物学研究而言，数据分析应该关注比例而非数量上的差异。如果少数类群在土壤群落中占主导，则可能导致显著的偏倚，这需要更先进的统计分析方法来适应数据的这种特性。

微生物组数据也具有稀疏性的特征，零计数可能是真实缺失或未检测到所致，分为三种类型：

（i）未被检测到的采样零点；

（ii）DNA提取方法、PCR或测序平台的限制等引起的技术零点；

（iii）样品中确不存在的生物学零点。

丢弃少于某一阈值的OTU可能会带来偏差，特别是对于alpha多样性的分析。总之，通过高通量测序检测到的微生物丰度远非真实情况。考虑由均值-方差关系所衡量的类群间异质性时，异质性随着数据稀疏性的增加而增加，产生所谓的过度分散，影响样品点在分布图中的坐标位置和结果的解读。因此，均值-方差关系在多元分析中十分重要，若被忽略，对微生物群落的分析将出现假阴或假阳性的误判。图1为内生微生物组数据集均值-方差间的显著正相关关系。

图1 植物内生微生物组数据集均值-方差间的显著正相关关系。

3.2 先验假设

基于距离的度量都隐含了关于数据变异性的假设，这增加了选择合适度量方法的难度。例如，只有当各类群变异性相似时，Bray-Curtis距离才适用。因为这种方法最初是为相似区域或相似容量的调查而设计的，需要均匀采样。若采样不均，就如同测序深度不同，那么Bray-Curtis距离将无法识别微生物群落间的差异，除非舍弃大量信息。不恰当的差异度量方法将无法得到数据真实的均值-方差关系，而这是多种群数据的一个关键的统计学属性。基于距离的分析方法倘若选择不当，会影响对样本间差异的区分。如果几组样品的变异性不同（方差不齐），那么非度量多维尺度分析NMDS就难以将真正相似的样品聚类到一起以区分组间差异，而主坐标分析（PCoA）往往展示的是虚假差异。因此，有必要明确选用的分析方法是否适合所研究的数据

另一个问题是分析方法对于假设条件的灵活性，比如一些非参数的统计学方法，这种灵活性可能使作为重要前提的检验方法的假设条件被忽略。例如，PERMANOVA被用来比较两个或两个以上处理组之间的差异，这种方法利用置换和随机过程来提供不依赖于数据正态分布的稳健的分析。然而，许多差异度量方法对样本组间的差异很敏感。因此，为避免造成有偏的PERMANOVA，应首先通过同质性检验来鉴定数据组间离散度。由此发展而来的PERMDISP检验可用来评估各组样本内的离散性或同质性。在许多已发表的文献中，这一关键步骤经常被忽略。进行相关性分析前也需要检验样本独立性的假设，由于测序所得物种丰度的数据不具备独立性，难以用Pearson相关分析确定两组分间的相关性。相关性分析还涉及如何将实际生物作用从其他影响因素中分离出来的问题，这将在4.4小节中进行讨论。

3.3 数据转换

开始处理数据时可能需要对数据进行转换，以减少均值-方差关系不一致的影响。考虑到微生物群落数据的稀疏性，可将数据进行对数转换，但这种方法存在效率问题，数据转换不见得能稳定数据的异质性并满足线性模型的假设。因此，当数据转换或归一化成为误差来源而非助力于公平比较时，就会强化分析方法的偏差。最常用的归一化方式—相对丰度（TSS），无法控制过度分散，且破坏了数据的独立性，增加了相关性是假阳性的可能。其实高通量测序数据本身的特性就造成了数据独立性的缺失，由于测序存在一定偏好性，需要通过将原始数值除以总数得到特定百分比来消除偏差（CSS）。通过模型还可以区分真正的生物学零点和技术零点，当一些样品中序列条数很少时很可能是各实验环节中测定技术的问题造成的技术零点。CSS的主要优势在于对数变换，若对TSS数据进行一个简单的对数变换也可以产生类似的效果（图2）。尽管将比例数据对数化后线性相关，但仍存在稀疏性的问题，或可考虑矩阵补全等方法进行进一步完善。此外，为应对数据较高的异质性，广义线性模型（GLM）也是常用的数据处理方法。GLM会选出数据的实际分布类型，而非通过调整数据改变其分布情况来适应分析方法。但GLM仍然依赖于模型中提出的假设，而这些假设并非广为人知，所以用的时候出错的风险更高。

图2 对不同植物内生微生物组数据进行MDS降维分析，（a）用CSS进行标准化的结果（b）用TSS对数转换标准化的分析结果。

4 土壤微生物组的无偏分析方法

合适的分析方法应该能够规避各种数据特性可能引入的偏差，应满足：

（i）基于物种组成分析微生物群落间的差异，

（ii）明确测序偏好性对数据变异性的影响，

（iii）适应数据的均值-方差关系，

（iv）在框架内加入环境因子等额外信息，

（v）考虑稀有物种的影响，

（vi）从环境响应的影响中分离出生物间相互作用。

群落分析不应过度依赖于相似性或差异性，而应关注分析方法能否有效处理高通量测序数据并检验模型充分性。有效的微生物群落分析方法应不受数据均值-方差关系的影响，目前来看，广义线性模型可以通过线性预测和分布假设将真实效应从随机噪声中分离出来，是处理测序数据的最佳选择。

4.1 基于模型的方法

广义线性模型将数据分为观察变量（如OTU）和解释变量（或预测因子，如温度、pH等环境参数），GLM需要三个基本条件：

（i）观测变量的概率分布，

（ii）两个或两个以上的解释变量产生的综合效应（如在养分充足的条件下，加热可能使某些微生物增多），

（iii）拟合模型提供解释变量作用下的期望值。

综上所述，GLM对观测变量与解释变量之间的线性关系进行定量描述，变量数据遵循某种特定分布，由此可将观测数据和解释变量联系起来。GLM允许对数据的概率分布进行选择，从而能够在特定的均值-方差关系下进行分析，进而可以检验GLM是否适合要分析的数据。为了进一步进行说明，这里用泊松分布对过度分散数据进行模型拟合并进行残差分析。如果数据过度分散，那么泊松分布模型的残差将随微生物类群增多而持续增加（图3a），其他辅助分析允许我们根据样本（图3b）和物种分类单元（图3c）检查残差同质性。分位数图（图3d）中，如果不是直线，数据偏离假定分布，说明先前的选择可能不适合给定的数据集。在基于模型的方法中，假设检验应该是强制性的；分析应运用先验知识和诊断工具先抓数据关键属性，再考虑如何用这些属性进行合理推断。将GLM结合矩阵回归扩展到差异分析，广义差异模型（GDM）可以将采样点和环境条件间差异进行关联，或可解决基于距离度量差异性的位置分散偏差。

图3 假定数据遵循泊松分布时微生物群落的残差图，（a）线性解释变量（b）样品（c）物种类别（d）分位数图。

4.2 基于模型的排序分析

排序分析（如PCA、PCoA和NMDS）从具有许多变量的高维数据中提取出两个或三个响应变量。基于模型的方法也可以进行降维，如带有潜变量的GLM或潜变量模型（LVM），或可克服测序数据本身的局限性。与经典排序方法相似，潜变量模型构建了潜因子的坐标轴来模拟数据变异性，且使用能有效捕捉数据变异性的特定概率分布。同样地，后面通过残差分析判断何为最佳选择。广义潜变量模型可应用于非限制性排序和带有环境因子的限制性排序分析。尽管输入数据为序列条数，LVM可以通过整合位点效应来衡量样品中物种相对丰度，所以不需要进行数据转换。为了进行解释说明，这里用已发表的模拟数据进行LVM分析，图4a中数据变异性未被潜变量充分捕获，而考虑位点效应后（图4b），排序偏差消失，组间差异更加明显。LVM还可以通过在预测因子（协变量）中包含多项式来处理复杂关系：例如土壤pH过高或过低都不利于微生物生存，而是存在一个最优值，LVM在模型方程中加入一个二次项（x²）来整合这种关系。与其它统计方法一样，LVM也是建立在前提假设的基础上的，建模前应确保满足假设条件。如果选择的概率分布不能正确捕获数据变异性，在分析期间就能看出来。即便有时候考虑和不考虑位点效应之间（图4 c-d）差异不大，在微生物组的分析中仍要考虑这个问题。更重要的是，基于模型的排序方法为许多不同数据类型提供了统一的分析框架。LVM可以整合物种间相互作用及其他相关性信息，从而可以实现节点物种分布模型的开发。

图4 潜变量模型（LVM）。

4.3 节点物种分布模型

节点物种分布模型（JSDM）不仅考虑了物种间相互作用，还确定了协变量对物种丰度的影响。然而，JSDM在微生物生态学领域还未被充分应用。最近，层次贝叶斯模型（HBM）进一步通过残差项来分析物种间相关性并引入结构误差对微生物群落数据进行解读。残差项（未被实验处理或协变量解释的变异性）不仅可能是随机噪声，还可能反应出未考虑到的影响因素。HBM的应用构建了物种群落层次模型（HMSC）的框架，将层级节点物种分布与多种影响因子相结合，探讨它们之间的相互作用关系。由此，可以度量微生物群落中可被（土壤理化性质、空间距离、时间等）解释的变异的百分比。图5为利用HMSC评估植物内生真菌数据中能被宿主和温度等其他协变量解释的变异，其中宿主对微生物起主要选择作用，温度也解释了约22.6%的变异，此外也有部分未能解释的变异，这些说明了不同的生态过程是如何塑造微生物群落的。

图5 利用HMSC评估植物内生真菌数据中能被宿主和温度等其他因子解释的变异。

4.4 微生物共现性

土壤中的微生物群落作为一个集体进行物质代谢和能量交换，从而对环境产生影响，所以大尺度上对物种分布和群落结构的研究需要考虑多种微生物间相互作用。变量间的相关性并不意味着内在的因果关系（除非解释变量为实验处理）。统计学的相关性分析有助于进行相关推断和假设，但实际的生物学相关性还需要通过进一步研究（如全基因组测序或培养方法）来验证。在进行网络分析时不应该涵盖那些稀有物种，数据稀疏性越强（0值占比越大），出现假阳性的可能性越大，关联分析和网络分析不适用于稀有物种。因此，在进行分析前需要对数据进行过滤，以适应数据分析方法对数据的要求并得到可靠的分析结果。不管采用何种方法，相关分析本质上都会混淆真正的生物相互作用和两微生物对同一环境因子的共同响应。针对这个问题，JSDM可以区分到底观测到的微生物共现性是源于响应环境变化，还是由于环境因子解释范围外的未知变量（如微生物间相互作用），这样可以避免由共同响应环境因子带来的微生物共现性混淆视听。最后，要将基于模型的方法应用于微生物组数据分析，应注意两个问题：

（1）模型能否较好地拟合数据，如何评估? 

（2）如果模型拟合良好，又是否有效地总结了原始数据?

对数据进行生物学解释前要明确这些问题以保障结果解读的可靠性。

5 展望与结论

测序数据本身和分析方法都是有局限性的，测得的物种丰度信息存在不确定性，低检出的常见物种和高检出的稀有物种不应赋以相同权重。基于模型的方法可以通过从观测丰度中剔除噪声干扰，优化出实际的物种分布情况。这需要对样品进行技术重复，例如对一份样品多次提取DNA、对一份DNA进行多次PCR等，并考察该过程引入的变异性，通过检测记录可获知各物种的检出情况。检测过程中各种来源的偏差可以整合到层次模型中，确定各物种的检出率后可以对实验设计进行分析。未来可通过模拟方法来探讨分析模型的缺陷，并与其他方法进行比较，以期完善测序数据对微生物群落的解读。土壤微生物组的高维数据结构带来了新的科学问题和挑战，解决方法的选择和应用要依赖假设检验作为先决条件。基于模型的方法可以应对数据的关键统计学特性，如缺乏独立性等。通过明确均值-方差关系，潜变量模型可以整合位点效应并评估样品中物种的相对丰度。然而，基于模型的方法通常比基于距离的方法需要更多的计算资源和时间。为提高群落分析的准确性，数据分析模型既要能测算物种丰度，也要考虑检测偏差、环境因素的影响、及物种间相互作用等。概率统计方法在分析微生物组数据时可以指导科学推理，明确结果解读的不确定性。

你可能还喜欢

1. 2019年度回顾 | 微生态环境微生物类微文大合辑
   

2. 2019年度回顾 | 微生态人体/动物微生物类微文大合辑

3. 2019年度回顾 | 技术贴合辑大放送
   

微生态科研学术群期待与您交流更多微生态科研问题

（联系微生态老师即可申请入群）。

了解更多菌群知识，请关注“微生态”。

阅读原文

微信扫一扫
关注该公众号