微生物生态学课程作业

使用Illumina MiSeq平台通过多通道16s rRNA测定减少错误的鉴定

DNA测序的费用由于下一代测序技术的进步而降低了。通过Illumina MiSeq平台可以测得大量的序列，而且测序的费用远比454磷酸测序法低。但是也存在一些问题：大量序列的分析以及需要减少3’端错误的序列。这些错误可能导致对微生物群落的错误鉴定。因此准确的分类识别要求我们必须减少这些错误的序列。已经有一些研究通过Illumina MiSeq 构建宏基因组提高通道来增加准确性。在本文中我们会评估通过Misep平台使用250bp配对序列进行微生物鉴定的错误率以及如何通过提高通量来减少错误的鉴定。我们比较了不同条件：改变phix的添加量、改变基因库的大小、通过已知序列选择不同的16sRna基因区域。我们得到的方法可以纠正错误的核苷酸和提高识别精度。总体上来说99.5%的所测序列与相应的序列有着95%的相似性，而且93.6%的所测序列与相应序列有97%以上相似度。这表明MiSeq平台在基因水平分析微生物群体具有高准确性。改进后的方法可以应用于多种环境下通过Misep平台获得大量序列。

关键词:16s RNA基因；Misep；鉴定

1 简介

下一代测序技术已经广泛应用在不同环境微生物的研究中而且提供了这些微生物大量的信息。人们已经开始应用16s rRNA来分析微生物的生物学种类，16S rRNA是细菌和真菌的特征因子。454焦磷酸测序技术由于他可以测得较长序列而在微生物鉴定中得到了广泛的应用。然而,最近,Illumina公司和IonTorrent平台在微生物测定的应用越来越被多。Illumina 公司平台可以测得大量的序列;然而,它的长度相对较短的序列对微生物群落的研究是是一个约束。这个限制可以通过Misep平台延生配对序列到250bp或者300bp得到解决。更长的序列可以对微生物进行更准确的分类。所以Illumina公司平台现在可以应用无微生物群落研究。此外，使用Illumina公司平台得到的大量序列可以减少成本和增加每个样品的研究深度。

将MiSeq平台应用于微生物群落研究需要生物信息学的发展来处理大量的配对序列。一些研究已经对Illumina公司平台获得的基本序列进行了不同方式的分析。在一些研究中，通过大量序列的检测对β多样性进行分析。在另一些研究中,通过质量分数来对获得的序列进行删减。众所周知Illumina公司平台在测序反应中会产生错误，改正这些错误的防范也已经有过相关的报道。此外由于PCR扩增漂移、非特异性引物杂交、较低的退火温度和模板再次退火也会产生诸如二元序列和PCR偏移等错误序列。修正这些错误的序列对于准确的环境微生物群落的分析是必须的。重叠片段中不配对序列的改正方法已经有人提出。但是这种方法不能改正非重叠区域出现的错误序列。使用250bp配对的Misep平台，对虚拟序列样品的重新检测和对这些序列的收集可以减少总体错误率。不幸的是,在删除嵌合体和错误序列后在V4和V5序列得到了98 - 271操作分类单元，而原来只有20个分类操作单元。这样会导致对微生物群落的错误鉴定。因此，改善方法，在最终序列中减少错误的测序要求我们使用Misep平台对微生物群落进行准确的分析。

在16S rRNA中有九个变异度高的区域，这些区域对于不同微生物有着很好的区别能力。然而没有研究证明那个区域对于微生物鉴定来说是最好的，而且对于不同样品的研究会针对不同的区域。比起DNA提取和PCR偏移来说，群落的鉴定更受这些区域的影响。这些区域在使用Misep平台增加序列长度对群落进行研究时同样重要。在以前的研究中，人们对不同区域在100bp和150bp配对序列上用Misep平台进行了研究。V6和V3区域在100bp长

度上进行了研究，V3/V4和V4/V5区域在150bp长度上进行了研究，还有一个研究针对V3/V4、V4和V4/V5区域。

在这篇文章中，我们通过改善对16S RNA的分析方法来减少错误的鉴定。此外，我们会比较不同的变化区域来确定使用Misep平台时的最佳区域，还会确定Phix（醋酸苯基汞）的最佳百分比以及测序库的最佳浓度。改进后的方法对于未来使用Misep平台对微生物群落鉴定将会非常有帮助。

2 材料和方法

（1）制备虚拟序列库和DNA的提取

虚拟数据库由以前从猪粪便中得到的47个已知序列的克隆体组成，这47个克隆体和序列在表1中给出。质粒DNA在47个克隆体中用质粒提取试剂盒提取并配成不同的浓度。提取的DNA浓度用DNA浓度分析盒确定。为了评估改进后的方法和生物分类过程，从猪基因样本中提取粪便基因。猪基因样本由韩国动物研究远提供，我们得到了相关的使用权。

（2）扩增子库的测序

通过硅片测试对不同区域结合后的序列进行分析从而确定扩增子测序，在这个测试中，每个引物得到的可检测序列的长度和比例使用EzTaxon-e数据库确定。序列库和猪粪便中的DNA使用ExTaq聚合酶C1000热循环扩大。扩增条件以前已经确定。放大后的样本使用QIA快速PCR试剂盒纯化并且使用DNA分析盒确定数量。然后用1 μg纯化后的样品根据操作手册组成扩增子库。扩增子库的浓度使用实时定量PCR仪、特定引物、SYBR绿色染料确定。然后用根据手册添加Phix和NaOH。我们通过通过不同浓度的扩增子库和phiX 比例来确定使用Misep平台进行16S rRNA扩增子测序条件。

（3）合并配对序列的测序结果

序列包含不止一个不确定的碱基而且低质量分数的序列会被除去。为了对质量进行过滤，我们比较了每个序列中4种不同长度序列的质量分数。一般来说，使用Misep平台时低质量分数会在3’端出现。所以，我们使用整个序列的平均质量分数来进行过滤。把从Misep 平台得到的序列合并，确定每一个序列中的重叠部分，两端对齐在引物和重叠区域确定不配对序列。引物序列被移除进行其他研究。计算重叠序列的相似性，分析不配对序列和相应质量分数之间的关系。Q分数随着3’端序列的延长而下降，与此同时错误率上升。有两种方

法可以解决这个问题，本文会对这这两种方法进行比较。第一种方法是删减3’端的低质量片段，根据Q分数修正重叠区域的不配对序列。当重叠区域出现不配对序列，优先选择离5’端近的序列，吐过离5’端一样近，那么就选择质量高的那一个。第二中方法是从3’端除去一段特定长度的序列。按以下步骤使用序列库计算测得序列的错误率：1)收索已知的47个序列构成的数据库，找出与合并序列对应的一段序列；2)将找到的序列与合并序列从头对齐；

3)将序列分成十份进行比较；4)计算每一份的错误率。

错误率的公式是：错误率=En/Tmr

其中En表示与相应序列不匹配的的核苷酸总数，Tmr合并序列核苷酸总数

实验将会比较不同浓度序列库相同条件下和同一浓度序列哭不同条件下的错误率，这些错误率将会用来修正重叠序列。

（4）改正错误的核苷酸来进行分类

使用散点图来分析错误序列对微生物鉴定的影响。合并序列和相应序列的相似性用来分析错误序列的影响。降低3’端序列的质量分数与重叠区域的异构序列有关。重叠区域不同的核苷酸和合并序列与相应序列的相似性会放在一张散点图上。不配对的序列会使用Burrows-Wheeler Aligner (BW A)与相应序列比较然后使用Integrative Genome Viewer(IGV) 表示。为了改正错误的核苷酸，相似性在97% 的序列组成一个USEARCH 系统，每个集合都会产生一致的序列。一致序列可以通过将不配对序列与集合中最长的序列从头对齐比较后筛选出来。一直重复这个步骤直到这个集合不在变化。这个集合中还包括前面移除的引物序列。改正的序列用EzTaxon-e 数据库确定，错误用UCHIME系统检测。在集合前和集合后要对错误序列进行两次检测。在第一次就把错误序列移除可以减少后续工作。第二次检测和移除错误序列对于正确鉴定微生物非常重要的。使用97%相似的USEARCH 系统可以检测操作分类单元。这篇文章得到的序列可以在EMBL SRA 数据库中查到。

（5）统计分析

通过统计分析可以确定序列库的不同对错误率的影响和质量分数与错误率的关系。细菌群落的不同用Fisher’s exact test 检验，多种比较采用Bonferroni correction实验检验。错误率和质量分数的关系用散点图研究。错误出现的频率用不配对序列在相应参考序列中出现的百分比表示。

3 结果和讨论

（1）硅片分析来确定引物组

为了确定哪个引物组用来在Misep平台中使用250bp配对扩增子测序，对由两个引物组结合后产生的可检测序列的长度和比例进行分析。结合区域的不同，扩增子的大小也不同。最短的序列是V6 / V7区域结合(193.7±14.0个bp),而最长的序列是V3和V4区域结合(416.8 ±11.2 bp)，然而有一些报道称V3和V4区域结合会产生误差，所以V3和V4区域不是我们所要的区域。在以250bp为基础的Misep平台上，低于250bp的的序列会失去配对优势和鉴定的准确性，因此扩增子低于250bp的也就被排除了。V4与V5区域结合得到的序列有86.05%的序列可以被检测到，这在所有结合后的序列中是最高的。与此相反，V1与V2区域结合后产生的序列可以被检测到奥的是最低的，只有43.87%，这种现象肯能是因为缺乏特定的引物序列，大多输公共区域不包含V1区域的相关序列。V8 / V9区域结合后的检测

率也比较低，因此,这个结合方式也被排除在外。因此我们选择V2和V3(388.8±19.2个bp),V4和V5(372.0±7.4bp)两个区域，这两个区域序列较长、在数据库中有较高覆盖率。

（2）比较phiX的添加百分比和不同的扩增子库浓度

在以前的Misep实验中扩增子库会与50% phiX混合以提高基因多样性。在最近的系统中MiSeq Control software (MCS)可以更好地进行实时数据分析。所以我们可以减少在序列库中添加的Phix比例，在基因多样性比较低的情况下得到更精确的样本数据，就比如16S rRNA 扩增子。在下表中比较了3中不同比例Phix添加量得到的结果。结果表明减少phix的添加量可以提高目标序列的数量和质量分数大于30%的序列。这表明地phix添加量可以产生更多的序列用于扩增子研究。这与先前的报告是一致的，他们添加8.0%phiX获得了9.0×106配对序列,其中80.1%是≥Q30而添加6.2%phiX 得到了10.5×106配对序列,其中74.6%是≥Q30。

在每个实验中我们比较了不同的扩增子库的浓度，每个样品得到的平均序列在表1C中表示。在添加50% phiX的情况下，不同扩增子库浓度得到的序列差不多。但是在添加10%和5% phiX 的情况下，得到的序列随着扩增子库浓度的增加而增加。使用4pm扩增子库可以得到超过6×105 序列，使用8pm扩增子库可以得到超过12×105序列。序列从6pm到8pm 的增加量小于从4pm到6pm的增加量。在以前一个10pm扩增子库得到的序列和本篇中8pm 的到的序列差不多。在以前的研究中表明10pm可以比5pm多得到2×105序列。合适的扩增子库有利于得到高质量分数的序列。因此我们采用8pm扩增子库和10% phiX添加量来进行本实验。

（3）合并配对序列

图2中给出了每个序列不同4中不同质量分数下得到的数据，大多数序列中在使用低质量分数过滤时被去掉了，因为低质量分数的序列大多出现在3’端。所以我们使用整个序列的平均质量分数进行过滤，经过质量过滤后，合并配对序列可以延长扩增子库的长度，预期的扩增子库序列长度是389bp (V2/V3) and 372 bp (V4/V5)，每一个引物序列的5’端被除去来进行分类。重叠序列前后共有80bp，假设目标长度是380bp。

两种不同方法合并后的序列与已知序列进行比较，在3’端减去40bp的序列由于长度就不会产生重叠序列。通过在3’端修剪不同长度序列后得到的序列在表2中给出。在重叠区

域中得到的确定序列随着3’端减去序列的增加而增加，这也表明更多的错误序列在3’端生成。重叠区域中质量分数高的序列可以不经过修剪而直接用于数据分析。当在V2和V3、V4和V5结合区域减去30bp和40bp序列后冲得区域的确定序列会显著增加。这表明在V2和V3结合区域减去30bp的序列、在V4和V5结合区域减去40bp的序列可以出去重叠区域的大部分错误核苷酸。因为V2和V3结合区域长度是389bp所以在3’端减去40bp后会减少重叠区域的确定序列，所以减去40bp不会产生重叠区域。如果不修剪的话比较重叠序列和配对序列得到的可纠正序列就会增加，这表明在3’端删减序列会降低通过合并改正序列的可能。虽然修剪后的合并序列长度和不经修剪的合并序列长度差不多，但是纠正错误的序列可以提高测序的准确性。例如，在V4和V4结合区域3’端减去40bp后的序列(798,511）和不经修剪的序列(798,989)长度差不多，但是在重叠区域纠正之后的序列可以增加序列的可读性。

改正后的序列的准确性通过比较改正后的序列与序列库中对应序列而得到。在重叠区域，改正不同的核苷酸可以降低错误率，这表明在重叠区域改正错误的序列会增加测序的准确性。我们比较每个区域对应序列平均质量分数的错误率（图2A），最高的错误率出现在引物序列后30bp，然而这个区域的平均质量分数超过了35，这表明错误率与质量分数没有关系，每一等分的错误率也与质量没有关系（图2B）。这与先前一个报告中体出的错误率与低质量分数有关是相反的，这种不同可能与测序条件的不同有关，比如引物和序列库不同。然而还有一篇文章高估了移除错误的核苷酸和改正序列后会导致所有区域出现的分类操作单元。21哥分离体混合在一起后，有关出现20个分类操作单元，但是在V4和V5区域出现了超过98个分类操作单元这表明高质量分数的错误序列不经过错误改正步骤也会被检测

出来，这种高估的信息可能在对扩增子库分析是出现错误的鉴定。

(4)模型和错误序列的影响

我们比较使用不同浓度测序库和添加不同比列phix时候产生的错误率进而错误率的产生和序列区域之间的关系。同一次实验中使用6pm序列库和使用8pm序列库所产生的错误率是差不多的，然而在使用6pm序列库和8pm序列库的错误率模型是不一样的。使用6pm 序列库是十份中的第三份错误率是最高的，第九份次之；在不同的实验环境下得到的模型与前面提到的两个模型又是不一样的，这表面错误率不是在某哥特定区域产生，而是在整个序列上随机产生。使用V4和V5区域结合所鉴定的微生物群落和使用V2和V3区域结合鉴定的微生物群落是差不多的，因此我们对V4和V5结合区域进行了后续研究（图3）。

为了研究错误的核苷酸对微生物群落分析的影响，对合并序列与相应序列的相似性与配对序列中重叠序列的相似性在散点图中进行比较，对这两相似性进行研究是因为在以前的报告中重叠区域的不匹配数量对序列的确定又影响。图中的两条线用来确定每个轴的平均值。重叠区域的不匹配序列数量与合并序列的错误序列之间没有关系，重叠区域不匹配序列的平均值是2，这表明大部分序列的重叠区域错误的核苷酸数量小于2，然而合并序列与相应序列的相似性的平均值低于97%，而高于97%的合并序列的重叠区域错误的核苷酸数量却超过了两个，这是因为错误核苷酸是在序列的全部区域随机生成的。合并序列与相应序列相似度低于97%的序列可能是在微生物群落鉴定中的潜在错误信息，即使去掉重叠区域包含两个不匹配核苷酸的序列，错误率仍然会产生。在对生态环境中的微生物群落进行鉴定是没有参考的序列库，所以，改正错误的信息对于准确的微生物鉴定是艰难但是必须的。

（5）纠正错误的核苷酸

集合可以减少错误的序列，在分析的时候用集合已经有过报道。使用集合中一致的序列可以减少随机产生的错误，一致的序列是在集合中对齐端点后筛选出来的（图4），不同序列库条件下减少的集合见表3，一般来说在经过3次集合之后集合的数量就不会再改变(不

同目标地区、图书馆浓度和测序运行)。经过第一次集合之后会得到超过1000个集群，集群的数量进经之后的每个步骤都会减少。V4和V5区域扩增子集群扩增子数量低于V2和V3地区域扩增子集群数量，也低于在另一个测序条件下得到的集合数量。V4和V5区域添加5%phix得到的合并序列圆圆超过了添加10%phix得到的合并序列。这表明增加序列的数量会导致更多的错误；所以，我们有必要确定合适的phix添加量。

通过比较纠正后的配对序列在序列库中对应序列的比列可以评估集合和改正序列的作用。其他两种合并序列的方法也在图4中进行了比较，不经过集合时，合并序列中与对应序列相似度超过97%的占总序列的75%，相似度超过95%的占总序列的93%，这些相似性可以用来确定系统进化树。采用这种方法，可以对那些相似度低于97%的序列进行分析。丢去重叠序列中超过2哥错误核苷酸的序列会导致7%的序列是错误的。而采用集合改正方法之后，,占总数的99.5%读序列与参考序列的相似度大于95%,93.7%的序列与参考序列的相似度大于97%。尽管6.3%的序列不同于相应的参考序列,但是获得了95%的相似度的高度准确的信息。这种相似性已经可以在属间进行分类。因此使用集群和改正方法可以在250bp Misep 平台对微生物群落在属间进行准确的鉴定。

（6）比较使用不同引物组和不同分析方法得到的群落信息

使用两种引物组扩大后得到的扩增子基于他们的组合与序列库在门和属水平上进行比较。实验结果显示使用V4和V5区域结合后的扩增子比V2和V3区域结合后的扩增子更接近与序列库，扩增后的基因比例与原来的不同，这是因为每个基因的扩增效率不一样。然而V4与V5区域结合后产生的序列与初始序列库最接近而且具有很高的分类准确性。V2与V3区域结合后的扩增序列甚至在门水平上就不同。因此，我们认为V4与V5区域是基于250bp使用Misep平台测序的合适区域。

从猪粪便样本中使用不同的分析方法得到的细菌种类在图5中进行了比较，在门水平上，由合并序列确定的种类和由集合改正后的序列确定种类差不多，然而除去重叠区域有两个不匹配核苷酸后的序列确定的细菌种类与前面得到的种类是不同的，这可能是因为去掉重叠区域导致。去掉不匹配序列后得到的序列比合并序列和改正序列短，使用除去双不匹配序列原则可能导致错误的鉴定。集合和改正方法可以分析大部分序列而且可以改正错误的序

列，在合并序列中可以观察到19个门，而在集合和改正的序列中可以观察到11个门；合并序列中可以检测到537个属，而在集合和改正的序列中可以检测到173个属。这表明改进后的方法可以消除对生物多样性的过高估计。我们把我们的结果与先前的报告与序列库进行比较，在以前的报告中有参考序列库外的结果，相反的是，我们得到了参考序列库的结果。这表明本文中改进后的方法可以提高微生物鉴定的准确性。

在这篇文章，我们评估了使用Misep平台测序的不同条件并且提出了改正错误序列准确鉴定微生物的一种方法。使用Misep平台可以产生大量的序列，但是同时也会产生微生物鉴定的一些错误信息，因此为了进行准确的鉴定，我们必须改正错误的序列。V4和V5区域扩增后可以产生最精确的信息，用于硅片测试和序列库分析中。使用Misep平台时，8pm序列库和10%phix添加量可以得到最好的结果。采用集合和改正的分析方法，可以微生物多样性的高估，也可以减少对序列库错误的分析。我们希望其他人在使用Misep平台对微生物群落进行鉴定时候可以使用我们提到的集合和改正的方法。本文中提到的方法对于以后微生物的研究会有很大的帮助。