人类肠道微生物的新基因组蓝图

Alexandre Almeida1,2*, Alex L. Mitchell1

 中文由于文编辑完成翻译；

人类肠道微生物群的组成与健康和疾病有关，但要理解其生物学作用，需要学习单个微生物物种的相关知识。尽管进行了广泛的培养和测序工作，人类肠道微生物群的完整细菌库仍未确定。本文通过从11850个人体肠道微生物群中重建92143个宏基因组组装的基因组，鉴定出1952种未培养的候选细菌。这些未培养的基因组极大地扩展了人类肠道微生物群的已知物种组成，系统发育多样性增加了281%。尽管与参考分离基因组相比，新发现的物种在充分研究的种群中并不普遍，但它们将正在研究的非洲和南美样本的分类提高了200%以上。这些候选物种编码了数百个新的生物合成基因簇，并具有可以解释其复杂本质的独特功能能力。本研究揭示了未经培养的肠道细菌的多样性，为肠道微生物群的分类和功能特征提供了前所未有的分辨率。

过去十年来，对人体肠道微生物群的研究表明，微生物和宿主之间的相互作用与各种重要的医学表型有关【1-2】。鸟枪法宏基因组分析方法可以从复杂的微生物群落中推断出分类和功能信息，为旨在了解其在人类健康和疾病中的潜在作用的表型研究提供指导。

然而，用于宏基因组数据集分析的各种策略依赖于高质量的参考数据库【3】。这突出了收集广泛和特征良好的参考基因组的必要性，例如来自人类微生物组项目（HMP）【4-5】和人类胃肠道细菌基因组收集（HGG）的参考基因组【6-8】。虽然进行了新一轮的培养工作，但是肠道生态系统中未分类微生物多样性仍然显著但其多样性程度尚未确定【6,8-11】。虽然这些未知的群落成员可能由于各种原因（例如，由于生长介质中缺乏营养或肠道中营养不足）而无法采用当前的培养策略，但它们很可能发挥尚待发现的重要生物学作用。因此，获得具有代表性的基因组和肠道微生物群分离物的综合目录对于获得新的机理见解是至关重要的。

免培养法和无参考法已被证明是发现和表征物种的成功策略【12-16】。最常见的方法是将鸟枪法基因组reads序列从头组装成contig，并根据序列覆盖范围和四核苷酸频率将其放入不同的排序箱中【15】，该过程可恢复潜在基因组（称为“宏基因组组装基因组”（MAGs）【13,17-19】。一些研究利用这些方法重建了大量的MAGs，其中最显著的是恢复了数千个基因组，这些基因组揭示了对进化树的新见解【16】。

本研究中，从11,850个人类肠道宏基因组装配体中生成并分类了92,143个MAGs，以扩展对肠道相关微生物群多样性的理解。研究中，发现了1,952种未培养的细菌，并研究了它们与特定地理环境的联系以及它们独特的功能能力。这为了解这个未表征的细菌群落中哪些物种和功能在人类肠道环境中可能扮演着未被充分认识的角色提供了新的见解。

未培育物种的大规模发现

为了对人体胃肠道微生物群进行全面的表征，从75个不同的研究中检索了13133个人体肠道宏基因组数据集（补充表1和扩展数据图1）。样本主要来自北美（n =6869, 52%）或欧洲（n =4716, 36%），反映了当前人类肠道微生物组研究的地域偏差。大部分含有可用元数据的数据集都来自疾病患者（n = 4323, 33%）和成年人（n = 3053, 23%）。

在用SPAdes【20-21】组装之后，13133个宏基因组装体中的11850个产生了可以通过MetaBAT【15】进行基因组分箱的contigs，共产生242,836个排序箱。根据基因组完整性和污染程度，用CheckM 【22】对每个箱的质量进行评估（扩展数据图2）。基于这些指标，获得了40029个具有完整性> 90%和污染水平< 5%的MAGs（以下称为“接近完整”）【16】。此外，还生成了65671中等质量MAGs（完整性≥50%和污染水平< 10%）,其中52347 个MAGs的质量得分（QS）【16】 > 50（QS定义为完整性−5×污染水平）【23】。用两种独立的装配/分箱方法评估了MAGs的稳健性（见补充讨论和扩展数据图3），这表明MAGs具有高度的可重复性，其独立于用于装配或分箱的方法【24-25】。

由于CheckM无法评估非原核生物基因组，因此本文分别研究了代表已知真核生物或病毒序列的排序箱数量（参见补充讨论和补充表2）。然而，对于主要分析集，主要集中在39891个接近完整的MAGs上，CheckM将其分解为细菌谱系（补充表3），排除了其余的139个属于古菌的MAGs。为了确定从纯细菌培养中分离出的物种（即分离基因组）的MAGs数量，将每个MAG分配给一个人类特异性参考（HR）数据库，该数据库由2468个分离基因组组成，这些基因组来自HMP目录和HGG（图1）【8】。该数据集包括956个个体物种（553个专门从胃肠道培养的物种），这些物种根据先前报道的物种描述基因组阈值进行定义（≥60%的基因组的平均核苷酸一致性≥95%）【26，27】。为了扩大分类潜力，将MAGs与RefSeq中8778个完整的细菌基因组进行了比较（图1b）。在39891个MAG中，根据≥60%的MAG的平均核苷酸一致性≥95%（ANI）的标准，将26898个分配给HR数据集，12970个分配给RefSeq。HR中不同分类类群的覆盖度较广（扩展数据图4），其中最常见的3个基因组分别属于瘤胃球菌（n = 1255）、泽利斯腐败杆菌（n = 1142）和直肠真细菌（n = 839）。它们都是已知的人类肠道定植体，证实了这些物种是肠道微生物群的共同成员【28】。

随后，重点研究了11888个未分配给HR或RefSeq的接近完整的细菌MAGs（30%）（图1b）。在物种估计水平上（见补充讨论和扩展数据图5）对MAGs进行去重复，产生了1175个接近完整宏基因组物种（MGS），根据CheckM估计，这些宏基因组物种的中位完整性为96.5%（四分位距，IQR = 93.8%-98.4%），污染水平为0.8%（IQR = 0.00%–1.5%）。使用1175 个MGS的数据集，通过将分析扩展到QS> 50的接近完整和中等质量的细菌MAGs（n = 92143，扩展数据图2），评估了最初采集的人类肠道MAGs中仍未被分配的MAGs数量。因此，又鉴定出893种细菌物种，其中位完整性为77.8% （IQR = 68.9%–85.8%），CheckM污染水平为1.1%（IQR = 0.2%–2.0%），后文称之为中等质量MGS。因此，加上 1175个接近完整的MGS，本文共发现2068个MGS（扩展数据图6），代表了人类特异性和高质量参考数据库所缺乏的高质量细菌基因组（更多关于MAG质量评估的细节详见补充讨论）。

物种表征及分布

在确认了人类肠道中2068个MGS后，本研究对其进行了系统分类，并将分析扩展到更全面的参考数据库。通过用针对UniProtKB的蛋白质搜索来补充CheckM的系统发育推断方法，为每个MGS分配了最有可能的分类谱系【29】。这种方法利用多个标记基因和蛋白水平匹配，与各种分析工具所使用的方法类似，与传统的单标记基因分类（例如基于16S rRNA基因）相比，该方法提供了一种更可靠的分类分配方法【30-32】。利用物种水平阈值（≥60%的蛋白质具有≥96%的氨基酸一致性），可以看出94%的MGS（n=1952）与UniProtKB内的任何分离基因组都不匹配，因此代表未培育的候选物种【26,33】。截至2018年8月，这1952个未分类的MGS中（UMGS），74%对应于全新的基因组（见补充讨论和补充表4）。98%和94%的UMGS都能在门和纲级进行划分，91%的UMGS属于已知的目（图2a）。有趣的是，26%的UMGS无法在科水平上进行分配，而近一半（40%）的UMGS不能被划分为已知属，这意味着大部分UMGS可能属于新的科和/或属。最常见的三个科是红蝽菌科（20.6%）、瘤胃菌科（9.9%）和消化链球菌科（7.4%），而排在前三位的属是柯林斯菌属（17.7%）、梭菌属（7.3%）和普氏菌属（4.4%）。这些数据表明，尽管已知肠道微生物群的定植体，这些分支仍然含有相当多的未培养微生物的多样性。梭菌属被认为是高度多系的，最近的系统发育估计表明，该类群可能跨越29科121属【34】。因此，对该属的许多未培养物种的检测可能反映了当前的分类局限性，而不是生物信号。

为了确定每个肠道微生物组中未培养候选物种的流行度和丰度，将原始13,133宏基因组数据集的原始reads与UMGS收集数据进行了比较。通过考虑基因组覆盖水平、平均测序深度和均匀性（扩展数据图7），来确定每个基因组中样本数量进而估计流行度。半数UMGS至少在31个宏基因组样本中被发现（扩展数据图7c）。最常见的UMGS属于瘤胃菌科和粪杆菌属，其中的大部分成员来自梭菌纲（图2b）。

为了将这些未培养的物种与已知的人类肠道细菌定植体放在一起，我们随后将UMGS定位于HR数据库中的肠道特定物种中，以下简称人类肠道微生物参考（HGR）。基于使用specI提取的40个标记基因，建立了1,952个UMGS和553个HGR基因组的最大似然系统发育（图3a）【32】。系统发育分析显示，UMGS基因组在总分支长度的基础上，使人类肠道细菌谱系的多样性增加了281%，其中Firmicutes门内增幅最大（图3b）。从放线菌门（尤其是柯林斯菌属）中检索到一些具有高度系统发育相似性的未培养基因组。这表明，与人类肠道细菌的其他分支相比，该类群中物种和属之间基于基因组的界限更为脆弱。值得注意的是，UMGS包括属于蓝细菌（嗜胃腺菌）、糖类细菌、螺旋体和疣微菌的基因组。这些很可能与从人类肠道中更稀有或更难培养的分支，因为在HGR数据库中没有一个具有代表性的分离基因组。

随后，将每个UMGS和HGR基因组的流行度和丰度与样本的地理来源进行关联，以推断出任何关联（图4）。研究中，调查了在相对丰度> 0.01%的情况下，每个物种在不同大陆上发现的样本数量（图4a）。在大多数取样群体中，UMGS的流行度低于HGR基因组，这可能说明了在以前的基因组研究中没有检测到UMGS的原因。然而，与HGR基因组相比，UMGS在非西方生活方式的非洲和南美的正在研究样本中更为频繁（图4a）。这对于含75和120 个UMGS的子集尤为明显，它们分别以大于0.01%的丰度存在于非洲和南美20%以上的样本中（图4b）。这仅仅是6个和16个HGR基因组的情况，这表明新发现的一些UMGS更好地展现了这两个采样不足群体中少量样本中的肠道微生物多样性。

为了进一步评估UMGS对完整宏基因组数据集分类的贡献，评估了可分配给HR、RefSeq和UMGS数据集的reads百分比。利用所有可用基因组（HR, RefSeq，加上所有UMGS），观察到中位数分类为72.8%（IQR = 65%–81.1%）。这比仅使用HR组成的数据库提高了23%，比使用HR和RefSeq组成的数据库提高了17%。由于所采集的UMGS是HR数据库中肠道物种数量的三倍多，这一温和的增长数据再次表明，与肠道分离基因组相比，大多数未培养的微生物在大多数样本中以较低丰度存在。

根据地理来源对数据进行划分后，来自非洲（n = 21）和南美（n = 36）的少量数据集的reads分配分别提高了215%和278%（图4c）。这证实了一些UMGS在这些特定的肠道群落中更为丰富。为了推断可能未被检测到的物种多样性，根据UMGS数量（以每个大陆样本数量的函数进行检索），构建了累积曲线（图4d）。欧洲和北美群体的覆盖度最高，趋于饱和。相反，在北美和欧洲以外的样本中，仍以一致的覆盖度检测到新的未培养物种。这些结果强调了对采样不足的区域进行取样在继续揭示人类肠道微生物群的全球多样性方面的重要性。

独特功能库

通过对2,505种人类肠道物种 （1952种UMGS和553种HGR）的研究，对集体肠道菌群进行了全面深入的功能表征。使用antiSMASH，筛选了同时在UMGS和HGR中编码的次生代谢产物生物合成基因簇（BGCs）（补充表5）【35】。本文检测到编码囊肽、非核糖体肽合成酶（NRPSS）和细菌素的200个以上BGCs（扩展数据图8a）。值得注意的是，在UMGS和HGR中检测到的总BGCs中，分别有85%和70%是新的簇（即在生物合成基因簇数据库MIBiG的最小信息中没有阳性匹配；扩展数据图8b）。这表明了肠道微生物群可能产生了许多未被发现的天然化合物，可能在未来的研究中用于抗菌和/或生物技术。

接下来，本文应用互补的方法来识别UMGS和HGR基因组之间最显著的特征。首先，从预测的蛋白质编码序列中，我们使用InterProScan【36】生成注释，这些注释被翻译成1199种基因组特性【37-38】（GPs）和115个元基因组基因本体（GO）【39-40】slim术语，这是对来自元基因组数据的GO注释的总结分类【41】。根据检测到的参与该属性的蛋白质数量，判断每个GP（预测将被编码到基因组中的功能属性）是否存在、部分存在或不存在。同时，使用GhostKOALA生成KEGG Orthology （KO）注释来跟踪UMGS和HGR集合中特定功能类别的差异丰度【42】。在全球范围内，通过根据分类组成对GPs库进行分析，可观察到GPs在门级分离良好（ANOSIM R = 0.42, P < 0.001），其中拟杆菌门和变形杆菌门显示非常独特的功能特征（图5a）。我们进一步研究了每个门内UMGS和 HGR基因组的分离，发现放线菌、厚壁菌门、变形杆菌和特内酯菌之间有很强的分化（相似性分析R≥0.30，扩展数据图9a）。特别是，检测到来自放线菌、厚壁菌门、变形杆菌和特内米克菌的UMGS基因组中分别富集了182、207、115和68个GPs（经过卡方检验，调整后P值<0.05），类杆菌组中只有8个功能富集。在这四个最具特色的门中，UMGS中有21种功能持续富集，其中包括与铁代谢和转运有关的特性（扩展数据表1）。

随后，通过评估GO和KO注释的频率，能够应用一种定量方法来比较HGR和UMGS功能库。总的来说，UMGS和HGR基因组间参与碳水化合物代谢的KEGG通路差异最大，表明培养物种和未培养物种之间存在明显的代谢亲和性（扩展数据，图9b）。在GO方面，UMGS中含量较少的基因（Wilcoxon秩和检验，调整后P值< 0.05）与抗氧化和氧化还原功能的相关性特别强（图5b），对活性氧的耐受性较低。如果UMGS对应的是对环境氧气更敏感的严格厌氧菌，那么它们很可能更难分离和培养。相反，根据GP结果，还观察到了UMGS基因组中编码铁硫和离子结合的基因的富集以及其它各种功能。在缺氧条件下，同时有利于硫和氮配体的亚铁形式（Fe²⁺）最为丰富【43】。铁硫结合基因的富集再次表明，UMGS可能更适应于在低氧压力或高铁浓度下的胃肠道特定生境（niches），上述两种环境下都会产生高水平的铁离子【43】。总的来说，这些数据表明，这里描述的未培养物种具有特殊的可解释其复杂本质的功能，同时，也提高了人们对目前从纯细菌培养中获得的参考基因组集合中未充分表征生物特征的认识。

讨论

人类肠道微生物群是研究最多的微生物环境之一，但技术和实践的限制阻碍了研究人员分离和测序每个组成物种的能力。宏基因组学方法提供了获取未培养微生物多样性的途径，本文使用这些方法发现了1952种未培养的候选细菌。这些推定种中几乎有一半不能划分入现有的种属中，这表明代表了相当程度细菌多样性的微生物仍未培养。这一资源进一步扩展和补充了最近的一项研究，该研究发表于我们的论文审阅期间，旨在调查人体微生物群落未被探索的多样性【44】。

本研究获取了92143个MAGs的代表性基因组，这些基因组是由人类肠道组件重建的，人们也能够对73%的潜在reads数据进行分类。然而，培养和从头分析方法都固有地偏向于最丰富的生物体，这意味着可能仍然会遗漏数量一贯较少的物种。此外，在目前的研究中，非洲和南美洲等地理区域的代表性严重不足。因此，将这一分析扩大到世界范围内的大群体将是完全了解人类肠道微生物群景观的必要条件。此外，因为有更全面的参考数据库和完善的标准和工具，研究工作主要集中在细菌基因组上。然而，正如这里所展示的，从肠道微生物群中产生的宏基因组组合包括许多需要进行更彻底研究的其他生物体，如古生菌、真核生物和病毒。

获得全面的细菌基因组样本能够进行精确且高计算效率的基于参考的基因组分析，从而实现微生物生态系统组成的详细分类。本文研究旨在生成从纯培养物到MAGs的高质量参考基因组，这将成为人类微生物群宏基因组分析的蓝图。在未来的关联和机制研究中利用近2000个额外物种的能力将为研究微生物群对人类健康和疾病的影响带来前所未有的动力。

在线内容

所有方法，额外参考资料，自然研究报告摘要，源数据、数据可用性声明和相关的登录码可在https://doi.org/10.1038/s41586-019-0965-1上查询。

图1. 成千上万的宏基因组组装的基因组与独立基因组不匹配。a. 接近完整（完整性>90%，污染水平<5%）的与人类特定（HR）数据库匹配的宏基因组组装的基因（MAGs）用绿色表示（至少60%的基因组的平均核苷酸一致性≥ 95%），不能分类的宏基因组组装的基因组用灰色表示。比对分数至少为60%的MAGs被放大，如右图所示，并根据与其匹配最佳的HR基因组相关的平均核苷酸一致性（ANI）进行着色。b. 与HR（蓝色）和RefSeq（粉红色）匹配的接近完整的MAGs数量，以及那些与两个数据库中的任何参考基因组都不匹配的MAGs数量。

图2. 最常见的未培养肠道细菌种类的分类。a. 1952种未分类宏基因物种（UMGS）的分类组成，按其在UMGS中所占比例的增加从上到下排列。图例中仅表示了五个最常见的分类群，其余的谱系被归为“其它分类过的分类群”。b. 在13133个宏基因组数据集中，根据基因组覆盖水平、测序深度和均匀度推断出的最普遍UMGS基因组的前20名。每一个物种都是根据类别着色的，括号中是预测的分类单元。

图3. 参考和未培养的人类肠道细菌基因组的系统发育。a. 最大似然系统发生树，包含553个属于人类肠道微生物参考（HGR）基因组，1952个属于未分类宏基因组的物种（UMGS）。根据基因组类型（HGR，接近完整和中等质量的UMGS）标记各演化分支，并在第一外层描述相应的门。第二层的蓝点表示在所有六大洲（非洲、亚洲、欧洲、北美、南美和大洋洲）至少一个样本中发现的基因组，红点表示在非欧洲和非北美样本中唯一检测到的基因组。最外层的绿色条代表了13133个宏基因组数据集中基因组的流行度。b. UMGS提供的系统发育多样性增加水平，其相对于每个门的完全多样性（左），并表示为总绝对分支长度（右）。括号中描述了属于各个门的HGR和UMGS基因组的数量。

图4. 样本和未培养物种的地理分布。a. 存在于至少一个样本中、相对丰度在0.01%以上的每个肠道分离基因组（HGR）或未培养物种（UMGS）的样本量分布（以log10转化）。HGR基因组包括n=33（非洲），348（亚洲），386（欧洲），418（北美），88（南美）和128（大洋洲）个样本。UMGS基因组包括n=232（非洲）、1202（亚洲）、1519（欧洲）、1385（北美）、484（南美）和289（大洋洲）个样本。b. 每个地理区域至少20%的样本中发现的物种数量（丰度>0.01%）。c. 按样本地理位置（非洲：n=21；亚洲：n=1447；欧洲：n=4716；北美：n=6869；南美：n=36；大洋洲：n=24）划分的reads比例的增加百分比，分配给人类特定参考（HR）、RefSeq和UMGS，仅与HR和RefSeq相关。d. 累积曲线，描述了以每个大陆宏基因组样本数的函数检测到的UMGS的数量。数据点为10次自展分析重复的平均值。用渐近回归得到的最优拟合曲线表示每个地理区域。在小图a和c矩形盒的长度代表数据的四分位距（IQR），须分别为自表示第一和第三四分位数的边出发延伸至1.5倍四分位距内的最小值和最大值。

图5. 未经培育的物种具有独特的功能能力。a.按照门进行着色的人类肠道参考（HGR，n=553基因组）和未分类宏基因组物种（UMGS，n=1952基因组）的基因组特性的主成分分析（PCA）。b. 在放线菌、厚壁菌门、变形杆菌和黄粉虫的HGR和UMGS基因组之间差异显著的基因本体（GO）功能。给出了五个具有最大和最小丰度差效应值、错误发现率小于5%的功能。正效应大小表示UMGS基因组的过度表达。与氧化还原功能相关的GO术语以粗体突出显示。

方法

宏基因组数据集。围绕75个不同的研究（补充表1），在欧洲核苷酸档案（ENA）中提取了13133个分类为人类肠道宏基因组的序列。使用mg工具包（https://pypi.org/project/mg-toolkit/）通过ENA API检索每个个体样本的元数据（位置、年龄、健康状况和抗生素使用情况），并通过检查链接到每个项目的出版物（如有）进一步整理这些元数据。只有在原始研究中明确说明，才将样本归类为从健康个体中获得的样本。

从头组装和分箱。每次运行的原始reads都首先使用含有meta 【21】选项的SPADES v3.10.0 【20】进行组装。此后，使用MetaBAT 2 （v2.12.1）【15】以2000 bp的最小contig长度阈值（选项--mincontig 2000）和默认参数对组件进行分箱。通过使用BWA-MEM v0.7.16 【45】将原始reads映射回其程序集，然后使用Samtools v1.5 【46】（Samtools view-sbu，后使用smtools sort）以及MetaBAT 2中的jgi_summary_bam_contig_depth函数计算每个单独contig的相应测序深度，推断出分箱所需的覆盖深度。每个宏基因组组装基因组（MAG）的质量得分（QS）是根据使用lineage_wf的 CheckM v1.0.7工作流程估计的，并计算表示为“完整性−5×污染水平”【22】。利用细菌5S、16S和23S rRNAs的Rfam协方差模型，使用cmsearch方程从INFERNAL v1.1.2（选项-Z 1000 --hmmonly --cut_ga）中检测出核糖体RNAs （rRNAs）【47】。通过所有非重叠采样数的总和推断出总比对长度。如果MAG中包含超过80%的预期序列长度，则认为每个基因都存在。使用细菌的tRNA模型（选项-B）和默认参数，用tRNAscan SE v2.0 [删除的超链接字段]识别转移RNAs（tRNAs）【48】。根据关于宏基因组组装基因组（MIMAG）标准的最低信息定义的标准（高质量：完整性>90%和污染水平<5%，存在5S、16S和23S rRNA基因，以及至少18个tRNAs；中等质量：完整性≥50%和污染水平<10%），将其分为高质量和中等质量的MAGS【23】。由于已知的rRNA区域组装问题【16，,49】，只有240个完整性> 90%和污染水平< 5%的MAGs通过了关于rRNA和tRNA基因存在的MIMAG阈值，因此将最高质量的MAGs称为“接近完整”的【16】。VirFinder v1.1被用来预测使用SPAdes产生的 13133个人类肠道组织中病毒contig的存在【50】。该工具使用基于k-mer的机器学习方法检测病毒和宿主（原核生物）序列之间的区别特征。对于长度≥5 kb的每个contig，计算出病毒序列存在的预期P值，随后使用错误发现率（FDR）阈值为10%的Benjamini-Hochberg方法对用于多重测试的P值进行校正。

将MAGs分配给参考数据库。使用4个参考数据库（HR、RefSeq、GenBank和公共数据集中的MAGs集）对从人体肠道组织中恢复的MAGs集进行分类。HR由从HMP目录（https://www.hmpdacc.org/catalog/）和HGG检索产生的2468个高质量基因组组成（完整性> 90%，污染水平为5%）【8】。从RefSeq数据库中，使用了截至2018年1月的所有完整细菌基因组（n=8778）。就GenBank而言，使用了截至2018年8月存储的共153359个细菌基因组和4053个真核生物基因组（3,456个真菌基因组和597个原生动物基因组）。最后，调查了于2018年8月【13,16-19】从最大公共数据集中可获取的18227个MAGs，包括存储在综合微生物基因组和微生物组（IMG/M）数据库中的数据集【51】。对于每个数据库，使用Mash v2.0中的函数mash sketch将参考基因组转换为具有默认k-mer和草图大小的MinHash草图【52】。然后，用Mash dist计算出每个MAG与参考集合之间的Mash距离，以找到最佳匹配（即Mash距离最小的参考基因组）。随后，每个MAG及其最近的亲属与来自MUMmer 3.23的dnadiff v1.3比对，以比较每对基因组与MAG一致的部分（比对查询，AQ）和平均核苷酸一致性（ANI）【53】。

基因组去重复。为了对未分类的细菌MAG（AQ<60%或ANI<95%）集合进行去重复化，首先用MASH生成高水平的相似聚类【52】。简言之，为这些基因组创建了一个MinHash草图，以便对所有基因组进行全面比较。然后，根据Mash距离关系建立层次聚类，并在0.2的分界点定义单个聚类。随后，每个聚类使用dRep去重复以提取出质量最好并代表单个宏基因组物种（MGS）的MAGs【54】。使用选项-sp 0.9 （90%水平上的初级集群），-sa 0.95（95%水平上的二级集群），-cm larger （覆盖方法：完整性更高）以及 -con 5（污染阈值为5%）进行dRep操作。对于接近完整的MAG，-NC参数设置为0.60（覆盖阈值为60%），而对于Qs>50的中等质量MAGs，该参数更改为0.30（覆盖阈值为30%）。根据上述对接近完整MAGs的标准，使用dRep将2468个HR基因去重复到956个代表性物种中。这些包括专门从人类肠道（称为HGR）中采集的553种物种。

系统发育和分类分析。使用Prodigal V2.6.3 【55】（默认单模）预测基因，并使用specI V1.0 【32】从每个基因组中提取40个通用核心标记基因。通过将标记基因与肌肉v3.8.31连接并比对，建立了系统发育树。仅将特定基因组中缺失的标记基因作为缺失数据保存在比对序列中。使用RAxML V8.1.15和选项-m PROTGAMMAAUTO构建最大似然树【56】。所有的系统发生树都在iTOL中可视化【57】。系统发育多样性根据phytools 软件包根据分枝长度总和进行量化【58】。

使用CheckM和UniProtKB，对每一个MGS进行系统分类【29】。首先，利用CheckM中的函数tree_qa来推断MGS基因组在CheckM内部参考树（包括2052个完成的基因组和3604个草图基因）中的大致系统发育位置。然后，使用Diamond v0.9.17.118的blastp函数，将至少在纲级的分类与根据与UniProtKB（2018_04发行）蛋白质比对推断出的分类分配进行比较【59】。根据先前报道的阈值，如果60%以上的蛋白质中的80%以上序列的氨基酸一致性≥ 96%，则推测在物种水平上有一个积极的影响【26,33】。如果使用不含以下任何术语（“未培养”、“sp”和“细菌”）的完整物种名称标记，则假定UniProtKB内的基因组代表已培养物种。对于那些没有指定物种（UMGS）的MGSs，按照以下标准设置属级边界，如前所定义的：至少50%的蛋白质e值小于1×10^-5，序列一致性大于40%，查询覆盖率大于50%【60】。如果用UniProt预测的分类单元在CheckM参考数据库中丢失，则手动检查完整的谱系以确定最可能的注释。由于可能贴错UniProt条目的标签，如果两个分类之间存在不一致，则保留CheckM分类谱系。最后，利用UMGS基因组在HGR系统树中的定位来进一步解决不一致或错误分类问题。

技术重现性和聚类质量。采用两种额外的方法，对1000个宏基因组的随机子集（补充表1）进行了测试以评估由此产生的MAG的重现性。其中一种方法是用MEGAHIT v1.1.3 【24】组装宏基因组，然后利用MetaBAT 2, MetaBAT 1 和 MaxBin v2.2.461对这些宏基因组进行分箱。然后使用MetaWRAP v1.0中的bin_refinement模块执行优化步骤，以合并并改进三个binner生成的结果【25】。第二种方法涉及一种改进的联合装配方法，其中来自同一研究的单个装配体通过CD-HIT v4.762（选项-c 0.99的cd-hit-est定义了99%的序列一致性阈值）合并且进行去重复操作。然后，宏基因组数据集被映射到使用BWAMEM进行的合并的、非冗余的组装中，以获得MetaBAT 2（选项minContig 2000）的分箱共丰度信息。使用上述混合的Mash和MUMmer工作流程，将用每种方法获得的QS>50的MAGS与用主管道（使用SPAdes单独装配，同时利用MetaBAT 2进行分箱）回收的用于相同1000个数据集的MAGS进行比较（单独装配，带SPAdes，并使用MetaBAT 2分箱）。

为了进一步评估报告的MGS潜在污染水平，使用Matthews相关系数（MCC）分析了包含每个MGS的Mash簇的质量。首先，使用CompareM v0.0.23（https://github.com/dparks1134/CompareM）分析Mash簇内和簇间specI标记基因的平均氨基酸一致性（AAI）。为了能够估计MCC，根据90%、95%和97%三个不同的AAI阈值确定真阳性、假阴性、假阳性和真阴性。对于每一个成对比较，当两个MAGs属于同一个集群并且AAI等于或高于阈值时，认为是真阳性；当他们属于同一个集群，但其AAI小于阈值时，认为是假阴性；当基因组分属于不同的集群但其AAI等于或高于阈值时，认为是假阳性；当基因组分属于不同的集群但其AAI低于阈值时，认为是真阴性。然后，使用mltools 软件包中的mcc函数计算MCCs【63】。MCCs可能的值在-1到1之间，1表示Mash聚类和标记基因AAI之间完全一致。

功能表征。对1952个UMGS基因组和553个HGR基因组的去重复集进行了功能预测分析。利用InterProScan v5.27-66.0（选项-goterms 和-pa）对预测基因进行首次功能性表征【36】。利用antiSMASH 4，选项-knowclusterblast推断出微生物次级代谢产物生物合成基因簇（BGCs）的存在【35】，以确定与生物合成基因簇（MIBiG）储存库最低信息匹配的BGCs的数量。基于InterPro（IPR）条目推导出每个基因的基因本体（Go）【39-40】注释，并使用http://github.com/ebi-pf-team/genome-propertie 中的assign_genome_properties.pl脚本转换为基因组属性（GPs）【37-38】。GhostKOAL被用于生成蛋白质编码序列的KEGG本体（KO）注释【42】。利用成分数据分析工具ALDEX2对UMGS和HGR基因组之间的GO slim和KO项频率进行了差异丰度分析【64】。当我们评估具有不同长度和完整度的基因组时，使用这种方法来考虑总基因计数的差异。aldex.clr函数与从Dirichlet分布中抽取的128个蒙特卡罗实例一起使用，以生成与观测数据一致的每个Go-slim/KO项的概率分布。这些数据随后被转换成对数比的分布，以解释数据的组成性质。利用aldex.effect函数计算各组分布差异的期望值（log2差异的中值）、合并组方差的期望值（log2分散的中值）以及各GO/KO分类丰度差异的标准化效应大小。所用的影响大小度量在概念上与Cohen的d类似，但计算是基于分布本身，而不是基于这些分布的汇总统计数据，从而得出相对稳健和有效的度量【65】。最后，利用aldex.ttest对两个测试组（UMGS 和HGR）之间的GO/KO频率进行非参数Wilcoxon秩和测试。将分类为“是”、“否”和“部分”的GPs分别转换为2、0和1，并用双尾卡方检验检测UMGS基因组中更为普遍的GPs。使用Benjamini-Hochberg方法，对用于多重检测的所有统计测试的预期P值进行了校正。利用FactorMine软件包对HGR和UMGS基因组的GP分布进行了主成分分析（PCA）【66】。根据GP文件之间的Gower距离，采用ANOSIM试验评估基于门和基因类型的分离度。

物种的普遍性和丰度。使用Sourmash v2.0.0a4，对13133份人类肠道宏基因组数据集相对于HR, RefSeq 和 UMGS基因组集合进行了read分类【67】。利用sourmash compute——scaled 1000-k 31——追踪丰度，生成了用于参考（FASTA） 和查询（FASTQ）文件的签名文件。对于每一组参考，创建了一个最低共有祖先数据库（sourmash lca index --scaled 1000 -k 31），其中每个基因组代表一个独特的物种谱系。然后将原始reads与针对每个数据库的sourmash lca gather进行比较。利用BWA-MEM测定物种的普遍性和丰度，通过评估基因组覆盖度、平均测序深度和深度均匀度推断物种的存在。首先，分别计算了缺失覆盖率（100%的基因组覆盖率）乘以log₁₀（平均深度）或深度变异系数（定义为测序深度的标准差除以平均值）对应的深度和变异罚分。这些指标的存在使得可以同时测量覆盖度和深度，因为具有高平均深度（或高深度变化量）但未被很好覆盖的基因组比具有相同覆盖度且测序深度较低的基因组在样本中出现的可能性更小。确定基因组存在的阈值设定为至少60%的最小覆盖度，确定基因组存在的深度和变异罚分都设置为99百分位的最大值（扩展数据图7）。每个物种的相对丰度由唯一映射和正确配对的reads计数（使用samtools view -q 1 -f 2过滤）占总reads计数的比例决定。在每个采样间隔内，基于每个地理区域检测到的UMGS数量的累积曲线自展分析重复10次。使用基于R stats包的SSasymp和nls函数进行渐近回归【68】。

参考文献