NBT：人类微生物组千万基因的参考基因集
人类肠道中整合参考基因集
An integrated catalog of reference genes in the human gut microbiome
Nature Biotechnology [IF:31.864]
06 July 2014 Resource
DOI: https://doi.org/10.1038/nbt.2942
第一作者：Junhua Li1–3,19, Huijue Jia1,19, Xianghang Cai1,19, Huanzi Zhong1,19, Qiang Feng1,4,19
通讯作者：Peer Bork5,17 & Jun Wang1,4,6,11,18
其它作者：Shinichi Sunagawa, Manimozhiyan Arumugam, Jens Roat Kultima, Edi Prifti, Trine Nielsen, Agnieszka Sierakowska Juncker, Chaysavanh Manichanh, Bing Chen, Wenwei Zhang, Florence Levenez, Juan Wang, Xun Xu, Liang Xiao, Suisha Liang, Dongya Zhang, Zhaoxi Zhang, Weineng Chen, Hailong Zhao, Jumana Yousuf Al-Aama, Sherif Edris, Huanming Yang, Jian Wang, Torben Hansen, Henrik Bjørn Nielsen, Søren Brunak, Karsten Kristiansen, Francisco Guarner, Oluf Pedersen, Joel Doré, S Dusko Ehrlich, MetaHIT Consortium
主要作者单位：
1 中国深圳，华大基因(BGI-Shenzhen, Shenzhen, China)
5 德国海德堡，欧洲分子生物学实验室(5European Molecular Biology Laboratory,
Heidelberg, Germany.)
热心肠日报
https://www.mr-gut.cn/papers/read/1045745408
Nature子刊：在人类肠道中鉴定出近988万个微生物基因（2014）
创作：赵弘烨 审核：蓝灿辉 2016年09月27日
原标题：一个人类肠道微生物组全面的参考基因集数据库
许多人类肠道微生物的相关研究都是基于已构建的参考基因集数据库进行的，而已有数据库样本来源不具代表性，构建方法不统一，不利于相互比较；
作者集合了249个新测序样本，和之前已公布的1018个来自三大洲人体肠道微生物样本，构建了一个包含9879896个基因的高质量人类肠道微生物基因集数据库；
通过用新建数据库对比中国和丹麦人群样本发现肠道菌群在物种组成和功能组成上都存在显著差异；
这些差异可能反映宿主不同的饮食、生活环境或遗传因素。
主编评语：2014年，以MetaHIT为班底的团队，特别是华大基因团队，在Nature Biotechnology上公布了截至当时最全的肠道微生物参考基因集，包含9879896个基因，这是令人震惊的数字！
摘要
人体肠道微生物组的许多分析依赖于参考基因集。人肠道微生物组的现有基因集基于来自单个群体或参考基因组或蛋白质序列的样品，其限制了全球微生物组多样性的覆盖范围。在这里，我们将人类肠道宏基因组学项目（MetaHit）的249个新测序样本与1,018个先前测序的样本相结合，创建了一个来自三大洲的队列，其数量至少比以前基因集的队列大三倍。由此我们建立了包含9,879,896个基因的整合基因集（IGC）。该集合包括大多数肠道微生物的接近完整的基因组，其质量也比以前的数据高得多。使用该集合分析来自中国和丹麦个体的一组样本，揭示了特定国家的肠道微生物特征。这个扩展的集合应该有助于定量表征来自肠道微生物组的宏基因组，宏转录组学和宏蛋白质组学数据，以了解其在人类健康和疾病中的变异。
Many analyses of the human gut microbiome depend on a catalog of reference genes. Existing catalogs for the human gut microbiome are based on samples from single cohorts or on reference genomes or protein sequences, which limits coverage of global microbiome diversity. Here we combined 249 newly sequenced samples of the Metagenomics of the Human Intestinal Tract (MetaHit) project with 1,018 previously sequenced samples to create a cohort from three continents that is at least threefold larger than cohorts used for previous gene catalogs. From this we established the integrated gene catalog (IGC) comprising 9,879,896 genes. The catalog includes close-to-complete sets of genes for most gut microbes, which are also of considerably higher quality than in previous catalogs. Analyses of a group of samples from Chinese and Danish individuals using the catalog revealed country-specific gut microbial signatures. This expanded catalog should facilitate quantitative characterization of metagenomic, metatranscriptomic and metaproteomic data from the gut microbiome to understand its variation across populations in human health and disease.
要点 Main
众所周知，我们肠道中的微生物群称为人类肠道菌群，对人类肠道内外的生理和疾病很重要。但是，我们对肠道微生物的遗传和功能多样性的了解还远远不够。通过靶向16S rRNA基因焦磷酸测序来分析越来越多的粪便样品，并且通过宏基因组鸟枪法测序在较小程度上进行粪便样品分析，因为与后者相关的成本更高且数据分析更复杂。与基于培养的方法或单细胞方法相比，短测序读段的宏基因组组装可实现功能分析，并且是获取环境微生物基因组信息的更方便且无偏的方法。但是，来自不同研究的数据分散（最显着的是在MetaHIT和人类微生物组计划（HMP）基因目录中），并且还没有一个能够代表世界各地人类肠道菌群的全面且经过统一处理的数据库。随着测序数据量的增加，尚不清楚肠道微生物组中发现的物种和基因的数量将以何种速度继续增长，以及我们目前的采样和数据分析在多大程度上捕获了肠道微生物组中的常见和稀有的实体(物种和功能基因)。
人肠道微生物组中参考基因的目录对于功能性宏基因组学分析至关重要。测序读长可以映射到目录中，以分析样品的种类和基因含量；在全基因组关联研究中，具有丰富度水平变化的基因可以聚类以揭示疾病标记物；基因含量分析可指导从粪便样品中分离菌株，并在培养过程中可能发生变化之前记录菌株在原始栖息地的基因组信息；随着宏转录组学[8,9]和宏蛋白质组学[10]的普及，基因目录将极大地促进RNA或蛋白质数据的分析。MetaHIT和HMP基因目录分别基于欧洲国家或地区的124个样本（以下简称“欧洲样本”）和美国136个样本（“美国样本”）。但代表性有限，可能包含部分或嵌合基因，可以通过更多的测序数据和最新的处理算法进行扩展或消除。
在这项研究中，我们通过处理来自MetaHIT，HMP的249个新测序样品和1,018个已发表样品以及来自中国的一项大型糖尿病研究，以及511个与肠道相关的基因组测序，建立了人类肠道微生物基因目录细菌和古细菌。可通过我们的网站（http://meta.genomics.cn）免费访问这个9,879,896个基因的非冗余参考目录，并将数据保存在GigaScience数据库中。除了为将来的分析提供扩展的资源外，对目录的研究还表明我们**可能已经达到核心基因内容和功能的饱和范围，但是随着采样的增加，稀有基因将继续被发现**。我们还展示了使用该目录发现肠道菌群的特定种群特征的能力。
结果Results
构建整合基因集
Construction of the integrated gene catalog
在这里，我们通过对丹麦或西班牙成年人的249份粪便样本进行测序，完成了MetaHIT队列研究，从而收集了760份欧洲样本和8,096,991个非冗余基因的目录（图1）。为了创建人类肠道微生物组的洲际基因目录，我们将欧洲样本微生物基因目录与来自368个中国样本和139个美国样本的数据进行了整合（图1a和补充表1）。因为我们使用的标准化和自动化工作流程先前已显示出来，以提高装配，基因预测和冗余去除的质量（在线方法），所以来自美国样本的基因平均要长32.8％，总数要少41.5％。与更新的HMP目录（2013年8月下载）进行比较。与我们相比，更新后的HMP目录具有更多的碎片集（尽管与原始研究相比有所改善）（图1b）。队列的合并产生了一个目录，包含9750788个基因（这里称为三个队列非冗余基因目录（3CGC）），基于来自三大洲的1267个肠道宏基因组（1,070个个体），总计6.4 Tb的宏基因组测序数据（图1和图2）。补充图1a，b），比以前的队列要多得多。
图1. IGC的构建
Figure 1: Construction of the IGC.
（a）数据处理和整合的流程（详见在线方法）。来自欧洲、中国和美国队列的宏基因组测序数据，用MOCAT流程处理以产生它们各自的基因集，再合并为3CGC。根据IMG系统（红色），来自HMP（绿色）的16S rRNA基因序列OTU或3CGC（蓝色）覆盖的序列原核基因组或草图基因组被认为可能是人类肠道来源。通过我们的宏基因组测序数据比对覆盖度，过滤该初始的983个基因组，其产生511个基因组，其基因命名为SPGC。最后，3CGC与SPGC合并以生成IGC。
（b）基因集的基本统计信息。*原始HMP研究报告了一个包含5,183,353个基因的肠道微生物基因集，比2013年8月从HMP网站下载的（http://www.hmpdacc.org/HMGC/）数量多13％; 样本数量为139而不是原文所述的136。**，对中国样本的原始研究创建了基于145，现替代为368个样本的基因集。***，IGC还整合了511个原核基因组。ORF，开放阅读框; N50，此长度或更长的总长度的50％; N90，在此长度或更长的总长度的90％; NA，不适用。
丰富的肠道微生物在3CGC中得到了很好的体现，但是一些低丰度却又常见的微生物却没有被充分覆盖，这可能是因为这些物种的测序深度很低。例如，标记为梭菌D5的菌株。美国国家生物技术信息中心（NCBI）其命名为梭菌D5（但根据核糖体数据库项目（RDP）的16S分类器，我们发现它被分类为链霉菌科的梭菌XlVa而不是梭菌科的梭菌），被整合微生物基因组（IMG）系统列为分离自人类粪便的菌株，我们在HMP的16S rRNA基因数据（在线方法）中53％的粪便样本（n = 325）中检测到该菌株。但是，只有4.9％的基因组被3CGC基因覆盖（补充图1c）。为了确保代表这种低丰度但普遍存在的生物，我们从511个与人类肠道相关的细菌和古细菌菌株的基因组中提取了基因，并在我们的宏基因组测序队列中检测了其基因组（> 90％的基因组累积覆盖率）所有1,267个宏基因组；在线方法）。这导致了一组659492非冗余基因，我们将其称为测序原核基因目录（SPGC）（图1和补充表2）。
我们将SPGC与3CGC合并组成了IGC。IGC包含9,879,896个基因，分别比现有的MetaHIT和重新组装的HMP基因目录分别增加2倍和3倍（图1）。每个样本平均包含762,665个基因，平均贡献469个独特基因。任何两个样本平均共有250,382个基因（762,665个基因中的32.8％）。
整合基因集的质量和完整度
Quality and completeness of the integrated gene catalog
与MetaHIT 2010和HMP 2012基因目录（补充图2a，e）相比，IGC中分别有75.7％和74.1％的基因是新的。对于来自MetaHIT 2010研究的测序读长，我们对IGC的读长段比对MetaHIT 2010目录的读长多了约10％，平均映射率达到79.24％（图2a）。与未整合的欧洲，中国和美国基因目录相比，IGC允许更好地映射其构建群所使用的测序读长（平均80.54％）（补充图3a）。未在IGC的构建中在美国，日本和瑞典进行的三项研究的数据分别占IGC所代表的测序数据的73.67％，81.36％和76.15％（图2b，c和在线方法）。因为在所有原核基因组中基因编码区的百分比约为87％（补充表3），并且平均长度为77个碱基对的估计7.25％的测序读长无法可靠地定位，因为它们仅与基因部分重叠（在线方法），我们在IGC上观察到的映射读长的百分比接近最大可实现的映射速率。此外，基于Chao的丰富度估算（参考文献18）表明，IGC覆盖了采样的肠道微生物组中94.5％的基因含量（图2d），类似于使用基于发病率的覆盖率估算器估算的95.4％（图2d）。
图2. IGC覆盖度
Figure 2: Coverage of the IGC.
（a）MetaHIT 2010研究中总读长（n = 124个样本）比对到MetaHIT 2010和IGC的比对率分布图。绘制的是四分位数范围（IQRs;盒子），中位数（方框中的暗线），最低值和最高值即在第一和第三四分位数（盒子上方和下方的须线）的1.5倍IQR内，以及超出须线的异常值（圆圈） ）。（b）日本样本（Sanger测序，n = 13个样本）和美国样本（Roche 454测序，n = 18个样本）在MetaHIT 2010和IGC上的可比对率，具有同一性标准 ≥ 90％和映射长度 ≥ 100bp
（c）主要瑞典队列（图b中显示）与正常葡萄糖耐量（NGT），葡萄糖耐量降低（IGT）和II型糖尿病（T2D）的作图比率的分布。每个点代表一个样本，根据出生时的国籍进行着色。原始研究无法获得比对率
（d）基于使用Chao2估计的1,267个样品的基因数量的稀疏曲线
IGC中组装的基因与先前目录的比较表明，与IGC中存在的87.8％MetaHIT 2010基因相比，IGC中不存在的MetaHIT 2010基因的12.2％较短、更加分散、并且通常具有未知的分类和功能。IGC（补充图2c，d）。这种差异可能是由于用于生成IGC的方法所致，包括对测序读段进行更严格的质量控制（使用FASTX Toolkit；http：//hannonlab.cshl.edu/fastx_toolkit/），改进的组装软件（SOAPdenovo 1.06）、拼接修订版（在MOCAT流程中）、更具体的基因鉴定（MetaGeneMark）和用于合并基因目录（CD-HIT）的标准化超快聚类算法（图1和在线方法）。类似地，与IGC中存在的76.4％HMP 2012基因相比，IGC中不存在的23.6％HMP 2012基因要短得多，并且与一小部分测序读长对齐（补充图2g）。在目录中共享的基因中，大多数在IGC中更长（补充图2b，f）。
IGC中的物种
Taxonomic representation in the IGC
我们使用3,449个细菌和古细菌的参考基因组对IGC进行了分类注释（补充表3和在线方法）。与先前的研究相似[4,23]，IGC中21.3％的基因可以唯一可靠地分配给门，而16.3％的基因可以分配给属。可以归为属的基因占总测序读数的44.4％（在单个样品中占测序读数的5.3％至78.4％；补充图3b）。
与MetaHIT 2010相比，对于3CGC（不带SPGC的IGC），我们观察到平均而言3CGC的单个基因组覆盖率提高了32.26％（与MetaHIT 2010相比，3CGC的覆盖率提高了-0.71％至+ 80.44％;细菌基因组为87％）（图3a，补充表3和4）。随着队列规模从MetaHIT 2010中的124个扩展到IGC中的1,267个，基因组覆盖率的提高与该属的最大丰度的增加相关（图3a，b）。除两个菌株外，最丰富的细菌属（Bacteroides）的覆盖率不超过10％，而在当前队列中被采样到丰度更高的属（例如，普雷沃氏菌，乳杆菌，肽链球菌，肠球菌和幽门螺杆菌，特别是H. winghamensis）在我们的基因目录中显示出其基因组覆盖率的显着改善。在菌株水平上，由于采样增加，IGC中病原性菌株（如大肠杆菌O157：H7）和干酪起子菌株（如德氏乳杆菌）表现得更好（图3c，补充图4a和补充表4）。肠球菌的分析显示，大多数样品中该属的含量较低，但中国和欧洲样品中偶尔存在的高丰度，加上足够的测序深度，使得IGC中肠球菌的基因组覆盖率提高了70-80％（补充图4b，c）。因此，增加采样数量可能是更深层测序的一种更有效的替代方法，以提高稀有物种的覆盖率。
图3. 3CGC改善基因组覆盖度
Figure 3: Improved genome coverage in 3CGC.
（a）与MetaHIT 2010基因相比，3CGC涵盖的每个菌株基因组的百分比有所提高。仅显示了基因组被MetaHIT 2010或3CGC基因覆盖超过60％的593株（BLASTN v2.2.24，评分≥60且查询的映射长度≥80％）。完整的菌株列表显示在补充表4中。菌株按属分组。每个点代表一个属中的一个菌株。点的大小与同一属中的菌株数量成反比（即，只有一个菌株覆盖超过60％的菌种具有大点）。虚线显示理论最大值87％（细菌基因组的平均基因含量）。（b）MetaHIT 2010队列（n = 124）和IGC队列（n = 1,267）中所见属的最高相对丰度的差异。根据MetaHIT 2010中的相对丰度最大值对属进行排序，得到的x轴标记如a所示。（c）3CGC和MetaHIT 2010中不同乳杆菌菌株的基因组覆盖。
出现在大量样品中的属（出现频率高）往往是那些先前已知居住在人类肠道中的物种（补充图4d，以及补充表5和6）。一个值得注意的例外是葡萄酒发酵中使用的Oenococcus，尚未报道为肠道共刺激物（补充表6），但在当前队列中该基因的注释频率为13.5％。尽管根据IMG（补充表6），与人类肠道无关的属明显多于在肠道中发现的属，但它们仅贡献了相对低发生的基因（补充图4e）。
肠道微生物的功能组成
Functional representation of gut microbes
我们根据《京都基因与基因组百科全书》（KEGG）和基因非监督直系同源基因组（eggNOG）数据库的进化谱系对IGC中的基因进行了注释。我们确定了总共6,980个KEGG直系同源组（KOs）和36,489个eggNOG直系同源组，分别占总测序读数的51.6％和69.3％（补充图3b），分别涉及IGC基因的42.1％和60.4％。
IGC中存在876个KO，但MetaHIT 2010目录中不存在，但由于严格性增加，IGC中不存在36个KO。与Chao2（图2d）和ICE进行的丰富度估算相一致，并且如从MetaHIT 2010到IGC的代谢途径覆盖率的局部而不是全局改善所表明的那样（补充图4f），IGC可能会提供饱和的覆盖率。人类肠道原核生物的功能能力。尽管细菌是肠道微生物组中的优势生物，但与MetaHIT 2010相比，我们在IGC中获得了500个以上的真核生物KO，但是在更高水平的真核生物中的途径（如糖鞘脂生物合成，蛋白聚糖生物合成和二萜类生物合成）仍然不存在（补充图4f和补充表7）。
为了测试IGC在分析宏转录组和宏基因组学数据方面的有用性，我们将最近一项研究的宏转录组测序读长比对到参考基因目录中。与仅使用肠道细菌和古细菌的参考基因组（SPGC）相比，去除与rRNA，tRNA和信号识别颗粒RNA等非编码RNA相对应的基因后，可以将更高百分比的宏转录组读长定位到IGC。3c）。尽管有严格的比对标准（在线方法），但根据IGC比对到每个样品中蛋白质编码基因的读段数量与原始研究的值相关性很好（补充图3d）。同样，使用IGC代替参考基因组（SPGC）可以鉴定更多的KO，尤其是在诸如碳水化合物代谢，细胞过程和信号传导以及膜运输等途径中（补充图3e）。
国家特异的标记
Country-specific signatures
为了证明IGC在人群之间肠道微生物组定量比较中的作用，我们从1267例样本中选择了表型匹配组，包括60名中国南方和100名丹麦健康个体（补充表8和9），并通过与IGC比较。由于两组样本使用的DNA提取方法略有不同，因此在比较之前，我们随机选择了368个中国样本中的11个（补充表10），并使用两种方案提取了DNA来估计由这种差异。使用不同方法从同一样本衍生的宏基因组显示出很高的自相关性和相同的关键特征（补充图5a-d）。在后续比较（在线方法）之前，我们删除了剩余的差异。
通过基于基因（图4a），KO或属谱（补充图5e，f）的主成分分析（PCA），我们可以轻松地将中国人和丹麦人队列分开。与丹麦人群相比，中国人群在基因和属上显示出较低的α-多样性，但在KOs中却没有显着降低（P = 7.82×10−6，P = 1.90×10−6，P> 0.1，Wilcoxon秩和检验），甚至在提取方法和可映射测序读数标准化后（图4b和补充图5c，g，h）。从分类学上看，307个属中的151个显示出中国样本和丹麦样本之间的明显差异（P <0.01，错误发现率（FDR）为0.0048，功效= 0.7，Wilcoxon秩和检验;补充图5i，j和补充表11 ）。例如，丹麦样品通常富含厚壁菌门，包括Oenococcus和其他乳酸菌，而中国样品中的变形菌门丰富（图4c）。
图4. 中丹人肠道微生物差异
Figure 4: Differences between Chinese and Danish gut microbiota.
（a）PCA用于中国（n = 60）和丹麦（n = 100）队列中的肠道微生物基因组成。绘制前两个主要成分。
（b）中国或丹麦队列的样本内多样性（香农指数），根据测序数据等量抽样之前和之后的基因概况计算（在线方法）。无论是否缩小，队列之间的差异都是显着的（Wilcoxon秩和检验，P <0.01;双侧）。
（c）中国和丹麦队列中的前50个差异富集属（P <0.01，Wilcoxon秩和检验，双尾），由相应的门进行着色。
两个队列之间有3,491个KO显着不同（P <0.01，FDR = 0.003，功效= 0.7，Wilcoxon秩和检验;补充图5k和补充表12）。最显着的差异涉及与饮食有关的过程，例如能量代谢，碳水化合物代谢，氨基酸代谢以及辅因子和维生素的代谢，以及与异生质有关的功能，例如膜运输和异种生物的生物降解与代谢（补充图6和图7，补充表13,14,15,16,17,18,19,20,21,22,23,24,25和补充说明）。健康的中国和丹麦成年人之间肠道微生物组代谢潜能的这些差异可能受到饮食（也许是面包，乳制品和维生素）和环境因素（可能是芳香致癌物或氮氧化物）差异的影响（补充图8和补充说明） 。
个体特异的基因
Individual-specific genes
IGC中增加的基因数目无法通过测序和装配错误来解释，因为在基因目录的编制过程中，序列中的此类错误将作为冗余基因（具有> 95％相同性的那些）被消除（图1）。为了确定增加的基因数量的来源（图5），我们在样本大小从50变为1,267时模拟了目录中的基因内容（补充表26）。在超过5％的样本中检测到的基因数量仅略有增加，并达到约320万个基因的饱和度；与MetaHIT 2010（参见图2）相同（图5a），超过50％的受试者中存在的基因数量仍低于300,000。相反，随着样本量的增加，在不到5％的样本中发现基因，尤其是在不到1％的样本中，基因继续增加（图5a）。因此，出现在少数个体中的基因对IGC的扩大大小起了最大作用。尽管它们的出现频率很低，但在确实含有它们的样品中这些基因仍然丰富（图5b）。
图5. 低频基因的功能
Figure 5: Abundance and function of low-occurrence genes.
（a）样本基因稀疏曲线，模拟所鉴定的非冗余基因的数量，样本大小为50至1,267，并根据样本中指示的出现频率作图。低发生基因（即，<1％样本中出现的基因）随着采样的增加而主导基因数量的增加。
（b）包含它们的样品中基因的平均相对丰度在出现频率上作图。通过从所有1,267个样品（1,070个个体）估计的出现频率对基因进行分类，并且在每个箱内计算含有基因的样品中基因的平均相对丰度
（c）eggNOG功能类别中不同发生频率（OF;基于1,070个个体）的基因的分布。基因按出现频率分为三类<1%, 1-50%, >50%
（d）富含个体特异性（OF <1％）基因的功能。在OF <1％的基因和OF≥50％的基因之间比较注释到所列OG的基因的百分比。MGE，移动遗传元件。
从同一HMP个体在不同时间点（平均相距218 d）采集的样本中，低发生基因的丰度和组成在很大程度上是一致的，而来自不同个体的低发生基因则存在显着差异（图6）。这些基因在样品处理或肠道的短暂“乘客”过程中均未被污染。实际上，当用于区分不同个体的样品时，低发生率基因可能比高发生率基因更有效（图6a，b）。
图6. 低频基因的时间稳定性。
Figure 6: Temporal stability of low-occurrence genes.
（a，b）通过从所有1,070个体估计的出现频率对基因进行分类。在每个箱内，基于基因含量的Sorenson指数被估计用于从43个HMP个体（a）和来自94个不同HMP个体的第一时间点（粪便1）样品取得的多个样品的成对比较（b）。在所有出现频率中显示出比其余数据显着更低的Sorenson指数的两组点来自样品763536994-粪便2与相同个体的粪便1和粪便3样品之间的比较; 这个样本似乎是一个异常值。
（c）比较来自在两个时间点取样的43个HMP个体的样品中相同的基因丰度，以计算Spearman相关系数。
使用eggNOG数据库，我们将不到1％的个体与超过50％的个体所发现的基因的功能进行了比较，此处分别称为“个体特异性”和“共同”基因（最新Cell也有此类比较方法）。个体特异性基因在细胞壁/膜/包膜生物发生和DNA复制，重组和修复的类别中适度富集。常见基因的功能丰富，例如信号转导机制，能量产生，碳水化合物转运和代谢以及氨基酸转运和代谢（图5c和补充表27）。
当我们在eggNOG资源中查看确切的直系同源基团（不同生物体中具有相同序列和功能的基因组）时，负责细胞壁成分合成的基因，尤其是肽聚糖和脂多糖的合成在个体特异性集中过高（图5d和补充表28）。我们还观察到在这些个体特异性基因中，与噬菌体相关的蛋白质（包括尾部蛋白质，噬菌体阻遏物和末端酶）的比例提高了八倍。与DNA有关的功能（例如转座酶，内切核酸酶和DNA甲基化酶）富含个体特异性基因（图5d和补充表28），可能与肠道微生物暴露于外源DNA有关。此外，个体特异性基因编码更多的乙酰基转移酶，例如使氨基糖苷类抗生素失活的GCN5相关的N-乙酰基转移酶。这些结果表明，常见基因提供了生存所必需的功能，而个体特异性基因可能反映了对宿主免疫系统，病毒感染，抗生素治疗和肠道微生物组所面临的其他挑战的适应性。
讨论
IGC是用于进一步研究肠道微生物组的综合资源，涵盖了人类肠道中出现频率、丰度和传播持续时间各不相同的菌株。将来增强此目录的工作可能会更多地针对特定目标菌株的高丰度样品，这可能表明偏离健康状况或与特定环境因素有关。由于肠道可能会被食物和饮料中存在的微生物所播种，因此，关于微生物摄入和排泄的定量信息，肠道中菌株的半衰期等，对于可靠地定义肠道指标是必不可少的。与粪便采样相比，诸如结肠镜检查之类的侵入性技术也可能识别出更多的与黏膜相关的微生物。
我们根据IGC对健康的中国和丹麦成年人的两个表型匹配的队列进行分析，发现他们的肠道菌群在营养代谢和异生质毒素解毒的许多方面存在差异，这可能是饮食和环境所致（补充图8）和补充说明）。但是，其他影响，例如宿主遗传学，仍然可能。
低发生率基因对IGC中增加的总基因数起了压倒性的作用，并且可能反映了宿主中遗传，营养和医学因素的独特结合。尽管个体近期没有接受抗生素治疗，但我们观察到在人群水平上可能存在的抗生素抗性基因都富集（丹麦人的青霉素抗性和中文的多药耐药性；补充表21）以及个体特异性基因（例如乙酰转移酶）的富集。和肽聚糖的合成，这突出显示了需要密切监测直接和间接接触抗生素的情况。
肠道噬菌体被认为主要是温带的，但可被诱导进入裂解周期。我们确定了维持溶原性的基因，例如噬菌体阻遏物，以及参与复制和感染的各种基因。噬菌体还可能携带其他个体特异性基因，这些基因已知会改变细菌宿主的代谢并赋予其抗逆性，并且似乎与给定个体稳定相关。非肠道微生物组中的稀有基因是否也富集噬菌体或可能由噬菌体携带的适应性功能尚待探讨。
与人类遗传学领域类似，新的等位基因的搜索已从普通的变为稀有的，我们的数据表明，我们的“其他基因组”（人类肠道微生物组）的分类也进入了识别稀有或单个个体的阶段。特定的基因，而不是共同的和共享的基因。它也已经达到了全球人口之间进行定量比较的阶段。参考基因目录（例如IGC）可以对给定肠道宏基因组的遗传和功能库进行快速多基因组分析，并有助于对其地理，遗传，时间和生理特征的研究。
已经建立了一个网站（http://meta.genomics.cn，针对Safari优化），以更好地显示基因目录的注释信息，并指导对使用我们的数据集和下载特定数据集感兴趣的研究人员。
方法
样品收集和转移
Sample collection and transfer
在MetaHIT联盟下，将249份粪便样品收集到为此目的而提供给参与者的家中的容器中，然后立即转移到-20°C的冰箱中，并在第二天的冷藏箱中冷冻。然后将样品转移至-80°C，并保持在该温度或干冰上。丹麦首都地区道德委员会的所有丹麦志愿者以及医院Univeritari Vall d’Hebron的所有西班牙志愿者均获得了知情同意。所有其他样品以前都已报道过
样品DNA提取
Sample DNA extraction
如前所述从249个新的MetaHIT样品中提取DNA。
为了比较DNA提取方法，BGI的方案与MetaHIT方案相同，只是对于每个粪便样品（最大约1 g），在最初离心后添加25 mg溶菌酶和12.5 mg蛋白酶K以促进细胞裂解。在37℃下孵育1小时以符合溶菌酶的最佳反应温度。
为了评估不同DNA提取方案的影响，我们从中国队列中随机选择了11个粪便样本，并将它们发送到法国国家农学研究所（INRA）。按照MetaHIT方案，我们在INRA的MetaHIT合作者再次从这11个样品中提取了DNA。
HMP使用PowerSoil DNA分离试剂盒（MO BIO Laboratories），根据我们和其他机构的评估（数据未显示），其DNA产量较低。结合HMP样本中较低的α多样性（补充图5l）和非重叠年龄（HMP小于40，MetaHIT大于40），我们在洲际比较中未包括HMP数据。
DNA文库的构建和测序
DNA library construction and sequencing.
按照制造商的说明书（Illumina）进行DNA文库的构建。我们使用与其他地方2所述相同的工作流程来执行簇生成，模板杂交，等温扩增，线性化，封闭和变性以及测序引物杂交。
我们为欧洲队列的249个新MetaHIT样品构建了Illumina文库，插入片段大小为350 bp，然后进行高通量测序，获得约3600万对配对（PE）读数。每个末端的读取长度为90 bp。MOCAT流程从Illumina原始数据中提取了高质量的读长。这些样本中高质量数据的比例平均为89.5％。
我们为来自中国队列的11个随机选择的样本构建了Illumina文库，然后进行高通量测序，以获得约1400万个PE读长或1500万个单端（SE）读长。每个末端的读取长度为90 bp。MOCAT流程从Illumina原始数据中提取了高质量的读数。平均而言，这些样本中高质量数据的比例为87.9％，我们的PE库的实际插入大小范围为311 bp至326 bp。
使用与249个MetaHIT样本相同的方案，在BGI上构建了来自MetaHIT项目的511个欧洲粪便样本的Illumina文库和368个中国粪便样本的文库，并对其进行了测序。HMP测序中心使用相似的协议和平台处理了139个HMP样品。
公共数据重用
Public data used
本研究使用的公共肠道微生物基因组包括：（i）从粪便部位采集的139种HMP样品，可从http://www.hmpdacc.org/HMASM/下载；（ii）368份中国粪便样本，这些样本是从NCBI下载的（登录号SRA045646和SRA050230）；（iii）511个来自MetaHIT项目的欧洲粪便样本，从欧洲生物信息学研究所（EBI）下载，登录号为ERA000116，ERP003612和ERP002061，并在MetaHIT内部共享。所有这些公共宏基因组测序样品均使用MOCAT流程进行了处理，以提取高质量的读长。
用于验证IGC代表性的其他肠道宏基因组学数据包括：（i）来自美国个人的数据，这些数据是从NCBI下载的，登录号为SRA002775；（ii）来自日本个人的数据，这些数据是从EBI下载的，登录号为PRJNA28117；（iii）来自欧洲个人的数据，该数据是从NCBI下载的，登录号为ERP002469。
该项目中使用的两个先前发布的人类肠道微生物组基因目录包括：（i）MetaHIT从124个欧洲人建立的基因目录，可从http://gutmeta.genomics.org.cn/下载；（ii）由HMP建立的基因目录，该目录已于2013年8月从http://www.hmpdacc.org/HMGC/ 下载。
从Gene Expression Omnibus下载了59个样品的肠道转录组数据，登录号为GSE46761。所有这些公共的转录组测序公共样本均由MOCAT流程处理，以提取高质量的读长。
3449个测序原核基因组或原始基因组的收集和质量控制
Collection and quality control of 3,449 sequenced prokaryotic genomes or draft genomes
收集和过滤原核基因组。简要地，下载了2012年2月23日在NCBI和EMBL Bank获得的所有原核基因组，并删除了300个以上重叠群和N50 <10 kbp的基因组。此外，我们删除了40个通用单拷贝标记基因中不到30个的基因组[40,41]。最后，对于具有相同分类法标识符但项目标识不同的基因组，随机选择一个基因组，这导致本研究中使用了3449个基因组。
整合基因目录（IGC）的构建
Construction of the integrated gene catalog (IGC)
使用MOCAT对来自欧洲，中国和美国成年人的粪便样品的Illumina测序读长进行了独立处理（质量控制，人类序列去除，组装，装配修订和基因预测），该方法可以以标准化和自动化的方式处理宏基因组，同时提高质量与在运行时基于参数探索和数据驱动的参数优化的支持的程序使用默认参数相比，装配和基因预测的准确性。我们选择FASTX Toolkit（http://hannonlab.cshl.edu/fastx_toolkit/）进行质量控制，选择SOAPaligner2（参考文献42）来鉴定人的序列，选择SOAPdenovo v1.06（参考文献20）进行装配，并选择MetaGeneMark进行基因预测。我们在MOCAT中使用的配置文件已保存在GigaScience Database中。每个队列中的基因都使用CD-HIT进行聚类。然后合并基因目录，以基于所有1,267个样品（称为3CGC）生成人类肠道微生物基因目录。
收集了3,449个测序细菌和古细菌基因组或原始基因组，并分两步选择了与人肠道相关的原核生物。首先，包括满足这三个条件中的任何一个的菌株：（i）根据IMG（http://img.jgi.doe.gov/cgi-bin/w/main ），该菌株的栖息地是“人的胃肠道”。cgi）（2012年7月24日下载），即该菌株的“身体部位”为“胃肠道”或“分离物”为“人类粪便”。（ii）该菌株的16S rRNA序列与报道的OTU相同通过HMP从粪便身体部位提取。将粪便部位的485个非嵌合HMP OTU与每个菌株的16S rRNA基因进行了同源性分析（1.23.1版），其整体相似度阈值≥99.5％。（iii）标准较弱的3CGC覆盖的基因比例很高。使用BLAT将每个菌株的基因与3CGC进行比对，其标准是重叠≥10％，同一性≥95％。选择其基因的80％以上被3CGC覆盖的菌株。根据这三个标准，我们获得了983例肠道相关原核菌株（图1a和补充表3）。其次，我们的宏基因组测序数据对这983个原核生物基因组进行了累积覆盖率过滤（1,267个样本中90％以上的基因组），以确认它们是人肠道微生物组的一部分。首先，使用相似性 ≥ 90％ 的SOAP2，将每个菌株的基因组或原始基因组与100个样品的测序读长进行比对。将尚未覆盖90％以上基因组的菌株与来自所有1,267个样品的数据进行比对，以进行进一步选择。经过这样的过滤后，保留了511个原核生物，并用于构建与肠道有关的SPGC（图1和补充表3）。
最后，使用CD-HIT将基于宏基因组测序数据（3CGC）的基因目录和基于测序原核基因组（SPGC）的基因目录结合起来，生成IGC（图1）。
1,267个样本的基因目录，注释信息，丰度概况，装配体和预测的开放阅读框架已保存在GigaScience数据库中。
研究低丰度但普遍存在的人类肠道细菌
Investigation on the representation of a low-abundance but prevalent human gut bacterium
选择了最初从人粪便中分离出来的NCBI 556261.HMPREF0240_10201（根据RDP database，Lchnospiraceae | Clostridium XlVa，在NCBI中的Clostridiaceae | Clostridium sp.D5）的基因组作为参考。通过BLAST将3CGC的基因与NCBI 556261基因组进行比对，其标准是具有95％的相似性和90％的重叠。3CGC仅代表了4.9％的基因组。通过对来自HMP的325个粪便样品在16S rRNA基因可变区中的测序数据来评估菌株的发生频率。每个样本的长度超过150 bp的标签通过mothur（1.23.1版）与来自NCBI 556261基因组的16S rRNA基因进行比对，具有97％以上的同一性和90％以上的重叠。HMP粪便样本中有53％带有此物种（平均2.1个标签），这表明它是人类肠道环境中普遍存在但数量较少的物种。通过SOAP2（参考文献42）评估了从1,267个样品中进行的宏基因组测序读取的基因组累积覆盖率，其相似性超过95％。根据从样本组装的ORF覆盖的NCBI 556261基因的最大数目，选择覆盖率最高的样本。
基于参考基因组的系统发育注释
Phylogenetic annotation based on reference genomes
系统发育注释是使用内部流程执行的。（i）我们使用BLASTN（v2.2.24，默认参数，除了-e0.01 -b100 -K 1 -F T）将970万个3CGC基因比对到3449个原核基因组的数据库中。（ii）对于每个基因，仅保留覆盖基因长度≥80％和同一性≥65％的得分最高的前10％的比对。（iii）为每个基因分配的比对分类法，其分类学等级的相似性阈值高于50％或更高的共识（门类> 65％，属> 85％，物种> 95％）。（iv）对70万个SPGC基因进行了分类。
如前所述[4,23]，我们的系统发育注释方法可确保唯一的分配并最大程度地减少歧义。门类和属类的假阳性率分别为0.77％和1.84％23。
注：现在有很多又好又快的物种注释工具，有Kraken2
功能注释KEGG和eggNOG
Functional annotation (KEGG and eggNOG)
我们使用BLASTP（v2.2.24，默认参数，除了默认参数之外）将整合基因目录中翻译的推定氨基酸序列与eggNOG（v3.0）和KEGG数据库中的蛋白质或域进行了比对（版本59.0，不包括来自动物或植物的基因）。-e 0.01 -b 100 -K 1 -FT）。KEGG注释是使用内部管道执行的，当得分最高的带注释的匹配包含至少一个超过60位的比对时，将每种蛋白质分配给KO。eggNOG注释是使用Smash社区（v1.6，find_best_hit.pl＆og_mapping.py，具有默认参数）执行的。
MetaHIT 2010，HMP 2012和IGC基因之间的比较
Comparison between MetaHIT 2010, HMP 2012 and IGC genes.
将来自IGC的990万个基因和来自MetaHIT 2010目录的330万个基因集中在一起，并使用CD-HIT鉴定出≥95％相同且≥90％重叠的冗余基因。对于共享的（重叠的）基因，将其长度进行比较，并且将大于10％的差异视为在IGC中明显更长或更短。使用相同的工作流程比较了来自IGC的990万个基因和来自更新的HMP 2012目录3的460万个基因。
注：blast太慢，目前可使用diamond或eggnog-mapper等工具进行功能注释。
将人类微生物测序公共数据比对到基因目录上
Aligning public human microbial sequencing data onto gene catalogs
使用BLASTN（v2.2.24），根据比对长度≥100 bp的标准，将来自18个美国双胞胎的Roche 454读长以及来自其13个日本人17的Sanger读长与MetaHIT 2010和IGC进行比对。可以与MetaHIT2010或IGC进行比对的读长比例通过两个重叠阈值（与基因目录比对的读长比例）进行过滤，分别为1％和80％（图2b）。
来自145个欧洲个体（其中130个在瑞典出生）的Illumina读长使用SOAP2与MetaHIT 2010和IGC进行了比对，同一性标准≥95％。
将人类公共微生物的转录组数据与基因目录进行比对
Aligning public human microbial metatranscriptomic data onto gene catalogs
通过直接从肠道微生物组构建的基因目录，我们能够使用SOAP2（≥95％相同性）将59个宏转录组测序样品的转录序列分配至基因目录（IGC），以鉴定表达的基因并检索其注释功能。SPGC是从1267个宏基因组中的511个肠道相关细菌或古细菌基因组或原始基因组编译而来的基因目录，用于与IGC进行比较，因为原始研究中使用的539种与人类相关的微生物参考基因组不详（截至2014年3月31日，使用的人类微生物组数据库现在包含2673个基因组，http：//www.hmpdacc.org/catalog/）。原始研究中使用的基于散列的软件SSAHA2具有非常宽松的比对标准（参数：“-best 1 -score 20 -solexa”），并且可能不支持对基因功能的明确鉴定。此外，就单位时间的比对碱基而言，SSAHA2比诸如SOAP2之类的短读比对剂要慢得多，并且太慢而无法处理我们的基因目录。
在将511个人类肠道相关参考基因组整合到我们的目录中时，涉及了一些非编码RNA。为了计算映射到蛋白质编码基因（编码序列；CDS）的读段的比例，我们在SPGC和IGC中消除了注释为非编码RNA的基因，关键字为“ RNA”，但没有原始基因组注释中的“ -ase”，“ enzyme”、“蛋白质”。根据该搜索标准，SPGC和IGC分别包含923和866个非编码RNA基因（rRNA，tRNA，SRP RNA等）。
相对基因丰度的计算
Computation of relative gene abundance
SOAP2使用相似度 ≥ 95％的标准，将每个样品的高质量读数与基因目录进行比对。如前所述，进行了基于序列的丰度分析。
属，KO和酶谱的构建
Construction of genus, KO and enzyme profiles.
对于属谱，我们使用了原始基因目录中每个基因的系统发育分配，并对来自同一属的基因的相对丰度求和，以得出该属的丰度。样品中每个属的相对丰度构成了该样品的属谱。KO表征使用相同的方法构建，酶的相对丰度由其相应的KO的相对丰度之和计算得出。
估计在基因边界的可比对读长的丢失
Estimating loss of mappable reads at gene boundaries
当针对基因目录进行定位时，一部分短读会在基因的边界区域丢失（补充图9a）。
丰度的损失率，损失率可以计算为，公式见原文
欧洲和中国同类人群使用的BMI标准
BMI criteria used for European and Chinese cohorts
许多报告表明，亚洲人的BMI标准低于欧洲人，并且亚洲人倾向于积累腹部脂肪并发展与肥胖相关的疾病，而没有整体肥胖。因此，我们使用较低的BMI临界值来定义中文的肥胖状况。对于中国人，我们使用的BMI值<21 kg / m2表示瘦，而≥25 kg / m2表示肥胖。对于丹麦人，我们使用的BMI值<25 kg / m2（瘦）和≥30 kg / m2（肥胖）。
生物多样性和丰富度分析：α多样性
Biodiversity and richness analysis: α-diversity.
根据基因，属或KO谱，我们计算了α多样性（样本内多样性），以使用Shannon指数估算样本的基因，属或KO多样性：
其中S是基因的数目，ai是基因i的相对丰度。高α多样性表示样品中存在高均匀度或许多类型的基因。
通过线性回归进行调整
Adjustment by linear regression
基于11个随机选择的样本生成了线性回归方程，这些样本的DNA使用BGI和MetaHIT方案提取了两次（补充图5c）。为了与丹麦人群进行比较，根据线性回归方程调整了中国人群的样本内多样性指数（n = 60），以消除不同DNA提取方案引入的可能偏倚。
读长数量标准化
Read downsizing
为了消除测序量波动的影响，我们对比对结果进行了采样，并将每个样本的比对数重采样至1100万。
稀疏曲线分析
Rarefaction curve analysis
为了评估我们中国或丹麦队列中的基因丰富性，我们生成了稀疏曲线。对于给定数量的单个样本，我们在队列中进行了有放回的随机抽样1000次，并估计了可以通过Chao2丰富度估算器从这些样本中识别出的基因总数。为了最大程度地减少错误识别，样品中仅存在具有≥1对可比对读长的基因。
肠道宏基因组的统计分析
Statistical analysis of the gut metagenome
为了确定宏基因组学概况与人群之间的关联，在概况中使用了双尾Wilcoxon秩和检验。如果至少一个队列的P值 < 0.01 并且发生频 > 10％，我们确定了属标记；如果至少一个队列的P <0.01且出现频率> 30％，我们确定了KO / 酶标记。
用于检测提取方法中偏差的统计方法与上述方法类似，但由于样本量仅为11，因此我们没有考虑出现频率。
估计错误发现率和统计能力
Estimating the false discovery rate and statistical power
没有使用顺序P值拒绝方法，我们应用了先前研究中建议的“ qvalue”方法来估计FDR。统计假设检验针对属谱和KO谱的大量特征进行。假设通过q值方法获得了FDR
模拟基因目录大小对抽样规模的依赖性
Simulation for the dependence of gene catalog size on sampling scale
模拟前准备了两个数据表。第一个是“基因概况”表，其中包含有关每个样品的每个基因的相对丰度的信息。该表是使用参考文献中的方法生成的。第二个是“基因组装”表，其中包含有关基因是否由于存在单个样本而被组装的信息。该表是从CD-HIT的聚类输出文件生成的，该文件将对应于同一聚类的基因（代表基因）追溯到原始样本。
通过随机抽样而不替换选定的样本进行模拟，样本大小为50到1,267个样本，每步50个样本。对于每个模拟集，根据“基因组装”表计算非冗余基因目录的估计大小，并通过“基因概况”表计算出现频率范围内的基因分布。对于每个样本大小，进行1000次模拟，并绘制平均值。
样本之间低频基因的一致性
Concordance of low-occurrence genes between samples
Sorenson指数，也称为Sørensen-Dice指数，用于测量HMP样本中所有出现频率（根据所有1,070个个体）的基因是否存在相似性。
Sørensen的原始公式以以下形式应用于存在/不存在数据。
Spearman秩相关系数用于测量HMP样本中所有出现频率的基因的丰度相似性。
个体内相似性基于在不同时间（平均间隔218 d）从相同的43个HMP个体获取的样本。个体间的相似性基于来自94个不同HMP个体的任意两个粪便样本（第一个时间点）。
不同发生频率的基因的来源属。
Source genera for genes of diverse occurrence frequencies.
每个IGC基因的出现频率四舍五入到最接近的整数，例如1％代表1％±0.5％。对于每个来源属，计算每个出现频率百分位数中的基因数量（补充表5）。补充图4e源自该表，其基因编号截止值为≥10。
数据访问编号
Accession codes
欧洲生物信息学研究所序列阅读档案：ERP004605（249个欧洲样品和11个中国样品的宏基因组测序数据）。
Reference
Li, J. et al. An integrated catalog of reference genes in the human gut microbiome. Nature Biotechnology 32, 834-841, doi:10.1038/nbt.2942 (2014).
猜你喜欢
10000+：菌群分析 宝宝与猫狗 梅毒狂想曲 提DNA发Nature Cell专刊 肠道指挥大脑
文献阅读 热心肠 SemanticScholar Geenmedical
16S功能预测 PICRUSt FAPROTAX Bugbase Tax4Fun
生物科普: 肠道细菌 人体上的生命 生命大跃进 细胞暗战 人体奥秘
写在后面
为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外5000+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。PI请明示身份，另有海内外微生物相关PI群供大佬合作交流。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍未解决群内讨论，问题不私聊，帮助同行。
学习16S扩增子、宏基因组科研思路和分析实战，关注“宏基因组”