Nucleic Acids Research│PacBio CCS + DADA2→ASVs, not OTUs
文章核心点
利用PacBio CCS测序模式与新的DADA2鉴定分析方法，从菌群中鉴定出不同菌种的ASVs（ASVs不同于以往的OTUs，它具体到单个的SNP，而不是依据相似度聚类），可以将不同的菌种区分开，准确度与分辨率都非常高。背景与目的
16S rRNA基因片段靶向PCR扩增与高通量测序（扩增子测序）广泛应用于鉴定微生物种群。新的长读长测序技术（PacBio与Nanopore）能够测通16S rRNA基因全长，然而颇高的错误率限制了它们的应用。本研究介绍一种基于PacBio CCS扩增子测序模式与DADA2软件算法，以获得单碱基/近零错误率的全长16S rRNA基因序列。利用此方法对2个Mock communities测序分析，结果获得了预期种群数量的全长16S序列（零错误率），并且成功鉴定出的基因组内序列变异，将E.coli菌株准确地划分成O157:H7与K12亚群分支。作者又用此方法对人类的粪便中的复杂菌群进行了技术重复模拟，获得了若干菌株中16S rRNA基因，成功地鉴定到种水平。此外，此技术方法在其他测序分析领域也将会发挥重要的用途。
材料方法
样品——Zymo mock community；HMP mock community；人类粪便
引物—— 27F: AGRGTTYGATYMTGGCTCAG；1492R: RGYTACCTTGTTACGACTT
测序——PacBio Sequel / CCS；S/P1-C1.2， S/P2-C2/5.0， S/P3-C3/5.0
算法与软件：DADA2 R；uparse；mothur
DADA2 R package分析流程——
图1. 算法流程图
研究结果
1、16S全长扩增子测序准确度评估
本研究利用2个mock菌群评估了全长16S rRNA扩增子测序方法。菌群1——Zymo mock菌群包含8种进化关系较远的细菌菌株，且每种菌株DNA含量均衡；菌群2——HMP mock 菌群包含20种细菌菌株，具有不同的遗传进化相似度，且rRNA基因频率变化超过3个数量级。利用S/P1-C1.2 测序化学试剂进行了单个cell的Sequel测序，Zymo mock菌群产出77 453条CCS reads（Passes＞3；准确度＞99.9%），除去引物与过滤后，还剩余 69 367条reads (Supplementary Table S1)；HMP mock菌群产出78 328条CCS reads，过滤去除引物后，剩余69 963条reads。
本研究利用DADA2 R软件包最新版本1.12.1处理长片段扩增子reads获得精准的ASVs，接近标准的PacBio CCS测序准确度。朴素贝叶斯分类算法与SILVA v128数据库将ASVs注释到“属”以下水平，而在“种”水平利用BLAST比对NCBI数据库。在这两个数据集中，每个ASV注释到对应的一个属与种，符合mock种群中对应的身份。为了实现所期望的检测出多拷贝16S rRNA基因中存在的变异，根据分类学特点将ASVs划分到假定的菌株bins中。
图2. 全长16S rRNA基因 ASVs丰度比较（Zymo mock community）
在Zymo mock种群中，检测到29个ASVs，其中25个与已知的16S rRNA基因序列完全一致（100%覆盖度、100%一致率），3个发酵乳杆菌ASVs与1个金黄色葡萄球菌ASV跟已知的序列仅存在1个碱基的差异。在HMP mock种群中，共检测到51个ASVs，其中48个序列完全匹配，而1个表皮葡萄球菌ASV存在2个碱基的差异，2个拜氏梭菌ASV存在1个碱基的区别，并且ASVs频率变异范围跨越3个数量级，最低达到0.00019。
如果ASVs代表属水平的等位变异，那么同一菌株ASVs根据每个基因的拷贝数目而存在一定的比例。相反，如果ASVs所代表的变异来自不正确的扩增子测序样品，这个比例将不是所预期的。此方法所检测到的ASVs彼此间呈现一定的频率或比值，并与相应的菌株存在对应性。图2展示的是在Zymo mock种群中每个ASV的频率，所有ASV：基因组频率均是整数（例如，1,2,3,.....），最大偏差±0.2。在HMP mock种群中，除了蜡样芽孢杆菌ASV频率受菌株丰度影响，其它的ASVs频率也是整数（图3.）。
图3. 基因组与ASVs丰度（HMP mock community）
2、CCS扩增子序列错误率与降噪后ASVs
这些mock种群可利用的参考数据库是不完整的，因此，作者利用多种方法（全面）评估ASV准确度，其中整合了参考数据库比较、属间等位变异频率模式。在Zymo mock种群中，作者发现此方法检测到的ASVs不存在假阳性，并且利用DADA2降噪后，碱基错误率为零。在HMP mock种群中，作者发现检测到的ASVs不存在假阳性，每个碱基错误率也是零。然而，作者仍然发现了这种可能性，蜡样芽孢杆菌ASVs包含错误，并且错误率为2.6 × 10−6。
接着，作者利用Zymo mock种群研究未经矫正的PacBio CCS扩增子reads的错误率分布情况。排除嵌合体与污染物后，将剩余的reads与参考ASVs序列进行比对，序列间所有的差异以变异类型（替换、插入、缺失）、read的碱基坐标以及对应的质量值来表征与衡量。图4展示的就是错误率分布，200bp为一个可视化代表性窗口。
未矫正的CCS reads总错误率是4.3×10−4/碱基，其中插入错误率为2.2 × 10−4，替换为1.1× 10−4，删除为1.0 × 10−4。替换的错误率分布相对比较均匀一致，高错误率数量较少，且与质量值低相关（图4.）。插入型错误率在reads的不同位置差异较大，可能是由于context-依赖的错误模式，比如同聚物错误，但是高错误率的碱基与低质量值紧密相关。而删除型错误率分布变异系数介于插入与替换之间，因为删除型碱基没有相应的质量值，并且某些位置高频率的缺失与技术没有明显差异。
据之前报道，RSII产出的16S rRNA全长PacBio CCS数据在降噪前后存在较高的错误率，可能也与特定位点重复出现的系统性错误有关。作者利用升级的DADA2 R软件包分析了Wagner et al.的金黄色葡萄球菌测序数据（之前用于评估系统性错误），获得了5个非嵌合型ASVs。金黄色葡萄球菌基因组包含5个rrn操纵子拷贝，并且5个ASVs与已知测序的16S rRNA基因序列完全匹配。基因组内的变异体差异可能被误解为系统性错误，因为利用 short-read组装序列，会将5个rrn操纵子拷贝视为1个。此外，还发现当subread的覆盖度高于默认阈值3的时候，并不能提升ASVs的准确度，暗示默认阈值对于Sequel测序试剂与当前的软件版本是有效的。
图4. PacBio CCS扩增子插入/缺失/替换的错误率与质量值分布（Zymo mock community ）
3、DADA2-ASVs与长读长-OTU鉴定方法准确性比较
作者将DADA2计算方法与在PB扩增子测序数据中先前应用的两款软件算法mothur和uparse进行比较分析。所使用的数据均为Zymo与HMP mock数据集中预处理的fastq数据，分析结果见下表1。
表1：不同计算方法的准确度比较（PB扩增子测序数据）
所有算法表现得都相当不错，反映了PacBio CCS reads具有高的准确度，但是算法间的差异也是非常明显的。DADA2与uparse具有相当好的鉴定特异性，因为鉴定到的每个ASV or OTU真实存在与预期的mock种群中。Mothur算法分析得到16个与13个假阳性OTUs（对应于Zymo与HMP mock种群），大部分假阳性的OTUs以单体形式存在（例如，OTUs仅有1条read），因此，后续分析需要将其删除。这3种方法都在Zymo mock菌群中检测到8种菌株，而在HMP菌群中，DAD2算法检测到17个菌种，另外两种算法仅检测到16个菌种。DADA2可以区分金黄色葡萄球菌与表皮葡萄球菌，然而， mothur与uparse却将这2个菌种划归为同一个OTU。Mothur算法也不能检测到中等丰度(>100 reads)的淋球菌菌株，但是却是唯一能够检测到仅有3条reads的耐辐射奇异球菌的方法。
4、利用16S rRNA等位基因的全长互补序列进行亚种分类
系统分类典型的分类方法是依据个体的16S rRNA序列进行，但是利用细菌菌株中16S等位基因的全长互补序列可以获得比较高的分辨率。在Zymo mock菌群中，鉴定到5个ASVs，丰度比例为3:1:1:1:1；在HMP mock菌群中，鉴定到6个ASVs，丰度比例为2:1:1:1:1:1。作者考虑利用单个ad hoc完全互补序列分类策略将每个ASV的最优BLAST hits进行Nt比对，与基于菌种最大数量的ASVs reads选择最优的hits区分不同亚种。下表2即是基于长读长与短读长以及不同的分类策略进行的分类结果统计。
表2：基于全长与短读长16S rRNA基因扩增子测序的E. coli菌株分类
利用本研究中的方法鉴定到的16S rRNA基因序列（Zymo mock菌群）被归类为肠出血性大肠杆菌O157:H7，并且利用16S rRNA等位基因全长互补序列进行分类可以进一步提升准确度。有14条Hits与5个ASVs中的4个完全匹配，但是没有没有一个Hit与所有的5个ASVs匹配。其中，有12条被注释为血清型O157:H7，1条注释为O157而没有H抗原，还有1条没有血清型注释。如果我们忽略缺失值，Zymo菌株可以被明确地归类为肠出血性大肠杆菌O157:H7。相反，短读长V3V4与V4V5扩增子测序不能区分Escherichia 与Shigella两个属，更不用说区分亚种。
利用本研究中的方法鉴定到的16S rRNA序列（HMP mock菌群）将菌株归类为K12。有32条Hits完全匹配菌群中的6个ASVs，其中，27条注释为K-12, MG1655或源于MG1655菌种。另外5个没有被注释，但是最好的BLAST hits是K-12菌株。利用最大丰度全长序列进行分类不能区分所属的成员，仅注释为E.coli。同样，短读长V3V4与V4V5扩增子测序不能区分Escherichia 与Shigella两个属。
5、PacBio CCS & DADA2在人类粪便样品中的应用价值评估
为了评估新方法在复杂微生物菌群中的应用效果与价值，本研究采集了9份人类粪便样品进行全长扩增子测序，并设置技术重复，另外还有3份样品没有用于进一步分析。每个样品加上barcode，并将12个样品混合后，利用单个Sequel cell进行多通道测序。第一个重复使用S/P2-C2/5.0 Sequel试剂，产出177 691条 CCS reads，准确度达到99.9%，去除引物与过滤数据后，还有146 589 reads，12个样品中过滤掉的reads中位数为12 158。第二个重复使用预释放版S/P3-C3/5.0 Sequel试剂，产出289 644 CCS reads，准确度也均在99.9%，除去引物与数据过滤后，剩余 249 802 reads，12个样品中过滤掉的reads中位数为21 411，几乎是第一种测序试剂的2倍产出。
由于不知道人类粪便样品中微生物菌群的组成，作者考察了技术重复间的一致性，以评估数据的准确度。为了排除任何测序深度的影响，在利用DADA2软件分析之前，作者将每个重复样品筛选出10000条reads，共有7对样品在2次重复数据中数据超过10000条reads。DADA2从这些样品中鉴定出的ASVs在技术重复间是一致的（图5A），以每个样品为基础，两次重复中鉴定到881个ASVs，然而，分别有65与70个ASVs仅在重复1或重复2中被检测到。丰度评估也是高度一致的：样品间丰度Pearson’s相关系数为0.998，然而，仅在1个重复中检测到的ASVs具有比较低的频率(<50 reads, <0.5% frequency; 图 5A)。
创新点总结
百迈客优势
1. PacBio的CCS模式minPasses≥5，自我校正，超高准确性
2.三代全长序列中携带的特征信息更多，可分辨能力更高，可以支持“种”水平的分辨率
3. 三代全长多样性产品类型：16S全长；ITS全长；18S全长 功能基因多样性全长
4. 微生物多样性云平台全面开放，主流程个性化均免费分析
文：文风刺骨
排版：市场部