青岛大学苏晓泉组-微生物组16S扩增子功能校正算法Meta-Apo
Meta-Apo：微生物组16S扩增子测序的功能校正算法
Meta-Apo improves accuracy of 16S-amplicon-based prediction of microbiome function
BMC Genomics [3.594]
原文链接：https://doi.org/10.1186/s12864-020-07307-1
第一作者：Gongchao Jing （荆功超）
通讯作者：Xiaoquan Su （苏晓泉）
主要单位：青岛大学计算机科学与技术学院，中国科学院青岛生物能源与过程研究所单细胞中心
摘要
背景： 16S rRNA基因扩增子测序（以下简称16S扩增子）与鸟枪法宏基因组测序（Shotgun Metagenomic Whole Genome Sequencing；以下简称WGS）相比，前者在实验和计算上的成本低得多，因此目前广泛使用16S扩增子测序来预测微生物组的功能谱。但由于PCR扩增偏好性以及16S rRNA基因-全基因组关联信息的差异，同一微生物组样本基于16S扩增子的功能谱与WGS产生的结果之间存在偏差，从而导致相左的结论。
结果： Meta-Apo可以极大地减少甚至消除这种偏差，使两种方法得出结论更加一致。我们对来自4个身体部位超过5,000例人体微生物组的16S 扩增子样本进行测试发现，Meta-Apo仅使用15个WGS: 16S扩增子配对样本来进行训练，就可以显著降低两种测序之间功能解析的结果。另外，Meta-Apo还可以在WGS和16S扩增子样品之间进行跨平台功能比较，可以极大的改善基于16S扩增子的微生物组诊断。总而言之，利用Meta-Apo，可以让低成本的16S扩增子测序产生与WGS相近的、可靠的、高分辨率的微生物组功能图谱。
结论： 使用Meta-Apo，可以让大规模，以功能主导的微生物测序项目应用成本更低的16S rRNA扩增子策略，获得与WGS相同的准确率。在GitHub上提供了Meta-Apo优化的C ++实现（https://github.com/qibebt-bioinfo/ meta-apo）。它以少数WGS：16S扩增子配对样本的功能谱作为训练集，可以对大量16S扩增子样本校正功能谱。
关键词
微生物组(Microbiome),基因组(Metagenome),扩增子(Amplicon),功能(Function),校正(Calibration)
背景
微生物群落的高通量测序通常采用以下两种策略[3]：16S rRNA基因扩增子（以下简称16S扩增子）测序和鸟枪法全基因组测序（以下简称WGS）。 在扩增子测序中，一般使用系统发育标记基因特定区域（如16S rRNA基因）的引物进行PCR扩增，来揭示微生物群落的物种组成结构，但不能直接得到功能基因的图谱。而WGS虽然可以直接构建功能谱[4]，但其实验和计算有更高的成本和复杂性 [3，5，6]。 因此，出现了使用16S扩增子预测功能谱的计算工具如PICRUSt等[7-10]，由于16S基因在PCR中的扩增偏差[11，12]，以及16S rRNA基因-全基因组关联信息的差异，同一微生物组样本16S扩增子预测的功能谱可能与WGS产生的功能谱有很大差异（图1和图3a）。
为了解决这一难题，我们开发了Meta-Apo（Metagenomic Apochromat）。仅使用少量匹配的WGS：16S扩增子数据对用作训练集（即每一个样本都分别通过WGS和16S扩增子测序），Meta-Apo就可以为大规模的16S扩增子样本生成校正后的功能谱，结果与WGS更加一致（图1）。此外，由于WGS比16S扩增子测序在微生物组疾病检测[13，14]中更加稳定，因此使用Meta-Apo可以让WGS样本和16S扩增子样本进行跨平台的功能比较，极大改善16S的微生物组的诊断。例如，使用校正的16S产生的功能谱，牙龈炎的诊断准确性从65％提高到了95％。
图1. 通过对少量成对的WGS：16S扩增子样本进行训练来校正微生物组扩增子样本的预测功能图谱
结果
16S扩增子和WGS的功能谱：不一致，但同构。
为了评估两种测序策略之间微生物组功能图谱（使用KEGG Orthology [15]（KO）注释）的差异程度，我们首先比较了622例配对的人体微生物组功能谱（数据集1；来自四个身体部位：肠道，皮肤，口腔和生殖道；表1），每个数据集都通过鸟枪WGS和V3-V5区16S rRNA扩增子进行测序。对于WGS，功能谱由HUMAnN2[17] 分析生成。16S则使用PICRUSt 2[8]预测得出。通过比较从两种测序方法得出的功能谱，我们发现配对的WGS：16S扩增子之间差异显著高于WGS的内部差异（即来自同一部位的WGS样本之间的距离）。两种策略之间的差异十分显著，β多样性也表现出非常不同的模式（图3a； PC1 双尾配对 Wilcox秩和检验p < 0.01； PC2双尾配对 Wilcox秩和检验p < 0.01）并导致了一些错误的分类。例如，一些皮肤的16S扩增子的功能谱与口腔的WGS的功能谱被错误的分成一类。
然而，这两种方法得到的WGS与16S扩增子之间距离高度相关（图3b； Pearson相关性R = 0.86，p < 0.01），而且其β多样性之间的总体形状相似（图3a；蒙特卡洛检验p < 0.01）。因此，可以将16S扩增子预测的KO功能注释（
）通过等式(1)链接到WGS的KO功能注释（
）：
()
表1. 本研究中使用的WGS和扩增子数据集
通过线性回归建模减少WGS和16S扩增子数据集之间的功能谱差异
为了减少16S扩增子和WGS数据集之间的功能谱差异，我们开发了Meta-Apo。该软件主要包含两个部分：训练和校正。在训练部分中，Meta-Apo使用线性回归建模利用少量的WGS:16S配对样本来估算等式(1)中的
。在校正部分中，将WGS结果视为“黄金标准”，使用模型
校正16S扩增子样本的预测功能图谱。
为了定量评估Meta-Apo的效果，我们分别从数据集1中随机选择了N = 5、10、15、20、50和100个WGS:16S扩增子配对样本作为训练集，并使用Meta-Apo校正该数据集中其他16S扩增子样本。当使用N = 15个训练对建立模型
时， Meta-Apo校正效果变得稳定，并且在增加更多训练对之后（最多100个；图2b），校正效果也不会明显增加。在校正后（即N = 15个训练对），配对的WGS：16S扩增子距离（0.121±0.055）显著低于WGS样本的组内距离（0.136±0.056）。 经主坐标分析（PCoA）证实，Meta-Apo消除了两种测序策略产生的样本之间的总体功能分布差异（图3c；PC1双尾配对 Wilcox秩和检验p = 0.30，PC2双尾配对 Wilcox秩和检验p = 0.29；图3d）。此外，对KEGG BRITE所有代谢通路层级上的信息做进一步比较，校正后的16S扩增子与WGS得到了更一致的组成（图4）。而且Meta-Apo对于来自数据集2（表1）的V1-V3区16S rRNA序列也同样适用。
图2. Meta-Apo显着减少了数据集1中WGS和16S扩增子配对样本之间的功能谱的距离。a. WGS：16S扩增子配对样本之间的BrayCurtis距离（未校正，橙色条）高于WGS体内位点距离（来自同一部位的WGS样本之间的距离，蓝色条）。 b. 仅使用15个训练对，校正的16S扩增子样本与其配对的WGS样本之间的Bray-Curtis距离变得稳定，且显著低于WGS的组内距离。 两个图像共用X轴。 通过双尾Wilcox秩和检验计算p值，**表示p < 0.05，***表示p < 0.01。
图3. 数据集1的622个WGS：16S扩增子配对样本的beta多样性。 a. 16S扩增子和WGS方法的总体功能模式是同构的，但在PC1和PC2分布上存在明显差异。 b. 由WGS和16S扩增子计算的Bray-Curtis距离高度相关（Pearson相关R = 0.86，p < 0.01）。 c. Meta-Apo使用15个配对样本进行训练，将16S扩增子样本的预测功能谱与WGS样本的预测功能谱进行比对，从而使校正的功能谱的PC1和PC2比原始的未校正的16S扩增子样品更接近WGS样品。 d. WGS：16S扩增子对的ΔPC显著降低。 PCoA使用Bray-Curtis距离计算主坐标。 p值是通过双尾配对的Wilcox秩和检验计算得出，***表示p < 0.01。
大规模16S扩增子预测功能图谱的校正
为了评估大规模16S扩增子功能图谱的校正性能，我们将Meta-Apo样本扩展至5,350 个V3-V5 16S rRNA扩增子样本和与2,354 个WGS样本（数据集3，同数据集1一样从四个身体部位收集，并使用相同的方法处理序列；表1）。该数据集尽管是来自于相同的健康宿主队列，并由同一研究进行测序（人类微生物组计划[2]； HMP），但WGS和16S扩增子样品并未配对。另外我们发现，无论选择何种测序策略[18]，由WGS和16S扩增子所得出的物种结构组成是一致的，但在功能图谱上则有显著差异（图5a； PC1双尾 Wilcox秩和检验p < 0.01； PC2双尾 Wilcox秩和检验p < 0.01）。例如，在功能图谱上，肠道部位的16S扩增子与口腔中WGS聚类在一起，口腔等相同部位的样本会按照不同的测序策略分离，即身体部位在人类微生物组的功能格局中占主导地位[2，19]。之后，我们使用Meta-Apo对所有扩增子样本的预测功能图谱进行校正。经β多样性的分析证明，Meta-Apo校正后的16S扩增子和WGS样本之间功能谱的偏差大大降低（图5b； PC1双尾 Wilcox秩和检验p = 0.20； PC2双尾 Wilcox秩和检验p = 0.03）。
为了测试对不同可变区16S数据集的校正效果，我们也将Meta-Apo应用于表1中数据集4的2,186个V1-V3区16S扩增子样本[2]。使用数据集2的WGS：16S扩增子对做训练样本（训练对N = 15），Meta-Apo也可以有效地提高16S扩增子的功能谱重建的准确性。因此，Meta-Apo普遍适用于16S rRNA基因的多个可变区域。
图4. KEGG BRITE 第3级代谢通路层次上主要功能图谱的比较
图5. 来自数据集3的2,655个WGS样本和5,350个16S扩增子样本的功能beta多样性。 a . 16S扩增子和WGS方法获得的功能模式在PC1和PC2分布上有显著差异。 b. Meta-Apo使用15个配对样本进行训练，将扩增子样本的预测功能图谱与WGS样本的预测功能图进行比较，与原始的未经校正的扩增子样品相比，校正后的扩增子样本的功能谱的PC1和PC2更接近WGS样本。PCoA使用Bray-Curtis距离计算主坐标。通过双尾Wilcox秩和检验计算p值，***表示p <0.01< span="">
功能图谱的校正可实现WGS和16S扩增子测序之间的跨平台比较，并提高疾病状态分类的准确性
鸟枪WGS具有较高的分辨率和较低的序列扩增偏差，与16S扩增子相比，可以提供更稳定的微生物组的疾病检测和分类[13，14]。 由于WGS成本较高，难以被广泛用于商用或日常家用微生物组检测。利用数据集5，我们发现，基于WGS所构建的疾病检测数据模型，使用Meta-Apo校正的16S扩增子的功能谱也可以获得较高的分类精度，而未校正的16S扩增子功能谱则无法实现。
我们评估了Meta-Apo跨平台比较和状态鉴定的性能。尽管两种测序方法均能够揭示健康和疾病微生物组之间的差异，但WGS和未经校正的16S扩增子样本预测的功能图谱在β多样性上有显著差异（图6a）。这是因为测序类型的影响（Adonis R2）超过了疾病状态的影响（图6b，左图）。16S扩增子样本使用Meta-Apo进行校正可减少由测序策略的变化引起的功能图谱的偏差（图6c），并让疾病状况的影响占主导地位（图6b，右图）。因此，Meta-Apo可以进行跨16S扩增子和WGS平台的微生物组诊断。
为了定量评估Meta-Apo跨平台进行诊断的优点，我们使用微生物组搜索引擎（MSE）[20，21] ，采用不同数据对牙龈炎进行预测。首先通过18个WGS样本的功能特征构建数据库，然后分别使用123个原始16S及其校正后的16S扩增子预测疾病状态。未经校正的16S扩增子样本的疾病预测准确性较低，仅为65.04％（F1值 = 0.6446），主要原因是检测的敏感性较低（召回率= 0.4756；图6d）。相反，通过Meta-Apo校正后，疾病分类的准确性提高到95.12％（F1值 = 0.9570），同时检测疾病的敏感性也大大提高了（召回率= 0.9390）。 因此，对于同时可用16S扩增子和WGS类型数据的研究，Meta-Apo提供了一种跨平台微生物组分析策略，可以显著提高状态分类的性能。
图6. 用未校正和校正后的功能图谱跨平台比较口腔微生物组的健康和牙龈炎。 a. 16S扩增子和WGS方法获得的功能模式是不同的。 b. 使用Adnois检验比较抽样因素的影响大小。 c. Meta-Apo将16S扩增子样本的预测功能谱与WGS样品的预测功能谱对齐。 d. 基于MSE的WGS样本分类，对原始样本和Meta-Apo校正后的16S扩增子样本进行健康状态分类。 使用Bray-Curtis指标计算Adonis测试和PCoA的距离。
结论和讨论
微生物组结构和功能的大规模分析迅速地推动了微生物组的研究[20]。目前，基于16S扩增子样本的实验和计算成本较低，尤其是在生物量不足以用鸟枪WGS的情况下，基于16S扩增子的方法更加实用。因此，16S扩增子准确、可靠地重建微生物组功能的能力，将极大地加速大规模微生物组数据的挖掘，具有很高的应用价值。Meta-Apo能够综合16S扩增子测序的较低成本和WGS的较高精确度两方面的优势，使大规模的微生物组研究受益匪浅。
然而，16S rRNA的功能分析的准确性也取决于其物种注释的分辨率[23]和可用的参考全基因组数量[7]。目前，16S rRNA基因可用的参考全基因组数量大约为20,000 [24]，主要来自人体微生物群落（例如PICRUSt 2的参考数据库）。这比已知的16S rRNA基因的数量低一两个数量级[25]，而且在环境微生物群落中16S的功能图谱的校正更加困难。与此同时，大规模培养[26、27]和单细胞测序技术[28]的发展正在迅速提高新型微生物基因组（及其相关的16S rRNA基因）样本的数量。我们认为对于之前和新兴微生物组项目，借助Meta-Apo等新工具，16S扩增子的测序和分析策略将继续为微生物群的功能研究做出贡献。
参考文献
1.Integrative HMPRNC. The integrative human microbiome project.Nature. 2019;569(7758):641–8.
2.Human Microbiome Project C. Structure, function and diversity of healthy human microbiome. Nature. 2012;486(7402):207–14.
3.Knight R, Vrbanac A, Taylor BC, Aksenov A, Callewaert C, Debelius J, et al. Best practices for analysing microbiomes. Nat Rev Microbiol. 2018;16(7):410– 22.
4.Abubucker S, Segata N, Goll J, Schubert AM, Izard J, Cantarel BL, et al. Metabolic reconstruction for metagenomic data and its application to human microbiome. PLoS Comput Biol. 2012;8(6):e1002358.
5.Morgan XC, Huttenhower C. Chapter 12: human microbiome analysis. Comput Biol. 2012;8(12):e1002808.
6.Su X, Jing G, Zhang Y, Wu S. Method development for cross-study microbiome data mining: challenges and opportunities. Comput Biotechnol J. 2020;18:2075–80.
7.Langille MG, Zaneveld J, Caporaso JG, McDonald D, Knights D, Reyes et al. Predictive functional profiling of microbial communities using 16S rRNA marker gene sequences. Nat Biotechnol. 2013;31(9):814–21.
8.Douglas GM, Maffei VJ, Zaneveld JR, Yurgel SN, Brown JR, Taylor CM, et PICRUSt2 for prediction of metagenome functions. Nat Biotechnol. 2020; 38(6):685–8.
9.Asshauer KP, Wemheuer B, Daniel R, Meinicke P. Tax4Fun: predicting functional profiles from metagenomic 16S rRNA data. Bioinformatics. 31(17):2882–4.
10.Jun SR, Robeson MS, Hauser LJ, Schadt CW, Gorin AA. PanFP:
a)pangenome-based functional profiles for microbial communities. BMC Res Notes. 2015;8:479.
11.Walker AW, Martin JC, Scott P, Parkhill J, Flint HJ, Scott KP. 16S rRNA gene- based profiling of the human infant gut microbiota is strongly influenced by sample processing and PCR primer choice. Microbiome. 2015;3.
12.Bonnet R, Suau A, Dore J, Gibson GR, Collins MD. Differences in rDNA libraries of faecal bacteria derived from 10-and 25-cycle PCRs. Int J Syst Evol Microbiol. 2002;52:757–63.
13.Wirbel J, Pyl PT, Kartal E, Zych K, Kashani A, Milanese A, et al.Meta-analysis of fecal metagenomes reveals global microbial signatures that are specific for colorectal cancer. Nat Med. 2019;25(4):679–89.
14.Jackson MA, Verdi S, Maxan ME, Shin CM, Zierer J, Bowyer RCE, et al.Gut microbiota associations with common diseases and prescription medications in a population-based cohort. Nat Commun.2018;9(1):2655.
15.Minoru K, Susumu G, Yoko S, Miho F, TJNAR M. KEGG for integration interpretation of large-scale molecular data sets. 2012;40(D1):D109–14.
16.Huang S, Li R, Zeng X, He T, Zhao H, Chang A, et al. Predictive modeling of gingivitis severity and susceptibility via oral microbiota. ISME J. 2014;8(9): 1768–80.
17.Franzosa EA, McIver LJ, Rahnavard G, Thompson LR, Schirmer M, Weingart G, et al. Species-level functional profiling of metagenomes and metatranscriptomes. Nat Methods. 2018;15(11):962–8.
18.Rausch P, Ruhlemann M, Hermes BM, Doms S, Dagan T, Dierking K, et al. Comparative analysis of amplicon and metagenomic sequencing methods reveals key features in the evolution of animal metaorganisms. Microbiome. 2019;7(1):133.
19.Turnbaugh PJ, Hamady M, Yatsunenko T, Cantarel BL, Duncan A, Ley et al. A core gut microbiome in obese and lean twins. Nature. 2009; 457(7228):480–4.
20.Su X, Jing G, McDonald D, Wang H, Wang Z, Gonzalez A, et al. Identifying and predicting novelty in microbiome studies. MBio. 2018;9:6.
21.Su X, Jing G, Sun Z, Liu L, Xu Z, McDonald D, et al. Multiple-disease detection and classification across cohorts via microbiome search. mSystems. 2020;5(2):e00150–20.
22.Ravel J, Gajer P, Abdo Z, Schneider GM, Koenig SS, McCulle SL, et al. Vaginal microbiome of reproductive-age women. Proc Natl Acad Sci U S A. 2011; 108(Suppl 1):4680–7.
23.Yarza P, Yilmaz P, Pruesse E, Glockner FO, Ludwig W, Schleifer KH, et al. Uniting the classification of cultured and uncultured bacteria and archaea using 16S rRNA gene sequences. Nat Rev Microbiol. 2014;12(9):635–45.
24.Haft DH, DiCuccio M, Badretdin A, Brover V, Chetvernin V, O'Neill K, et al. RefSeq: an update on prokaryotic genome annotation and curation. Nucleic Acids Res. 2018;46(D1):D851–60.
25.QuastC,PruesseE,YilmazP,GerkenJ,SchweerT,YarzaP,etal.TheSILVA ribosomal RNA gene database project: improved data processing and based tools. Nucleic Acids Res. 2013;41(Database issue):D590–6.
26.Zou Y, Xue W, Luo G, Deng Z, Qin P, Guo R, et al. 1,520 reference genomes from cultivated human gut bacteria enable functional microbiome analyses. Nat Biotechnol. 2019;37(2):179–85.
27.Forster SC, Kumar N, Anonye BO, Almeida A, Viciani E, Stares MD, et al. A human gut bacterial genome and culture collection for improved metagenomic analyses. Nat Biotechnol. 2019;37(2):186–92.
28.Xu J, Ma B, Su XQ, Huang S, Xu X, Zhou XD, et al. Emerging trends for microbiome analysis: from single-cell functional imaging to microbiome data. Engineering. 2017;3(1):66–70.
29.Truong DT, Franzosa EA, Tickle TL, Scholz M, Weingart G, Pasolli E, et al. MetaPhlAn2 for enhanced metagenomic taxonomic profiling. Nat Methods. 2015;12(10):902–3.
30.McDonald D, Price MN, Goodrich J, Nawrocki EP, DeSantis TZ, Probst A,
a)et al. An improved Greengenes taxonomy with explicit ranks for ecological and evolutionary analyses of bacteria and archaea. ISME J. 2012;6(3):610–8.
31.Jing G, Sun Z, Wang H, Gong Y, Huang S, Ning K, et al. Parallel-META 3: comprehensive taxonomical and functional analysis platform for efficient comparison of microbial communities. Sci Rep-Uk. 2017;7:40371.
32.Chen IA, Chu K, Palaniappan K, Pillay M, Ratner A, Huang J, et al. IMG/ M v.5.0: an integrated data management and comparative analysis
a)system for microbial genomes and microbiomes. Nucleic Acids Res. 2019;47(D1):D666–77.
33.Dixon P. VEGAN, a package of R functions for community ecology. J Sci. 2003;14(6):927–30.
34.R-Core-Team (ed.). R: a language and environment for statistical computing. Vienna: R Foundation for Statistical Computing;2013.
35.Dray S, Dufour AB. The ade4 package: implementing the duality diagram ecologists. J Stat Softw. 2007;22(4):1–20.