CN110021338A - 一种基于Hill numbers与时间关系的幂法则模型预测微生物群落结构的方法 - Google Patents

Abstract

本发明公开了一种定量监测微生物群落结构变化趋势的方法，具体指基于微生物群落多样性Hill numbers指数与时间关系的幂法则标度模型来监测或预测群落多样性随时间变化的方法。根据被测区域内微生物群落多样性(Hill numbers)指标，通过建立“幂法则标度模型”和“指数截断的幂法则标度模型”，分别构建累积群落多样性指数与累积采样时间之间的数学模型，并根据模型参数获得一套评估体系，从而能够有效地评估或预测群落多样性随时间的变化。该方法可用于监测或预测自然界中各种(包括人体、食品、水质、土壤等)微生物生态群落多样性随时间的变化，特别是人体微生物菌群，从而为菌群相关疾病的个性化精准诊治以及健康菌群的监测和预测提供定量化指标。

Description

一种基于Hill numbers与时间关系的幂法则模型预测微生物 群落结构的方法

技术领域

本发明涉及计算生物学、医学及生态学领域，特别涉及一种监测或预测微生物群落结构随时间变化趋势的方法，具体指利用微生物物种多样性Hill numbers指数与时间关系的幂法则标度模型来监测群落多样性随时间变化的方法。

背景技术

微生物是生态系统重要的组成之一，它们无处不在，数量庞大，在生态系统中发挥着重要功能。特定区域内微生物群落的稳定性牵涉到该区域整个生态系统的稳定性，如河流、湖泊和温泉等水体内的微生物与水质密切相关，而土壤中微生物的组成会影响土壤的营养。微生物群落的物种数量会随着时间进行波动，这种时间动态广泛存在于各个特定位点的微生物群落当中。而不同群落的动态模型又不尽相同，该不同体现在微生物群落对外界环境或人为干扰的响应程度。如新生婴儿肠道和新叶表面部分微生物群落，或水体系统中随季节更替的微生物群落，会直接随着外界环境的改变而改变。然而，如废水处理系统中微生物群落却极少出现时间波动，特别是其中的核心种群表现出较强的时间稳定性。对微生物群落的时间动态进行监控，可以更好的预测外界干扰对群落的影响，从而进行相应的人工干预来维护生态系统的稳定。

人体共生着大量的微生物，除少数病原菌外，它们与人体互利共生，相互影响。大量研究表明，人体微生物菌群与人体的健康和疾病有着紧密的联系。宿主生活习惯、健康状况的变化直接或间接地影响着特定部位菌群的结构，如结肠癌、牙周炎、细菌性阴道病、皮炎及肥胖等疾病还被称为“微生物菌群相关疾病”。而不同时期的结肠癌患者肠道内微生物的物种组成及多样性也发生了改变。人体共生菌群的时间动态可以在一定程度上反应人体的健康和疾病状况，为疾病的预防、诊断和治疗提供全新的视角。

在生态学领域中，研究物种随时间变化的情况主要从物种丰富度(Richness) 与时间的关系着手，即物种-时间关系(Species-Time Relationship，STR)。其中，物种丰富度指一个群落内物种的数目。作为重要的生态学法则，STR描述了在观测的时间范围内，一个群落的物种丰富度随时间增长的情况。然而，传统的STR模型仅能反应一个群落物种丰富度随时间的改变，忽略了其他重要的群落结构特征，如每个物种的丰度信息。物种丰度是衡量物种自身个体数量的指标，物种丰度及物种丰度分布情况在一定程度影响着群落结构的稳定。STR的片面性限制了STR的应用范围，因此需要一个更为可靠、综合的评估体系对群落结构进行时间序列上的监测。

发明内容

本发明的目的在于针对现有技术方法的不足，提供一种基于“多样性Hillnumbers指数一时间”(Diversity-Time Relationship，DTR)的幂法则标度模型来监测自然界中各类微生物群落结构变化的方法。DTR模型系STR模型的全新扩展STR仅限于物种丰富度，而DTR则不但涵盖了对物种丰富度的检测，还可以用于群落中不同水平物种的多样性检测。该方法适用范围广，可用于检测各生态系统中微生物菌群多样性随时间的变化情况，尤其可用于人体微生物菌群相关疾病的监测、诊断和干预，为个性化精准医疗的诊断和防治提供可靠的理论基础和定量化指标。

为实现以上目的，本发明采用的技术方案为：

步骤【1】数据采集：一定时段内对特定区域的微生物生态群落进行监测，获得一组不同监测时间点下各个物种的丰度信息数据。例如，通过16s-rRNA宏基因测序获得各个样点所有物种的丰度。

步骤【2】计算多样性指数(Hill numbers)：物种多样性作为群落结构的重要属性之一，它除了包含物种丰度的信息外，还反应了群落物种个体数分布的均匀程度。物种多样性包括alpha和beta两种，前者度量单个群落的多样性，后者度量群落间的差异。

首先，按照采样时间点顺序对物种丰度信息进行累积。接着，利用Hill numbers多样性计算公式计算各累积采样时间下微生物群落的alpha物种Hill numbers指数，公式(1)所示：

在该式中^qD表示q阶下的alpha累积物种Hill numbers指数，S为各累积采样时间下微生物物种的数量，p_i代表各累积采样时间中每个物种的相对丰度(Relative Abundance)，q为Hill numbers的阶数，在本发明中选取了4个指数阶数：q＝0，1，2，3。当q＝0，2， 3时alpha累积物种Hill numbers可由公式(1)获得。当q＝1时，公式(1)无法定义，故用公式(2)进行计算：

其次，对微生物群落数据按照采样时间点(依次为前2个采样时间点，前3次采样时间点，......，前n次采样时间点)的顺序计算beta物种Hill numbers，计算公式如 (3)所示：

^qD_β＝^qD_γ/^qD_α (3)

其中，^qD_β表示物种的q阶下beta物种Hill numbers。^qD_α表示q阶下的alpha物种Hillnumbers，计算公式如(1)(2)所示，其中p_i代表各时间点中每个物种的相对丰度。^qD_γ为q阶下gamma物种多样性，gamma物种多样性：

其中，p_i+代表各累积采样时间中每个物种对应的相对丰度。

其中，当q＝0时，各物种的相对丰度并未参与计算，此时⁰D代表的是群落中物种的丰富度；当q＝1时，¹D相当于Shannon熵的指数，代表群落中典型或常见物种的 alpha物种多样性；当q＝2时，²D等于Simpson指数的倒数，此时物种的丰度在Hill numbers的计算中占据重要的地位，所得结果代表群落中占优势物种的多样性。

步骤【3】模型构建：基于幂法则标度(PL：Power Law Scaling Model)构建微生物DTR的数学模型，如(5)所示：

^qD＝cT^w (5)

其中，^qD为步骤【2】所得的q阶物种Hill numbers指数，T代表累积采样时间，c和w为模型参数。通过对数转换，可得：

ln(^qD)＝ln(c)+wln(T) (6)

公式(6)可用来评估PL-DTR模型中的各个参数，其中w即为评估体系中的w_pl。当w_pl＞ 0时，表示物种多样性随时间呈增长趋势，当w_pl＜0时，表示物种多样性随时间呈减少的趋势。

根据参数w_pl可以求出物种多样性在不同时间尺度间的重叠程度，本发明的评估指标g为时间尺度T与时间尺度2T间的多样性重叠程度，将参数w_pl代入公式(7)计算各多样性阶数下的g值：

g＝(2D_T-D_2T)/D_T＝2-2^w (7)

若g＝0，则两个时间尺度内的群落物种多样性没有重叠，若g＝1，则表明两个时间尺度内的物种多样性完全重合，当g不等于0或1时，则代表两个时间尺度的物种多样性均有重叠。

PL模型针对q＝0时的物种多样性效果较好，而面对群落中典型物种或常见物种的多样性时效果欠佳，为弥补该缺陷，本发明还扩展出第二个DTR模型。该模型基于指数截断的幂法则标度(PLEC：Power Law with Exponential Cutoff)，数学模型如(8) 所示：

^qD＝cT^wexp(dT) (8)

其中，^qD代表步骤【2】中计算得到的阶数q对应的物种Hill numbers指数，T代表累积采样时间，c，w和d为模型参数。通过对数转换，可得：

ln(^qD)＝ln(c)+wln(T)+dT (9)

公式(9)可用来评估PLEC-DTR模型中的各个参数，其中w即为评估体系中的w_plec。当w_plec＞0时，表示物种多样性随时间呈增长趋势，当w_plec＜0时，表示物种多样性随时间呈减少的趋势。

PLEC-DTR模型在面对典型物种和常见物种多样性时优于PL-DTR模型，此外，该模型还提供了一个重要的评估指标MAD(Maximal Accrual Diversity)，该指标预被监测群落在监测时间段内“群落最大累积多样性”，公式如(10)所示：

其中，Max(^qD)即为评估指标MAD。

根据上述计算，可获得被监测区域在监测时间段内群落结构随时间变化的 DTR模型。根据该模型的参数可获得完整的评估体系，如表1所示，具体包括：该体系包括依据PL-DTR模型参数获得的各多样性阶数(q)下的评估指标w_pl和g，以及依据 PLEC-DTR模型的评估指标w_plec和MAD。

表1.基于DTR模型监测或预测微生物群落结构变化的评估体系*

*PL：Power Law Scaling Model，即幂法则标度模型。

PLEC：Power Law with Exponential Cutoff，即指数截断的幂法则标度模型。

MAD：Maximal Accumulation Diversity，即“群落最大累积多样性”

w_pl：diversity scaling parameter from PL，即“幂法则标度参数”

w_plec：diversity scaling parameter from PLEC，即“指数截断幂法标度参数”

g：pair-wise diversity overlap(similarity)parameter，即“多样性重叠(相似度)参数”

本发明的优势和独特功能在于：用群落物种多样性Hill numbers指数(Community Diversity)代替常规的物种丰富度(Species Richness)来探索群落物种结构随时间变化的情况。突破了现有方法针对群落单一物种结构特征(物种丰度)的局限，从多个层面反应了群落物种多样性随观测时间变化的规律，不仅包含了现有方法涉及的群落物种丰富度，还涵盖了群落中典型物种和优势物种多样性等多个群落物种结构。本发明适用范围广，可用于监测或预测各种微生物生态菌群结构或多样性随时间的变化情况，尤其可用于人体微生物菌群的监测和预测，包括人体微生物菌群相关疾病的评估、诊断和干预，以及为个体话精准医疗的诊断和防治提供可靠的理论基础。

具体实施方式

下面结合具体实施例对本发明做进一步的说明，但并不局限于此。

实施例1：监测一名健康男性肠道微生物菌群结构在一年内的变化

步骤【1】数据来源(Citation of the data source)：

该肠道微生物菌群数据采自一名健康男性。在第0天到第364天期间，对该名男性的粪便进行采样，每个时间点获取一个样本，一年内共得到340个时间点的肠道微生物菌群样本。原作者采用了16s-rRNA基因标记测序技术获得了该名男性肠道微生物菌群随时间变化的丰富度(Richness)和丰度(Abundance)数据。

步骤【2】计算多样性指数：首先，按照采样时间顺序对该名健康男性的肠道菌群的物种丰度进行累积，获得累积采样时间对应的微生物群落物种累积丰度的信息数据，并计算340个采样时间点下群落的累积alpha物种Hill numbers指数。

接着，计算该肠道菌群的beta物种Hill numbers指数。计算过程为：依次按照前2个采样时间点，前3个采样时间点，......，前n个采样时间点的顺序分别计算群落的beta物种Hill numbers，最终获得该肠道菌群的beta物种多样性Hill numbers指数。其中， alpha和beta物种Hill numbers指数均选取4个多样性阶数：q＝0，1，2，3。

步骤【3】模型构建：基于PL构建肠道菌群物种Hill numbers指数与采样时间关系PL-DTR数学模型，模型如下公式所示：

^qD＝cT^w

其中，^qD为步骤(2)所得的q阶Hill numbers指数，T代表累积采样时间，c和w为模型参数。通过对数转换，可得：

ln(^qD)＝ln(c)+wln(T)

公式(6)可用来评估PL-DTR模型中的各个参数，其中w即为评估体系中的w_pl。根据参数w_pl可以求出该名男性肠道菌群物种多样性在不同时间尺度间的重叠程度g，公式如 (7)所示：

g＝(2D_T-D_2T)/D_T＝2-2^w

基于PLEC构建的DTR数学模型，公式如下所示：

^qD＝cT^wexp(dT)

其中，^qD代表步骤(2)中计算得到的阶数q对应的累积Hill numbers指数，T代表累积采样时间，c，w和d为模型参数。通过对数转换，可得：

ln(^qD)＝ln(c)+wln(T)+dT

公式(9)可用来评估PLEC-DTR模型中的各个参数，其中w即为评估体系中的w_plec。基于PLEC-DTR模型参数可求得评估指标MAD，公式如下所示：

其中，Max(^qD)即为评估指标MAD。

根据上述计算，可获得该名健康男性的肠道微生物菌群在一年内的结构或多样性变化，这里以alpha多样性为例示范评估体系的四项指标，如表1所示。其中，w_pl或 w_plec＞0时，表示物种多样性随时间呈增长趋势，w_pl或w_plec＜0时，表示物种多样性随时间呈减少的趋势。若g＝0，则两个时间尺度的群落物种多样性没有重叠，g＝1，则表明两个时间尺度的物种多样性完全重合，当g介于0～1之间时两个时间尺度的多样性有重叠。

表2.一名健康男性肠道微生物菌群一年内的结构变化的DTR模型评估体系*

*PL：Power Law Scaling Model，即幂法则标度模型。

PLEC：Power Law with Exponential Cutoff，即指数截断的幂法则标度模型。

MAD：Maximal Accumulation Diversity，即“群落最大累积多样性”

w_pl：diversity scaling parameter from PL，即“幂法则标度参数”

w_plec：diversity scaling parameter from PLEC即“指数截断幂法标度参数”

g：pair-wise diversity overlap(similarity)parameter即“多样性重叠(相似度)参数”。

Claims (7)

1.一种基于多样性Hill numbers指数和时间关系(Diversity-Time Relationship，DTR)的幂法则标度模型定量预测微生物群落结构变化的方法。

2.按照权利要求1所述的方法，其特征在于：该方法基于不同多样性Hill numbers指数阶数下DTR模型参数构成的评估体系，如表1所示：

表1.基于DTR模型定量监测或预测微生物群落结构变化的评估体系*

*PL：Power Law Scaling Model，即幂法则标度模型

PLEC：Power Law with Exponential Cutoff，即指数截断的幂法则标度模型

MAD：Maximal Accumulation Diversity，即“群落最大累积多样性”

w_pl：diversity scaling parameter from PL，即“幂法则标度参数”

w_plec：diversity scaling parameter from PLEC，即“指数截断幂法标度参数”

g：pair-wise diversity overlap(similarity)parameter，即“多样性重叠(相似度)参数”。

3.按照权利要求2所示的不同多样性Hill numbers指数阶数下DTR模型参数构建的评估体系，其特征在于根据如下方法计算各评估指标：

(1)数据：在一定时段内对特定区域的微生物生态群落进行监测，获得一组不同监测时间点下各物种的丰度信息数据，例如通过16s-rRNA标记的宏基因测序获得物种丰度信息；

(2)计算多样性Hill numbers指数：对物种丰度数据按照采样时间点的顺序进行累积，获得累积采样时间下群落物种累积丰度的信息数据；利用Hill numbers多样性计算公式，选取4个多样性指数阶数(亦可超过4个)：q＝0，1，2，3，分别计算不同阶数下各累积采样时间的alpha累积物种Hill numbers指数和beta物种Hill numbers指数；

(3)模型构建：根据幂法则标度模型(PL：Power Law Scaling Model)获取各Hillnumbers指数阶数下的评估指标w_pl，模型核心如公式如下所示：

^qD＝cT^w

其中，^qD代表步骤(2)中计算得到的阶数q对应的alpha或beta物种Hill numbers指数，T代表累积采样时间，c为模型参数，w为评估指标w_pl。

根据如下公式计算各Hill numbers指数阶数下的群落多样性重叠(相似)度g：

g＝(2D_T-D_2T)/D_T＝2-2^w

其中，w为各Hill numbers指数阶数下的评估指标w_pl。

根据指数截断的幂法则标度模型(PLEC：Power Law with Exponential Cutoff)获取各Hill numbers指数下的评估指标w_plec，模型核心如公式如下所示：

^qD＝cT^wexp(dT)

其中，^qD代表步骤(2)中计算得到的阶数q对应的alpha或beta物种Hill numbers指数，T代表累积采样时间，c和d为模型参数，w为评估指标w_plec。

将c，d和w带入如下公式计算评估指标MAD，即“群落最大累积多样性”(MaximalAccumulation Diversity)：

其中，Max(^qD)即为评估指标MAD。

综上可得完整的DTR评估体系，即表1所示，具体包括：w_pl(diversity scalingparameter from PL)，即“幂法则标度参数”；g(pair-wise diversity overlap(similarity)parameter)，即“多样性重叠(相似度)参数”；w_plec(diversity scalingparameter from PLEC)，即“指数截断幂法则标度参数”；MAD，即“群落最大累积多样性”。

4.按照权利要求1-3所述的方法，其特征在于：该方法用于监测自然界及人类生产生活中各种微生物生态群落多样性随时间的变化规律。

5.按照权利要求4所述的应用，该方法用于食品中及食品生产过程中微生物菌群结构的监测和预测，其特征在于：该方法用于发酵产品中及发酵过程中微生物菌群结构的监测和预测，发酵产品包括酒类、酱油、醋、腐乳、酸奶、泡菜、饲料和发酵型面点等。

6.按照权利要求4所述的应用，该方法用于水体和土壤中微生物菌群结构的监测和预测，其特征在于：水体包括河流、湖泊、蓄水池、水库、废水处理系统、温泉和饮用水等。

7.按照权利要求4所述的应用，该方法用于人体微生物菌群结构的监测和预测。其特征在于：该方法用于人体微生物菌群相关疾病的个性化精准诊断和治疗，以及用于病程的监测和预测，并为人体健康微生物菌群的监测和预测提供定量化指标。