微生物操作分类单元确定和序列辅助分离的制作方法

微生物操作分类单元确定和序列辅助分离的制作方法
【专利摘要】基于样品中微生物的系统发生信息基因的DNA序列定义样品中的微生物操作分类单元(OTU)的方法，其通过处理原始序列读出以获得合格序列片段，根据其各自的相对丰度值对其进行排序，并且仅使用具有高丰度值的合格序列获得暂定OTU来进行。重新划归具有低丰度的合格序列，并且只有当所述合格序列与OTU序列具有至少97％的序列相似性时才将所述合格序列分别分配到合适的暂定OTU。本发明还提供了序列辅助的微生物分离方法。
【专利说明】微生物操作分类单元确定和序列辅助分离
【背景技术】
[0001]基于序列相似性或同源性将扩增的DNA序列（例如，16S rRNA扩增子序列）分配到 操作分类单元(operational taxonomic unit，0TU)中是微生物群落研究中的基础方案。 OTU划分(OTU delineation)对于揭示微生物群落的结构和鉴定关键物种是关键性的U，其 可以在下游分析中指导对重要功能细菌的分离和表征 3,4。
[0002] OTU划分方法当作为软件包执行时被称为"流程(pipeline)"。三种常用的流程是 QnME9、M0THUR8和USEARCH7。它们在用经454Illumina测序生成的相同短标签测序数据估计 OTU数目方面显示出不同的结果。Chen等示出，10种已评估的OTU划分方法(Mothur、Muscle+ Mothur、ESPRIT、ESPRIT-Tree、SLP、Uclust、CD-HIT、DNAClust、GramCluster 和 CROP)通常过 高估计了包括43个物种的模拟数据中OTU的数目（1708.5 ± 1386.9)。不同方法也显示出宽 范围的不同：ESPRIT给出了的最大估计OTU数目（4397)，是期望的102.3倍高;而CROP产生最 小的估计OTU数目（133)，仍然是真实数目的3.1倍 5 Jonder等在OTU划分方法（Qiime Blast、CD-HIT、ESPRIT_Tree、Mothur furthest、Mothur average、Uclust、Uclust ref和 Uclust ref optimal)之前对序列进行了去噪和嵌合体检查，但是最低的OTU数目（25,通过 CD-HIT、ESPRIT-Tree和Uclust获得）仍然比具有15个物种的模拟数据中所期望的高出 66.6%6 Jdgar等提出，UPARSE可得到与具有22个物种的模拟数据中的真实计数非常接近 的 0TU，而其他方法(△111口1;[(301^0186、]\1〇1:11111'和(>);[;[1116)将具有1.1±0.8、2.1±1.7和103.0± 36.1倍高的0TU 7。但是Uparse仍然多出1个0TU，该OTU与模拟参照的同一性〈97%。
[0003] 来自相同短标签测序数据的OTU数目的过高估计也存在于Illumina测序中。当通 过模拟样品的三个子区域扩增子测序进行评估时，MOTHUR产生了 20个物种之期望值的2.0 ± 0.1、2.5 ± 0.1和10.1 ± 3.4倍8。通过对具有22个物种的模拟数据的正向末端读出进行 OTU划分，QHME得到高8.4倍的0TU(206与22)，而USEARCH给出2个假OTU(与模拟参照的同 一性〈97% )。此外，当通过Usearch分析经合并的双末端读出时，出现了4.3± 1.3个假0TU7。 因此，所有的三种常用流程QIME9、M0THUR8和USEARCH 7都过高估计了 OTU的数目。
[0004] OTU过高估计产生了许多假0TU，这使微生物群落的组成分布进一步扭曲。在后续 实验中阻碍了重要功能细菌的分离和鉴定。因此，找出这些流程为何产生高数目的假OTU并 开发解决该问题的解决方案很重要。

【发明内容】

[0005] 本申请公开了使虚假(Pseud0)OTU减到最少的改进方法。在该研究中，构建了具有 22种不同16S rRNA基因克隆的7组模拟群落，每个群落的克隆成员浓度不同。在Illumina Miseq平台上对这些群落的16S rRNA基因 V3V4高变区1(3的扩增子进行三次独立批次的测序， 进行批次内重复和批次间重复。先前的研究揭示了原始读出的错误类型 11>12,而我们集中于 经过质量过滤(quality filtration)的"合格序列"，这是OTU划分准确度的直接原因。然后 应用并评估这三种常用的流程:QnME9、M0THUR 8和USEARCH7。追踪每个OTU的详细来源以推 断这些流程为何不同程度地过高估计OTU的数目，并且设计了改进方法以使这些虚假OTU减 到最少。另外，利用具有不同目标区域（V4或V3V4)和测序长度（150bp、200bp、250bpS 300bp)的四个实际数据集，通过测量OTU数目以及alpha和beta多样性的改善证实了该改进 方法。
[0006] 在一个实施方案中，本发明提供了定义样品中的微生物操作分类单元(OTU)的方 法，所述方法包括:获得样品，所述样品包含各自含有系统发生信息基因的微生物;使用基 于PCR的高通量测序技术获得样品中微生物的系统发生信息基因的原始序列读出；处理原 始序列读出以获得合格序列片段;通过处理器获得每条合格序列的相对丰度值，其中所有 合格序列的总相对丰度为100% ;通过处理器根据合格序列各自的相对丰度值对所有合格 序列从高到低进行排序，并将合格序列分成高丰度组和低丰度组，其中所述高丰度组由丰 度值高于低丰度组中合格序列的丰度值并且总体占总丰度约75%的合格序列组成，所述低 丰度组由占总丰度约25 %的剩余的合格序列组成;通过处理器仅使用高丰度组中的合格序 列划分样品中的OTU从而获得暂定(Tentative)OTU;以及通过处理器重新划归(re-map)低 丰度组中的合格序列到暂定0TU，并且只有当合格序列与OTU序列具有至少97%的序列相似 性时才将所述合格序列分别分配到合适的暂定0TU，以实现对OTU的最终定义。
[0007] 在一个实施方案中，系统发生信息基因选自16s rRNA基因或18s rRNA基因。
[0008]在一个实施方案中，系统发生信息基因是16s rRNA基因的一个或更多个可变区， 例如其￥3,3,4、￥4、￥5,6、￥9高变区。
[0009] 在一个实施方案中，通过过滤(filtering)、质量修剪(quality-trimming)、去重 复(de-replicating)和去除PCR引物来获得原始序列读出，从而获得合格序列。
[0010] 在一个实施方案中，通过流程对OTU进行划分，所述流程选自VAMPS、USEARCH(例如 v4、v5、v6、v7 和 v8,如 ￥8.1.1861)、01^^(例如￥1.0^1.1、￥1.8和￥1.9，如￥1.9.1)和]\?)1'冊1? (例如 ￥1.0、￥1.1、￥1.8和￥1.9，如￥1.29.0)。
[0011] 可使用例如Il lumina?测序仪通过Illumina测序法来测定DNA序列，并且从样品中 分离总核酸，然后测序。
[0012] 本发明还提供了用于从环境样品中分离微生物的方法，其中所述微生物包含系统 发生信息基因，所述方法包括:如上所述确定环境样品中的0TU;选择具有其唯一系统发生 信息基因序列的OTU作为待分离微生物;培养样品中的微生物;测定每种经培养微生物的系 统发生信息基因的DNA序列；以及分离其系统发生信息基因的序列与待分离微生物的系统 发生信息基因序列同源的微生物。优选地，使用常规的微生物学、生理学或生物化学参数验 证所分离的微生物。通常，其系统发生信息基因的序列与待分离微生物的系统发生信息基 因序列具有99%或甚至95%或甚至更小同一性的分离物是满意的且被分离出来。
【附图说明】
[0013] 图1示出了通过四种方法获得的"合格序列"根据其与最接近模拟参照的同一性水 平的分布。在质量过滤后修剪对应于PCR引物的序列区域。S + BH + P:与错误校正 (BayesHammer)组合的质量修剪(Sickle)，然后进行读出重叠(PANDAseq)。
[0014]图2示出了"合格序列"的错误类型分布。根据单因素 ANOVA检验，标在上方的字母 不同的数据在p〈0.01下是显著不同的。插图示出了被错误鉴定(包括替换和缺失）的各碱 基的比例之和。
[0015] 图3示出了模拟数据中唯一序列的分布。（A)唯一序列与最接近模拟参照的同一 性。具有>3%错误的唯一序列的相对丰度都〈0.05% JB)唯一序列基于其相对丰度的分布。 大多数唯一序列具有低丰度。
[0016] 图4示出了通过不同方法获得的OTU获取准确度和效率。（A)通过每种方法获取的 总OTU数目。（B)由每种方法所获取的OTU根据其与最接近模拟参照的同一性水平的分布。对 通过Qiime获得的OTU进行另外的嵌合体检测方法(ChimeraSlay er和Uchime )，而Usearch和 Mothur具有嵌入式嵌合体检测。
[0017] 图5不出 了通过（A)Usearch、（B)Qiime、（C)Qiime+ChimeraSlayer、（D)Qiime+ Uchime和(E)Mothur划分为三种类型OTU的唯一序列的散点图。不同的形状各自示出归属于 与模拟参照相同的OTU (完美OTU，圆圈）、与参照具有彡97 %同一性的OTU (良好OTU，方块)或 者具有>3%错误的OTU(虚假0TU，三角形）的唯一序列。
[0018] 图6示出了低丰度"坏序列"对OTU划分的不同影响可通过仅考虑高丰度序列来消 除。圆点代表属于单一质粒的所有唯一序列。中空的不同形状(方块、圆和三角形)代表通过 不同流程获取的OTU的中心。
[0019] 图7示出了模拟数据中用不同的相对丰度阈值通过我们提出的方法获得的0TU。 (A-C) OTU数目和(D-F)重新划归至OTU中心（centro i d)的序列比例达到稳定。在OTU划分的 初始步骤中，仅使用超过阈值的唯一序列。在划分之后，所有的"合格序列"以97%的相似性 重新划归到0TU。
[0020] 图8示出了在(A)PWS、（B)Ultra、（C)Water、（D)River数据中，唯一序列基于其相对 丰度的分布。大多数唯一序列具有低丰度。
[0021] 图9示出了在(A)PWS、（B)Ultra、（C)Water、（D)River数据集中，用不同的相对丰度 阈值通过我们提出的方法获得的OTU。在OTU划分的初始步骤中，仅使用超过阈值的唯一序 列。分割25%的总"合格序列"的阈值适用于所有的真实数据集。不同的流程在这些阈值下 获得接近的OTU数目。
[0022] 图10示出了保留在(A)PWS、（B)Ultra、（C)Water、（D)River数据集中的"合格序列" 的比例。在OTU划分之后，所有的"合格序列"以97%的相似性重新划归到0TU。在分割25%的 总"合格序列"的阈值下，>90%的总"合格序列"可在之后被重新划归回来。
[0023] 图11示出了在(A)PWS、（B)Ultra、（C)Water、（D)River数据中的信噪比与相对丰 度。信噪比随着相对丰度的降低而快速降低并且在更低丰度水平下达到稳定。垂直虚线示 出了 OTU划分的丰度阈值。
[0024]图12示出了在PWS数据中，我们提出的方法得到更少的OTU，但alpha多样性相当。 (A)通过每种流程获取的总OTU数目。（B)各样品的OTU数目。（C)Chaol指数。（D)Shannon指 数。（E)Simpson指数。默认工作流程:执行各流程的默认参数。我们提出的方法:唯一序列的 相对丰度阈值设定为〇 · 00058%。*邙〈0 · 01，*#p〈0 · 001，##p〈0 · 0001 (Kruskal-Wallis检 验，η = 108)。在计算alpha多样性之前，将测序量稀释到4,000个读出/样品。
[0025]图13示出了在PWS数据中，我们提出的方法用不同流程得到了更一致的beta多样 性。通过比较用以下流程获得的每对分析流程之间的Beta多样性距离矩阵获得Mantel统计 量：（A)默认流程，（B)并入不同流程的我们提出的方法。在计算beta多样性距离矩阵之前， 将测序数目稀释到4，000个读出/样品。
[0026]图14示出了构建模拟群落所使用的多个克隆的序列。
[0027]发明详述
[0028] 我们的研究表明，三种常用的OTU划分流程Qiime、Mothur和Usearch提供了模拟数 据中不同的OTU数目和准确度。真实数据集中也出现了这种不同，导致了显著不一致的 alpha和beta多样性信息。
[0029]根据模拟数据的调查，发现这是因为通过常用的质量过滤(quality filtration) 方法不能去除测序错误。这些错误主要分布在具有较低丰度的唯一序列中。虽然这些"坏序 列"的总体丰度较低（总"合格序列"的约5% )，但是在OTU划分中包括这些"坏序列"不仅因 将这些"坏序列"中的许多视为中心而使虚假OTU数目增加，而且因将一些高质量读出分配 到虚假OTU而扭曲了 "真实0TU"的丰度分布。
[0030]其余的错误表明，每个碱基的质量得分不足以指示真实的错误率8，12,28。我们的 研究显示，在根据序列的质量得分"合格序列"中仅允许每条序列有0.5个错误的条件下，存 在具有>3%错误的"坏序列"。嵌合体检测方法同样被广泛实施以克服PCR引入的错误。但 是，无论在OTU划分之前进行嵌合体唯一序列过滤还是在之后对嵌合体OTU进行过滤，都不 能消除我们的研究中的虚假0TU。这些结果表明，目前的质量控制方法的效率不足以去除错 误>3 %的所有序列。
[0031 ]质量过滤同样不改变"合格序列"的错误类型。替换是主要的错误来源，而不是插 入和缺失。C碱基比G碱基明显更可能产生替换错误，并且还出现了替换的趋势:A->G、C->A、 G->T和T->C，与原始读出中所报道的相似12。
[0032]已证实，丰度最高的序列被认为是生物学上的真实序列n，并且其被由较低丰度序 列（大多数为单拷贝序列（singleton))构成的"错误云（error cloud)"包围7'29。因此， Nicholas等将较低丰度OTU过滤掉n，但是这并没有改善OTU的划分。去除较低丰度OTU还意 味着不管所述OTU各自的准确度如何都删除了属于该OTU的所有序列。Chen等不管其准确度 如何都丢弃了 454测序数据中所有的较低丰度序列5,但是根据我们的结果，在Miseq数据中 仅很少一部分的较低丰度序列实际上是"坏序列"。而Edgar在通过Usearch进行OTU划分时 没有考虑单拷贝序列以防止其变为OTU的中心，然后将其重新划归到OTU以实现更好的覆盖 范围 7。我们的方法扩大了在OTU划分期间暂不考虑的唯一序列的范围，原因是单拷贝序列 不是"坏序列"的唯一来源。
[0033] 根据模拟数据，在较低丰度区域中观察到所有的"坏序列"。虽然在真实数据中并 不知道"坏序列"的实际分布，但是我们可以设定阈值以确定避开OTU划分的低丰度唯一序 列。为了找到普遍适用的阈值，我们调查了四个真实数据集中的唯一序列分布和OTU划分结 果。虽然它们用不同的测序方案对16s rRNA基因的不同目的区域进行了测序，但是宿主相 关或自由生存的微生物群落的所有真实数据集始终包含大部分的低丰度唯一序列。当这些 唯一序列不参与初始OTU划分时，OTU数目大大减小并且到达稳定期，在此期间不同流程提 供相似的结果。当低丰度唯一序列占总"合格序列"的25%时，相对丰度阈值对所有的数据 集是合适的。而且，根据自举再抽样(bootstrap resampling)，证明这些阈值水平足以保持 所有的可靠唯一序列。相比之下，较低丰度序列的信噪比表明其丰度确实有很大偏差，并且 不应用于进一步分析。
[0034] 之后将"合格序列"重新划归到预先定义的OTU是另一个重要的过程步骤。该过程 将OTU划分分为两部分：（i)选择每个OTU的中心，和（ii)基于参照的OTU分配。虽然25%的序 列在OTU划分的初始步骤期间暂不考虑，但是最终仅丢弃了不能满足97%相似性阈值的序 列（〈10%)。这为选择OTU划分的中心给出了严格的标准，但是仍然允许将高质量的较低丰 度序列分配到对应的0TU。
[0035]我们的方法防止了较低丰度唯一序列中的假象(artefact)成为OTU的中心，使得 由大多数现有方法产生的OTU数目过高估计降低至合理水平。OTU结果在下游分析和实验中 更为可靠且更易重现，从而促进了重要功能细菌的检测、分离和验证。OTU划分方法的选择 不再是问题，原因是整合有我们的方法的OTU划分流程提供了相似的OTU数目，并且产生了 一致的alpha和beta多样性。此外，我们的方法应用起来简单，原因是其既不需要知道每个 错误的确切来源，也不需要对假OTU进行额外的过滤。其还因仅分析一部分高丰度唯一序列 而降低了对计算资源的需求。我们认为，这种准确、简单、快速且容易整合的方法在微生物 研究中具有潜在的用途。
[0036]本发明提供了定义样品中的微生物操作分类单元(OTU)的方法，所述方法包括：1) 获得样品，所述样品包含各自含有系统发生信息基因的微生物;2)使用基于PCR的高通量测 序技术获得样品中微生物的系统发生信息基因的原始序列读出；3)处理原始序列读出以获 得合格序列片段;4)获得每条合格序列的相对丰度值，其中所有合格序列的总相对丰度为 100 % ; 5)根据合格序列各自的相对丰度值对所有合格序列从高到低进行排序，并将合格序 列分成高丰度组和低丰度组，其中所述高丰度组由丰度值高于低丰度组中合格序列的丰度 值并且总体占总丰度约70%至80%的合格序列组成；所述低丰度组由占总丰度约20%至 30 %的剩余的合格序列组成;6)仅使用高丰度组中的合格序列划分样品中的OTU从而获得 暂定0TU;以及7)重新划归低丰度组中的合格序列到暂定0TU，并且只有在合格序列与OTU序 列具有至少90%序列相似性时才将所述合格序列分别分配到合适的暂定0TU，以实现对OTU 的最终定义。
[0037] 在一些实施方案中，步骤4)、5)、6)和/或7)通过处理器进行。
[0038] 在一些实施方案中，在步骤5)中，高丰度组由丰度值高于低丰度组中合格序列的 丰度值并且总体占总丰度约71%至79%、72%至78%、73%至77%、74至76%、74.5%至 75.5%、74.6%至75.4%、74.7%至75.3%、74.8%至75.2%、74.9%至75.1%的合格序列 组成；低丰度组由占总丰度约21%至29%、22%至28%、23%至27%、24至26%、24.5%至 25.5%、24.6%至25.5%、24.7%至25.3%、24.8%至25.2%、24.9%至25.1%的剩余的合 格序列组成。
[0039] 在一些实施方案中，在步骤5)中，高丰度组由丰度值高于低丰度组中合格序列的 丰度值并且总体占总丰度约75%的合格序列组成;低丰度组由占总丰度约25%的剩余的合 格序列组成。
[0040] 在一些实施方案中，在步骤7)中，只有当合格序列与OTU序列具有至少91 %、92%、 93%、94%、95 %、96 %、97 %、98 %或99 %的序列相似性时才将所述合格序列分别分配到合 适的暂定OTU。
[0041] 因此，在一个实施方案中，本发明提供了定义样品中的微生物操作分类单元(OTU) 的方法，所述方法包括：
[0042] -获得样品，所述样品包含各自含有系统发生信息基因的微生物；
[0043] -使用基于PCR的高通量测序技术获得样品中微生物的系统发生信息基因的原始 序列读出；
[0044] -处理原始序列读出以获得合格序列片段；
[0045] -通过处理器获得每条合格序列的相对丰度值，其中所有合格序列的总相对丰度 为 100%;
[0046] -通过处理器根据合格序列各自的相对丰度值对所有合格序列从高到低进行排 序，并将合格序列分成高丰度组和低丰度组，其中所述高丰度组由丰度值高于低丰度组中 合格序列的丰度值并且总体占总丰度约75%的合格序列组成;所述低丰度组由占总丰度约 25%的剩余的合格序列组成；
[0047]-通过处理器仅使用高丰度组中的合格序列划分样品中的OTU从而获得暂定0TU; 以及
[0048]-通过处理器重新划归低丰度组中的合格序列到暂定0TU，并且只有当合格序列与 OTU序列具有至少90%的序列相似性时才将所述合格序列分别分配到合适的暂定0TU，从而 实现对OTU的最终定义。
[0049] 适用于本发明方法的样品可以是环境样品，例如土壤、水;或大气样品；或来自对 象的样品，例如临床样品，尤其是用于研究肠道微生物区的样品，如粪便样品。
[0050] 操作分类单元(OTU)是旨在代表基于核酸序列的系统发生研究中的分类单元或物 种的个体的群。每个OTU代表系统发生信息基因序列的相似序列变体的群，并且每个OTU可 代表取决于序列相似性阈值的物种或属。
[0051] 系统发生信息基因对于基于基因的系统发生研究领域中的普通技术人员是公知 的，并且是可用于划分两种或更多种生物体的系统发生关系的生物体基因组中的基因或区 域。具体地，系统发生信息基因序列包含足够的随机突变，其数目是自两种或更多种生物体 具有共同祖先以来的时间的结果并且与该时间成比例，从而允许说明生物体的系统发生关 系。众所周知，生物体之间的突变(或差异)数目不应太大或太小，其中的任何一个都将妨碍 破解有意义的关系。
[0052] 许多系统发生信息基因在本领域中是已知且被广泛认可的，这主要通过经验确 定。除了基因的突变率之外，对用于特定系统发生研究的系统发生信息基因的选择还取决 于所研究生物体的系统发生关系。明显地，快速进化的基因仅适于确定密切相关的生物体 (其最近仅根据进化树进行分离)的关系;而突变相对缓慢的基因可适于关系更为疏远的分 类群。
[0053]系统发生信息基因序列的实例包括原核生物中的16s rRNA基因或真核生物中的 18s rRNA基因。具体地，可使用16s rRNA基因的多个高变区（例如，￥1、￥2、￥3 - ￥9)或其相邻 区域，以及ITS(内转录间隔区）或甚至整个16s rRNA基因。
[0054] 通常，在基于16s rRNA基因的宏基因组学研究中，OTU群通过16S基因序列变体的 97%同一性阈值来定义，而且建议使用99%同一性以进行物种分离。
[0055]如本公开内容的上下文中使用的，术语"原始序列读出"意指由自动测序机器的检 测器直接生成的核苷酸序列，与其对应的质量得分一起指示每种核苷酸检测的准确度。 [0056]许多基于PCR的高通量测序技术或"下一代"测序技术在本领域中是已知的并且是 市售的，例如商标名为454SEQUENCER?、IonTorrent?、11 lumina?和PacBio?的多种测序机 器。
[0057] 原始序列读出首先通过双末端测序(PET)方案进行组装，其中PCR引物被截掉。通 过过滤、质量修剪、去重复、去除PCR引物来处理经组装的序列，然后基于平均质量得分或通 过其对应质量得分计算的预计错误率对其进行评估以确定其是否具有期望的准确度。满足 某预先定义标准的经组装序列将被认为是"合格序列"，然后将其归入非冗余唯一序列中。
[0058] 合格序列的长度取决于所使用的区域或相邻区域。通常，合格序列应长到足以提 供有意义的序列信息并且允许测定所研究的系统发生信息基因序列。
[0059]然后计算每条唯一序列的"相对丰度值"，其是唯一序列的丰度除以所有序列的总 丰度。唯一序列阈值的相对丰度根据数据集通常在〇. 0005 %至0.01 %中变化。
[0060] -旦确定了相对丰度，就使用计算机根据唯一序列的丰度值对其进行排序。然后 将序列分成两个组:高丰度组和低丰度组，其中所述高丰度组由具有较高丰度值的合格序 列组成，其总和等于总丰度的约75%;低丰度组由具有较低丰度值的剩余的合格序列组成， 其总和等于总丰度剩余的约25%。虽然在实施例中使用了 75%-25%划分，但是本领域普通 技术人员知道，该分界线可以根据序列的分布进行调整，例如，只要大于90%的总序列可以 被分配到暂定OTU即可。知道不同的划分方法将产生略有不同的OTU数目。
[0061] 在一个实施方案中，本发明使用16s rRNA基因或18s rRNA基因作为系统发生信息 基因，尤其是rRNA基因的一个或更多个可变区。
[0062]在一个实施方案中，本发明结合使用广泛用于OTU划分的流程，例如USEARCH、 QHME 和 MOTHUR。
[0063]本发明还提供了基于如上所测定的OTU系统发生信息基因的序列信息从环境样品 中分离微生物的方法。本发明的"序列指导分离"方法包括：i)在适合样品中所有微生物的 多种条件下进行培养，以获得尽可能多的微生物的纯培养物；以及ii)测定分离物的系统发 生信息基因的DNA序列，并鉴定其相关序列与OTU序列相同或足够相似的分离物。如果待分 离微生物的分类或其他特征是已知的或者基于OTU信息可测定的，则菌落形态或其他常规 的微生物学特性可以用于并且应该用于使需要序列验证的可能分离物的库变窄。
[0064]本发明还提供了定义样品中的微生物操作分类单元(OTU)的方法，所述方法包括：
[0065] 1)获得样品中所包含的微生物的系统发生信息基因的每条合格序列的相对丰度 值，其中所有合格序列的总相对丰度为100% ;
[0066] 2)根据合格序列各自的相对丰度值对所有合格序列从高到低进行排序，并将合格 序列分成高丰度组和低丰度组，其中所述高丰度组由丰度值高于低丰度组中合格序列的丰 度值并且总体占总丰度约70%至80%的合格序列组成;所述低丰度组由占总丰度约20% 至30 %的剩余的合格序列组成；
[0067] 3)仅使用高丰度组中的合格序列划分样品中的OTU从而获得暂定0TU;以及
[0068] 4)重新划归低丰度组中的合格序列到暂定0TU，并且只有当合格序列与OTU序列具 有至少90%的序列相似性时才将所述合格序列分别分配到合适的暂定0TU，从而实现对OTU 的最终定义。
[0069]在一个实施方案中，如下获得合格序列:使用基于PCR的高通量测序技术获得样品 中微生物的系统发生信息基因的原始序列读出，并处理原始序列读出以获得合格序列片 段。
[0070] 本发明还提供了用于鉴定、表征或评估样品中的微生物群落或微生物区的方法， 所述方法包括：
[0071] 1)获得样品中所包含的微生物的系统发生信息基因的每条合格序列的相对丰度 值，其中所有合格序列的总相对丰度为100% ;
[0072] 2)根据合格序列各自的相对丰度值对所有合格序列从高到低进行排序，并将合格 序列分成高丰度组和低丰度组，其中所述高丰度组由丰度值高于低丰度组中合格序列的丰 度值并且总体占总丰度约70 %至80 %的合格序列组成;所述低丰度组由占总丰度约20 %至 30%的剩余的合格序列组成；
[0073] 3)仅使用高丰度组中的合格序列划分样品中的操作分类单元(OTU)从而获得暂定 OTU;以及
[0074] 4)重新划归低丰度组中的合格序列到暂定0TU，并且只有当合格序列与OTU序列具 有至少90%的序列相似性时才将所述合格序列分别分配到合适的暂定0TU，从而实现对OTU 的最终定义。
[0075] 本发明还提供了用于鉴定、表征或评估对象的健康状况的方法，所述方法包括：
[0076] 1)获得来自对象肠道的样品中所包含的微生物的系统发生信息基因的每条合格 序列的相对丰度值，其中所有合格序列的总相对丰度为100% ;
[0077] 2)根据合格序列各自的相对丰度值对所有合格序列从高到低进行排序，并将合格 序列分成高丰度组和低丰度组，其中所述高丰度组由丰度值高于低丰度组中合格序列的 丰度值并且总体占总丰度约70 %至80 %的合格序列组成;所述低丰度组由占总丰度约20 % 至30 %的剩余的合格序列组成；
[0078] 3)仅使用高丰度组中的合格序列划分样品中的操作分类单元(OTU)从而获得暂定 OTU;以及
[0079] 4)重新划归低丰度组中的合格序列到暂定0TU，并且只有当合格序列与OTU序列具 有至少90%的序列相似性时才将所述合格序列分别分配到合适的暂定0TU，从而实现对OTU 的最终定义。
[0080] 本发明还提供了用于定义样品中的微生物操作分类单元(OTU)的软件，所述软件 包括：
[0081] 第一模块，用于获得样品中所包含的微生物的系统发生信息基因的每条合格序列 的相对丰度值，其中所有合格序列的总相对丰度为1 〇〇 % ;
[0082] 第二模块，用于根据合格序列各自的相对丰度值对所有合格序列从高到低进行排 序，并将合格序列分成高丰度组和低丰度组，其中所述高丰度组由丰度值高于低丰度组中 合格序列的丰度值并且总体占总丰度约70 %至80 %的合格序列组成;所述低丰度组由占总 丰度约20 %至30 %的剩余的合格序列组成；
[0083]第三模块，用于仅使用高丰度组中的合格序列划分样品中的OTU从而获得暂定 OTU;以及
[0084] 第四模块，用于重新划归低丰度组中的合格序列到暂定0TU，并且只有当合格序列 与OTU序列具有至少90%的序列相似性时才将所述合格序列分别分配到合适的暂定0TU，从 而实现对OTU的最终定义。
[0085] 本发明还提供了用于定义样品中的微生物操作分类单元(OTU)的系统，所述系统 包括：
[0086] 第一装置，用于获得样品中所包含的微生物的系统发生信息基因的每条合格序列 的相对丰度值，其中所有合格序列的总相对丰度为1 〇〇 % ;
[0087] 第二装置，用于根据合格序列各自的相对丰度值对所有合格序列从高到低进行排 序，并将合格序列分成高丰度组和低丰度组，其中所述高丰度组由丰度值高于低丰度组中 合格序列的丰度值并且总体占总丰度约70 %至80 %的合格序列组成;所述低丰度组由占总 丰度约20 %至30 %的剩余的合格序列组成；
[0088]第三装置，用于仅使用高丰度组中的合格序列划分样品中的OTU从而获得暂定 OTU;以及
[0089] 第四装置，用于重新划归低丰度组中的合格序列到暂定0TU，并且只有当合格序列 与OTU序列具有至少90%的序列相似性时才将所述合格序列分别分配到合适的暂定0TU，从 而实现对OTU的最终定义。
[0090] 通过以下实施例说明本发明，所述实施例不旨在以任何方式限制。通篇所使用的 范围用作描述在该范围内的每个值的速记表达。在该范围内的任何值都可以选择作为范围 的终点。应理解，当描述配方时，如在本领域中常见的，其可以以其成分进行描述，尽管在实 际配方中这些成分在其制备、储存和使用时可彼此反应，并且这样的产物旨在包括在所述 配方中。此外，本申请引用的所有参考文献通过引用整体并入本申请。 实施例
[0091] 材料和方法
[0092] 模拟群落的构建
[0093] 选择总计22种16S rRNA基因克隆来构建7个模拟群落。表1示出了克隆的克隆ID、 分类和序列信息(结合序列表和图14)。
[0094] 每个模拟群落具有不同的克隆组成（表2)。这些克隆的V3-V4高变区彼此具有彡 97%的相似性以避免OTU的混乱估计。每个群落在同一测序批次中具有3次重复。4个群落进 行额外2个批次的测序。
[0095] 获得真实数据集
[0096] PWS数据:获得由诊断患有Prader-Wi Ii综合征或单纯性肥胖的儿童在膳食干预期 间收集的110例人粪便样品4。通过2*300bp双末端测序由同一Illumina Miseq机器对V3-V4 尚变区和t旲拟样品进彳丁测序。
[0097] Ultra数据：已公开的数据集包括来自宿主相关环境或自由生存环境的微生物群 落，以150bp单末端对V4区进行测序3Q。
[0098] Water数据：已公开的数据集收集自荷兰的饮用水系统，以2*200bp的读取长度对 V4区进行扫描31。
[0099] River数据：已公开的数据集包括沿多瑙河中游的水样品，通过2*200bp测序对V3- V4区进行测序32。
[0100]表1:每个克隆ID的序列和分类信息

[0106] http://res.illumina.com/documents/products/appnotes/16s-metagenomic-library-prep-guide .pdf 所述，通过 11 lumina Miseq 对 16S rRNA 基因的高变区扩增子 V3-V4进行测序，其中具有以下修改。使用Platinum Pfx DNA聚合酶（C11708021，Invitrogen, USA)进行两步扩增。将Amplicon PCR(16S rRNA V3-V4区扩增）的PCR循环减少至21以减小 PCR偏差。根据方案进行指数PCR(Index PCR)和PCR产物纯化。所使用的引物对为：5-0_ Bact-0341-b-S-17，5 ' -CCTACGGGNGGCWGCAG-3 ' 和5-0-8&(^-0785-&-八-21，5'-GACTACHVGGGTATCTAATCC-3 '1Q'33。
[0107]质量过滤
[0108] 质量过滤使用Usearch7、Mothur13、Fastq-join 14(在Qiime9中实施），以及最近说明 的流程12来进行，所述最近说明的流程包括质量修剪(Sickle 15)、错误校正(BayesHammer16) 和读出重叠(PANDAseq17)(缩写为S+BH+P)。每个序列对需要彡50bp的重叠，产生彡400bp的 合并序列，并且没有不明确碱基。Usearch进一步过滤掉具有多0.5个预期错误的序列。随后 从"合格序列"中将PCR引物截掉。
[0109] OTU 划分
[0110] 在Usearch7流程中，对"合格序列"进行全长去重复以得到唯一序列，并根据丰度 降低进行排序，丢弃单拷贝序列。然后通过Uparse默认程序获取非嵌合OTU的代表序列。使 用UOHME 2q针对RDP分类训练数据库34(v9)进行进一步的基于参照的嵌合体检测。通过用 Usearch18全局比对算法将"合格序列"划分到剩下的OTU，完成OTU表。
[0111] 按照开发中的Mothur SOP(http://www.mothur.org/wiki/MiSeq_SOP)，对"合格 序列"进行去重复以得到唯一序列，并与SILVA参照数据库 35进行比对。保留在6430位或6430 位之前起始并在23439位或23439位之后结束的序列，并以至多两个差异进行预群集(precluster)。序列按照样品分开，并使用高丰度序列作为参照通过 υ〇ΠΜΕ2() 检查嵌合体。非嵌 合体序列按照RDP分类训练集ν934的Mothur格式化版本进行分类，并进一步过滤掉非细菌序 列。然后用平均邻居算法以>97%的相似性获取0TU。
[0112]在Qiime流程中，使用UCLUST18以>97%的相似性将"合格序列"群集到从头开始的 0TU。使用针对Greengenes核心数据库36的ChimeraSlayer19或者使用针对RDP分类训练集数 据库 34(v9)的UOHME2q，完成嵌合体OTU的额外鉴定。
[0113] 结果
[0114] 用模拟数据进行质量控制过程的评估
[0115] 三个Mi seq测序批次分别完成了每个样品平均15017.4 ± 999.6条(平均值土标准 差）、16247.3 ± 1856.4条和34060.0 ± 3923.9条序列的测序。使用四种质量控制方法来进行 质量过滤，包括Usearch7、Mothur 13、Fastq-jo in14和Sickl e15+BayesHammer16+PANDAseq17 (简 写为S+BH+P)。在各项质量控制和进一步截掉PCR引物之后，通过Usearch全局算法18将留下 的"合格序列"与模拟参照进行比对。整体测序准确度以测序读出与最接近参照的同一性 (图1)表示。
[0116] 如图1所示，质量过滤之后的大多数"合格序列"具有足够的准确度。Usearch、 1〇让1^、？&8七9-」〇丨和3+服+?分别提供了与模拟序列具有100%同一性的68.5±8.9%、 65.9±8.7%、68.9±6.1%和77.1±7.5%的序列。此外，多达94%的"合格序列"与最接近 模拟参照的同一性不小于97%。此外，还存在3.4±1.6%、3.5±1.6%、3.3±1.5%和5.3± 3.3%的"合格序列"具有大于3%的错误，其中的一些与参最接近模拟参照的同一性甚至〈 90% 〇
[0117]另一方面，对于同一模拟群落，虽然四种方法中"合格序列"的错误分布相似，但是 "合格序列"的绝对量在不同测序批次和过滤方法之间有显著不同（表3) Aastq-join和S+ BH+P的稳健性最小，其获得最少的"合格序列"。
[0118]表3.通过使用不同方法的质量过滤的序列数目
[0120] * 通过 PrimerTrim 修剪 PCR 引物
[0121 ] **与错误校正（BayesHammer)组合的质量修剪（Sickle)，然后进行读出重叠 (PANDAseq)
[0122] 模拟数据中"合格序列"的错误类型分布
[0123] 平均地，每条"合格序列"包含1.8±0.8个错误，包括替换(碱基被错误鉴定）、插入 和缺失。当考虑到相似的错误分布（图2)时，根据单因素方差分析(ANOVA)检验，替换的发生 率显著(p〈0.01)高于插入和替换。四种核苷酸具有不同的错误率，其中C比G显著地更可能 被错误鉴定(P〈〇. 01，单因素 ANOVA检验），并且A和T碱基为中等。每种类型的核苷酸还显示 出替换的特定趋势，例如，A更可能被替换为G，T更可能被替换为C，C更可能被替换为A，G更 可能被替换为T。
[0124] 模拟数据中唯一序列的分布
[0125] 对通过Usearch或Mothur流程提供的"合格序列"进行去重复以分别得到25564.7 ±6152.6和35219 ·3± 12133.6 条唯一序列。
[0126] 以Usearch流程获得的结果作为实例，具有>3%错误的唯一序列的丰度为总"合格 序列"的〈0.05 % (图3Α)。通常，对于相对丰度〈0.05 %的唯一序列，相对丰度越低，不同唯 一序列的数目越大，形成了L形丰度分布曲线（图3Β)。以0.05%的相对丰度作为阈值，可将 唯一序列分为相对较低丰度组和相对较高丰度组。发现大于90%的唯一序列在较低丰度区 域中，而其仅占总"合格序列"的27.8 ± 7.8 %。
[0127] 此外，通过〇1；[1116抑31&5^1'19("合格序列"的5.8±1.1%)和1](311；[1116 2()("合格序列" 的3.9±1.8%)检测到的嵌合体序列几乎都属于该较低丰度组。
[0128] 模拟数据中UsearctuQiime和Mothur以默认参数进行的OTU划分
[0129] 在模拟参数集下，Usearch对测序错误表现出最佳耐性，并且将"合格序列"分配到 26.3 ±0.6个OTU Wiime和Mothur获取了比22个物种的实际数据多得多的0TU(799.3±74.5 和429.0± 143.0)(图4A)。由于之前的章节已表明嵌合体不能忽视，所以对通过Qiime获得 的OTU进行了两种额外的嵌合体检测方法(ChimeraSlayer和Uchime )，而Usearch和Mothur 在OTU获取期间或之前实施嵌合体过滤。嵌合体过滤的确改善了 Qi ime的OTU获取，OTU数目 从 799.3±74.5降低至85.0±9.5(〇1址6瓜51&76〇和78.7±4.2〇](^址6)。01'1]数目仍然被 过高估计。
[0130]通过将每个OTU的代表序列与模拟参照进行比对来评估OTU划分的准确度（图4B)。 OTU被定义为"完美"（与模拟参照100%相同），"良好"（97%彡同一性〈100% )和"虚假"（同 一性〈97 % )。所有的方法都可以提供与22个模拟参照的同一性>99 %的OTU，显示出与22个 "真实"物种--对应。但是，Usearch、Qi ime加上嵌合体检查(ChimeraSlayer或Uchime )、以 及1〇让1^分别也获得了4.3±0.6(16.4±1.8%)、61.7±10.1(72.2±4.2%)、52.0±6.0 (66.0±4.4%)和381.7±130.4(88.7±2.1%)个假01'1]。
[0131]然后根据唯一序列所分配到的OTU类型（完美、良好和虚假)追踪唯一序列（图5)。 结果显示，所有的OTU划分方法都保留了"坏"的唯一序列（同一性〈97%)，其中的一些成为 虚假OTU的中心。Usearch在OTU划分期间丢弃了单拷贝序列（没有重复的唯一序列）和潜在 的嵌合体序列，因此明显减少了所保留的低同一性唯一序列。但是只丢弃单拷贝序列并不 足够，因为仍有非单拷贝序列的"坏"唯一序列被保留下来，并成为虚假OTU的来源（图5)。此 外，高同一性的唯一序列（同一性彡97%)同样可以被分配到虚假0TU，即使在经过嵌合体检 查之后也是如此（图5B-E)。这表明，这些"坏序列"不仅通过自身而且通过抓取高同一性序 列引入了额外的虚假OTU，产生了 "真实0TU"的扭曲分布。
[0132]还注意到，保留下来的低同一性唯一序列的相对丰度不超过"合格序列"的 0.05%，并且在嵌合体过滤后进一步降低至〈0.01 %。
[0133]通过我们的方法改善的模拟数据的OTU划分
[0134] 通过模拟数据，认识到具有相对较低丰度的唯一序列是"坏序列"和虚假OTU的主 要来源。应理解，在一条序列中出现的错误越多，另一条序列包含完全相同错误的可能性越 小。因此，通过避免所有的低丰度唯一序列参与OTU划分可以简单地消除"坏序列"。属于单 一质粒的大多数唯一序列群集在一起（图6)，而几条具有相对较低丰度的唯一序列散开分 布，与来源质粒的相似性〈97 %。获取这些序列中的一些作为虚假OTU的中心，使得OTU数目 远远大于一个物种的实际数目。然而，如果在OTU划分期间仅考虑高丰度唯一序列，所有流 程都将得到真实0TU。
[0135] 相对丰度值确定:将序列归并到非冗余唯一序列中，唯一序列的丰度是与原始数 据中该唯一序列完全相同的重复序列的数目。唯一序列的相对丰度是该唯一序列的丰度除 以总丰度。
[0136] 因此，我们提出了三步法来改进现有的分析流程：（i)设定唯一序列的相对丰度阈 值，（ii)仅将超过该阈值的较高丰度唯一序列输入到初始的OTU划分步骤中，（iii)如果较 低丰度唯一序列符合97%的相似性阈值则将其重新划归到0TU。
[0137] 我们设定了一系列的相对丰度阈值来测试我们的方法（图7A-C)。在模拟数据中搁 置"合格序列"的相对丰度〈0.01 %的唯一序列对于Usearch和Mothur已经足够，并且其对应 于每个模拟参照都一致地获得了 22个0TU。当将相对丰度阈值提高到0.05%或更高时，大多 数的结果都获得了与22个"真实"物种相同或非常接近的22个OTU，Qiime+ChimeraSlayer除 外，其仅产生了21个OTU 13OTU数目在进入在阈值高于1 %之前保持稳定，在该水平下真实物 种也不参与OTU划分。
[0138] 当丰度阈值没有超过1%时，由低丰度唯一序列中捕获的最多25%至38%的总"合 格序列"不参与OTU划分的初始步骤，但是之后通过将其重现划分回预先定义的OTU对其进 行了重新考虑。在Qi ime、Mothur和Usearch的结果中进行重新划归之后，至少93.9 %的"合 格序列"被最终保留下来。对经Qiime获取的OTU结果进行额外的嵌合体过滤明显影响重新 划归比例，产生了大幅度的波动（图7D-F)。
[0139] 丰度阈值的确定
[0140] 使用4个已公开的真实数据集来进一步评估我们的方法并找出普遍适用的阈值。 虽然不知道真实数据集中的实际准确度信息，但是在所有的四组数据中在对"合格序列"进 行去重复之后都存在相似的唯一序列L形分布（图8)。每组数据都包含大量的低丰度唯一序 列。将我们的方法并入不同的流程中并改变相似性阈值，对于每个数据集获得了一系列的 OTU划分结果（图9)。所有的结果都显示出，OTU数目在开始时就显著下降，并且当更多的序 列从OTU划分的第一步中搁置时，OTU数目的下降趋势保持缓慢。执行不同算法的不同流程 显示出不同的行为，但是在稳定期其都获得相似的OTU数目。在真实数据集中，这些稳定期 在约0.0005%的相对丰度水平下开始出现（图9)，这比模拟数据集中约0.05%的水平低得 多（图7)。与仅由22个物种构建群落的模拟数据相比，真实数据集要复杂得多。结果表明，不 能选择相对丰度作为普遍适用的阈值，原因是其可能是数据特异性的。然而，如果考虑所搁 置的总"合格序列"的比例，则25%水平都出现在每个数据集的平台期中。在该水平下，至少 90%的"合格序列"可以重新划归到预先定义的OTU(图7，图10)。这意味着，避免最低25%的 总"合格序列"参与OTU划分被认为是潜在的阈值。
[0141]然而，存在这样的担忧:可能丢失具有较低丰度的真实0TU。因此，我们应用自举程 序来估计唯一序列的不确定性水平。对于每组数据，根据唯一序列的原始分布，进行10,000 次有放回的自举重采样。计算每条唯一序列估计的标准误差以及对应的信噪比（丰度/估计 的标准误差）。信噪比随着相对丰度快速降低并在较低丰度下达到稳定（图11)。表4列出了 在四个真实数据集中搁置25%"合格序列"的相对丰度阈值及对应的信噪比值。OTU划分中 包括的唯一序列的最低丰度已小于0.0006%并且信噪比不大于4,这表明所有可靠的唯一 序列都参与了该划分。此外，所有具有非常低丰度的唯一序列仍然可以重新参与OTU的定 量，只要其与预先定义的OTU的相似性为至少97%即可，从而使丢失稀少且真实的OTU的风 险减至最小。
[0142]表4.四个真实数据集中OTU划分中包括的唯一序列的丰度阈值及其对应的信噪比 值
[0144] *平均值土标准偏差
[0145]我们的方法在真实数据集中更一致的alpha和beta多样性
[0146]以PWS数据为例，通过对7,798条唯一序列而非278,160条唯一序列进行OTU划分， 我们的方法大大节约了计算资源和计算时间。其还使该真实数据集的总OTU数目显著降低， 分别为从430降到272 (Usearch )，从7,979到493 (Qiime )，从1 ,671 到302 (Qiime + 〇1：[1116瓜31&5^1')，从1，621到327((>);[;[1116+1](311；[1116)，从4,419到328(]\1〇1:11111')(图12八）。
[0147] 为了比较alpha多样性，计算了每个样品的OTU数目、Chao I21、Shannon22和 Simpson23指数（图12为PWS数据）。前两个指数直接反映样品的丰富性，而后两个指数反映整 体的多样性信息。由于总OTU数目的巨大差异，默认流程与我们的方法之间在对每个样品的 OTU数目（图12B)和Chaol指数（图12C)估计方面存在显著差异(Kruskal-Wallis检验24)。这 两个指数在使用默认参数的OTU流程之间也不同（图12B-C)。但是默认流程与我们的方法之 间在Shannon(图12D)和Simpson(图12E)指数方面没有显著差异。Qiime使用默认参数时的 结果总是显著高于其他结果，这主要是因为其OTU数目大得不合理。
[0148] 为了检验OTU划分的这些差异可如何影响生物学解释，测量了四种beta多样性距 离矩阵，包括欧氏(EU)距离、Bray-Curtis(BC)距离 25、加权归一化Unifrac(WU)距离和未加 权Unifrac(UU)距离26。通过Mantel检验 27比较基于通过不同OTU划分获得的OTU表格计算的 不同矩阵，其相似性用Mantel r统计量来表示（图9为PWS数据，图9-11为其他数据集）。观察 到，在PWS数据中，UU在默认流程中提供了不同的beta多样性估计（图13A，ManteI r = 0.73 ~0.86)，但是使用我们的方法彼此变得更为相关（图13B,Mantel r = 0.85~0.91 KEU和BC 距离在默认流程之间的相关性非常高(Mantel r>0.95)，并且在使用我们的方法的所有方 法之间进一步增加 （Mantel r>0.99)，Qiime+ChimeraSlayer除外（EU的Mantel r = 0.92~ O. 93，BC的Mantel r = 0.92)。方法之间WU距离的相关性不因将默认流程变为我们的方法而 改变。
[0149] 参考文献
[0150] I.Goodrich,J.K.et al.Conducting a Microbiome Study.Cell 158,250-262 (2014).
[0151] 2.Schloss,P.D.&ffestcott,S.L.Assessing and Improving Methods Used in Operational Taxonomic Unit-Based Approaches for 16S rRNA Gene Sequence Analysis.Appl.Environ.Microbio1.77,3219-3226(2011).
[0152] 3.Fei,N.&Zhao,L.An opportunistic pathogen isolated from the gut of an obese human causes obesity in germfree mice.ISME J 7,880-884(2013).
[0153] 4.ZhangjC.et al.Dietary modulation of gut microbiota contributes to alleviation of both genetic and simple obesity in children.EBioMedicine doi: 10.1016/j.ebiom.2015.07.007
[0154] 5.ChenjW.,ZhangjC.K.,Cheng jY.,ZhangjS.&Zhao,H.A Comparison of Methods for Clustering 16S rRNA Sequences into OTUs.PLoS ONE 8,e70837(2013).
[0155] 6.Bonder,M.J.,AbelnjS.,Zaura,E.&Brandt,B.ff.Comparing clustering and pre-processing in taxonomy analysis.Bioinformatics 28,2891-2897(2012).
[0156] 7. Edgar ,R.C.UPARSE: highly accurate OTU sequences from microbial amplicon reads.Nat.Methods 10,996-998(2013).
[0157] 8·Kozich，J·J·，Westcott，S·L·，Baxter，N·T·，Highlander，S·K·&Schloss， P. D.Development of a Dual-Index Sequencing Strategy and Curation Pipeline for Analyzing Amplicon Sequence Data on the MiSeq 11lumina Sequencing Platform.Appl.Environ.MicrobioI.79,5112-5120(2013).
[0158] 9 . Caporaso , J . G . et al.QIIME allows analysis of high-throughput community sequencing data.Nat Methods 7，335-6(2010)·
[0159] 10.KlindworthjA.et al.Evaluation of general 16S ribosomal RNA gene PCR primers for classical and next-generation sequencing-based diversity studies.Nucleic Acids Res.41,el(2013).
[0160] 11.BokulichjN.A.et aI.Quality-filtering vastly improves diversity estimates from Illumina ampI icon sequencing. Nat Methods(2012).doi:10.1038/ nmeth.2276
[0161] 12.SchirmerjM.et al. Insight into biases and sequencing errors for amplicon sequencing with the Illumina MiSeq platform.Nucleic Acids Res. (2015).doi:10.1093/nar/gkul341
[0162] 13.Schloss,P.D.et al. Introducing mothur: open-source,platform-independent, community-supported software for describing and comparing microbial communities.Appl.Environ.MicrobioI.75,7537-7541(2009).
[0163] 14.Aronesty,E.Comparison of sequencing utility programs . Open Bioinform J 7,1-8(2013).
[0164] 15.Joshi,N.A.&Fass,J.N.Sickle:A sliding-window,adaptive,uality-based trimming tool for FastQ f iles(Version 1.33).(2011).at<https://github. com/ najoshi/sickle>
[0165] 16.Nikolenko,S. I., Korobeynikov,A.I.&AIekseyevjM.A.BayesHammer: Bayesian clustering for error correction in single-cell sequencing. BMC Genomics 14，S7(2013).
[0166] 17.Masella，A·P·，Bartram，A·K·，Truszkowski，J·Μ·，Brown，D·G·&Neufeld， J.D.PANDAseq:paired-end assembler for illumina sequences.BMC Bioinformatics 13,31(2012).
[0167] 18.Edgar,R.C. Search and clustering orders of magnitude faster than BLAST.Bioinformatics 26,2460-2461(2010).
[0168] 19.Haas,B.J.et al.Chimeric 16S rRNA sequence formation and detection in Sanger and 454-pyrosequenced PCR ampI icons.Genome Res.21,494-504(2011).
[0169] 20.Edgar，R.C.，Haas，B. J.，Clemente，J.C.，Quince，C.&Knight，R.UCHIME improves sensitivity and speed of chimera detection.Bioinformatics 27,2194-2200(2011).
[0170] 21. Chao , A.Nonparametric Estimation of the Number of Classes in a Population.Scand.J.Stat.11,265-270(1984).
[0171] 22 . Shannon , C . E . A Mathematical Theory of Communication . Be 11 Syst.Tech.J.27,379-423(1948).
[0172] 23.Simpson,E.H.Measurement of Diversity.Nature 163,688(1949).
[0173] 24.Hollander,M.&ffolfe,D.A.NonparametricStatistical Methods .(Hoboken, NJ:John Wiley&SonSjInc.,1999).
[0174] 25.Bray,J.R. &Curt i s , J . T . An ordination of the upland forest communities of sOTUhern Wisconsin.EcoI.Monogr.27,325-349(1957).
[0175] 26. Lo z upone,C. &Knight, R. UniFrac : a new phylogenetic method for comparing microbial communities.Appl Env.Microbiol 71,8228-35(2005).
[0176] 27.Mantel,N.The detection of disease clustering and a generalized regression approach.Cancer Res 27,209-20(1967).
[0177] 28.DePristojM.A.et al. A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nat. Genet.43,491-498 (2011).
[0178] 29.Edgar,R.C.&FIyvbjergjH.Error filtering，pair assembly，and error correction for next-generation sequencing reads.Bioinformatics btv401(2015).
[0179] 30.Caporaso,J.G.et al.Ultra-high-throughput microbial community analysis on the Illumina HiSeq and MiSeq platforms.ISME J 6,1621-4(2012).
[0180] 31.Roeselers,G.et al. Microbial biogeography of drinking water : patterns in phylogenetic diversity across space and time: Microbial biogeography of drinking water systems.Environ.Microbio1.17,2505-2514(2015).
[0181] 32.SaviojD.et al.Bacterial diversity along a 2600 km river continuum: River bacterioplankton diversity.Environ.MicrobioI.n/a-n/a(2015).doi:10.1111/ 1462-2920.12886
[0182] 33.Bertilsson,S.Transitions in bacterial communities along the 2000 km salinity gradient of the Baltic Sea.Isme J.5,1571-1579(2011).
[0183] 34.Cole,J.R.et al.Ribosomal Database Project:data and tools for high throughput rRNA analysis.Nucleic Acids Res.42,D633-642(2014).
[0184] 35. Quast，C .et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools.Nucleic Acids Res.gksl219(2012).
[0185] 36. DeSantis ,T.Z.et al·Greengenes，a chimera-checked 16S rRNA gene database and workbench compatible with ARB.AppI.Environ.Microbiol.72,5069-5072(2006).
【主权项】
1. 定义样品中的微生物操作分类单元(οτυ)的方法，所述方法包括： 1) 获得样品，所述样品包含各自含有系统发生信息基因的微生物； 2) 使用基于PCR的高通量测序技术获得所述样品中所述微生物的所述系统发生信息基 因的原始序列读出； 3) 处理所述原始序列读出以获得合格序列片段； 4) 获得每条所述合格序列的相对丰度值，其中所有合格序列的总相对丰度为100%; 5) 根据所述合格序列各自的相对丰度值对所有合格序列从高到低进行排序，并将所述 合格序列分成高丰度组和低丰度组，其中所述高丰度组由丰度值高于所述低丰度组中合格 序列的丰度值并且总体占总丰度约70%至80%的合格序列组成;所述低丰度组由占总丰度 约20 %至30 %的剩余的合格序列组成； 6) 仅使用所述高丰度组中的合格序列划分所述样品中的0TU从而获得暂定0TU;以及 7) 重新划归所述低丰度组中的合格序列至所述暂定0TU，并且只有当所述合格序列与 0TU序列具有至少90%的序列相似性时才将所述合格序列分别分配到合适的暂定0TU，从而 实现对0TU的最终定义。2. 根据权利要求1所述的方法，其中所述系统发生信息基因选自16s rRNA基因或18s rRNA基因。3. 根据权利要求2所述的方法，其中所述系统发生信息基因是16s rRNA基因。4. 根据权利要求3所述的方法，其中所述系统发生信息基因是16s rRNA基因的一个或 更多个可变区。5. 根据权利要求4所述的方法，其中16s rRNA基因的所述一个或更多个可变区选自V3、 ￥3,4、￥4、￥5,6、￥9高变区。6. 根据权利要求1所述的方法，其中通过过滤、质量修剪、去重复和去除PCR引物来获得 所述原始序列读出，从而获得合格序列。7. 根据权利要求1所述的方法，其中通过公开的流程进行0TU划分，所述公开的流程选 自 USEARCH、Q ΠΜΕ和 M0THUR。8. 根据权利要求1所述的方法，其中使用11 lumina?测序仪通过11 lumina测序法测定 DNA序列。9. 根据权利要求1所述的方法，其中从所述样品中分离总核酸，然后测序。10. 根据权利要求1所述的方法，其中步骤4)、5)、6)和/或7)通过处理器进行。11. 根据权利要求1所述的方法，其中在步骤5)中，所述高丰度组由丰度值高于所述低 丰度组中合格序列的丰度值并且总体占总丰度约71 %至79%、72%至78%、73%至77%、74 至76%、74.5%至75.5%、74.6%至75.4%、74.7% 至75.3%、74.8%至75.2%、74.9%至 75.1 %的合格序列组成；所述低丰度组由占总丰度约21 %至29 %、22%至28%、23 %至 27%、24至26%、24.5%至25.5%、24.6%至25.5%、24.7%至25.3%、24.8%至25.2%、 24.9 %至25.1 %的剩余的合格序列组成。12. 根据权利要求1所述的方法，其中在步骤5)中，所述高丰度组由丰度值高于所述低 丰度组中合格序列的丰度值并且总体占总丰度约75%的合格序列组成;所述低丰度组由占 总丰度约25 %的剩余的合格序列组成。13. 根据权利要求1所述的方法，其中在步骤7)中，只有当所述合格序列与0TU序列具有 至少91%、92%、93%、94%、95%、96%、97%、98%或99%的序列相似性时才将所述合格序 列分配到合适的暂定OTU。14. 用于从环境样品中分离微生物的方法，其中所述微生物包含系统发生信息基因，所 述方法包括： 根据权利要求1确定所述环境样品中的OTU; 选择具有其唯一系统发生信息基因序列的OTU作为待分离微生物； 培养所述样品中的微生物； 测定每种经培养微生物的系统发生信息基因的DNA序列；以及 分离其系统发生信息基因的序列与所述待分离微生物的系统发生信息基因序列同源 的微生物。15. 根据权利要求14所述的方法，其中所述微生物是细菌。16. 根据权利要求14所述的方法，其中分离其系统发生信息基因的序列与所述待分离 微生物的系统发生信息基因序列具有至少95%同一性的微生物。17. 根据权利要求16所述的方法，其中分离其系统发生信息基因的序列与所述待分离 微生物的系统发生信息基因序列具有至少99%同一性的微生物。18. 根据权利要求17所述的方法，其中分离其系统发生信息基因的序列与所述待分离 微生物的系统发生信息基因序列相同的微生物。19. 用于定义样品中的微生物操作分类单元(OTU)的方法，所述方法包括： 1) 获得所述样品中所包含的微生物的系统发生信息基因的每条合格序列的相对丰度 值，其中所有合格序列的总相对丰度为100% ; 2) 根据所述合格序列各自的相对丰度值对所有合格序列从高到低进行排序，并将所述 合格序列分成高丰度组和低丰度组，其中所述高丰度组由丰度值高于所述低丰度组中合格 序列的丰度值并且总体占总丰度约70%至80%的合格序列组成;所述低丰度组由占总丰度 约20 %至30 %的剩余的合格序列组成； 3) 仅使用所述高丰度组中的合格序列划分所述样品中的OTU从而获得暂定OTU;以及 4) 重新划归所述低丰度组中的合格序列至所述暂定OTU，并且只有当所述合格序列与 OTU序列具有至少90%的序列相似性时才将所述合格序列分别分配到合适的暂定OTU，从而 实现对OTU的最终定义。20. 根据权利要求19所述的方法，其中如下获得所述合格序列：使用基于PCR的高通量 测序技术获得所述样品中所述微生物的所述系统发生信息基因的原始序列读出，并处理所 述原始序列读出以获得合格序列片段。21. 用于鉴定、表征或评估样品中的微生物群落或微生物区的方法，所述方法包括： 1) 获得所述样品中所包含的微生物的系统发生信息基因的每条合格序列的相对丰度 值，其中所有合格序列的总相对丰度为100% ; 2) 根据所述合格序列各自的相对丰度值对所有合格序列从高到低进行排序，并将所述 合格序列分成高丰度组和低丰度组，其中所述高丰度组由丰度值高于所述低丰度组中合格 序列的丰度值并且总体占总丰度约70%至80%的合格序列组成;所述低丰度组由占总丰度 约20 %至30 %的剩余的合格序列组成； 3) 仅使用所述高丰度组中的合格序列划分所述样品中的操作分类单元(OTU)从而获得 暂定οτυ;以及 4)重新划归所述低丰度组中的合格序列至所述暂定0TU，并且只有当所述合格序列与 0TU序列具有至少90%的序列相似性时才将所述合格序列分别分配到合适的暂定0TU，从而 实现对0TU的最终定义。22. 用于定义样品中的微生物操作分类单元(0TU)的软件，所述软件包括： 第一模块，用于获得所述样品中所包含的微生物的系统发生信息基因的每条合格序列 的相对丰度值，其中所有合格序列的总相对丰度为1 〇〇 % ; 第二模块，用于根据所述合格序列各自的相对丰度值对所有合格序列从高到低进行排 序，并将所述合格序列分成高丰度组和低丰度组，其中所述高丰度组由丰度值高于所述低 丰度组中合格序列的丰度值并且总体占总丰度约70%至80%的合格序列组成;所述低丰度 组由占总丰度约20%至30%的剩余的合格序列组成； 第三模块，用于仅使用所述高丰度组中的合格序列划分所述样品中的0TU从而获得暂 定0TU;以及 第四模块，用于重新划归所述低丰度组中的合格序列至所述暂定0TU，并且只有当所述 合格序列与0TU序列具有至少90%的序列相似性时才将所述合格序列分别分配到合适的暂 定0TU，从而实现对0TU的最终定义。23. 用于定义样品中的微生物操作分类单元(0TU)的系统，所述系统包括： 第一装置，用于获得所述样品中所包含的微生物的系统发生信息基因的每条合格序列 的相对丰度值，其中所有合格序列的总相对丰度为1 〇〇 % ; 第二装置，用于根据所述合格序列各自的相对丰度值对所有合格序列从高到低进行排 序，并将所述合格序列分成高丰度组和低丰度组，其中所述高丰度组由丰度值高于所述低 丰度组中合格序列的丰度值并且总体占总丰度约70%至80%的合格序列组成;所述低丰度 组由占总丰度约20%至30%的剩余的合格序列组成； 第三装置，用于仅使用所述高丰度组中的合格序列划分所述样品中的0TU从而获得暂 定0TU;以及 第四装置，用于重新划归所述低丰度组中的合格序列至所述暂定0TU，并且只有当所述 合格序列与0TU序列具有至少90%的序列相似性时才将所述合格序列分别分配到合适的暂 定0TU，从而实现对0TU的最终定义。
【文档编号】C12Q1/68GK106055924SQ201610333530
【公开日】2016年10月26日
【申请日】2016年5月19日
【发明人】赵立平, 王景, 张梦晖
【申请人】完美（中国）有限公司