A geographically-diverse collection of 418 human gut microbiome pathway genome databases.

Advances in high-throughput sequencing are reshaping how we perceive microbial communities inhabiting the human body, with implications for therapeutic interventions. Several large-scale datasets derived from hundreds of human microbiome samples sourced from multiple studies are now publicly available. However, idiosyncratic data processing methods between studies introduce systematic differences that confound comparative analyses. To overcome these challenges, we developed GutCyc, a compendium of environmental pathway genome databases (ePGDBs) constructed from 418 assembled human microbiome datasets using MetaPathways, enabling reproducible functional metagenomic annotation. We also generated metabolic network reconstructions for each metagenome using the Pathway Tools software, empowering researchers and clinicians interested in visualizing and interpreting metabolic pathways encoded by the human gut microbiome. For the first time, GutCyc provides consistent annotations and metabolic pathway predictions, making possible comparative community analyses between health and disease states in inflammatory bowel disease, Crohn's disease, and type 2 diabetes. GutCyc data products are searchable online, or may be downloaded and explored locally using MetaPathways and Pathway Tools.

Background & Summary

The myriad collections of microorganisms found on and in the human body are known as the human microbiome[1]. Changes in microbiome structure and function have been implicated in numerous disease states including inflammatory bowel disease, cancer, and even cardiovascular disease[2,3]. Increasingly, researchers are using high-throughput sequencing approaches to study the genes and genomes of microbiomes and characterize diversity and metabolic potential in relation to health and disease states[4], opening new opportunities for prevention and therapeutic intervention at the interface of microbial ecology, bioinformatics and medicine. The most densely colonized human habitat is the distal gut, inhabited by thousands of diverse microorganisms, as differentiated at the strain level. Despite providing essential ecosystem services, including nutritional provisioning, detoxification and immunological conditioning, the metabolic network driving matter and energy transformations by the distal gut microbiome remains largely unknown. Several large-scale metagenomic datasets (derived from hundreds of microbiome samples) from the Human Microbiome Project (HMP)[5], Beijing Genomics Institute (BGI)[6], and Metagenomes of the Human Intestinal Tract project (MetaHIT)[7] are now available on-line, creating an opportunity for large-scale metabolic network comparisons.

While the studies cited above provide the sequencing data, they do not provide the software environment used for generating their annotations. In contrast to these proprietary pipelines, over the past few years a number of metagenomic annotation pipelines available to third parties have emerged including IMG/M[8], Metagenome Rapid Annotation using Subsystem Technology (MG-RAST)[9], SmashCommunity[10] and HUMAnN[11]. Differing pipelines used to process sequence information between studies introduces biases based on idiosyncratic formatting, and alternative annotations or algorithmic methods. Specifically, support for metabolic pathway annotation varies significantly among pipelines due to differences in reference database selection with resulting impact on metabolic network comparisons. The most common metabolism reference database currently in use is Kyoto Encyclopedia of Genes and Genomes (KEGG)[12]. Although extant pipelines often provide links to KEGG module and pathways maps[12] (using KEGG ontology (KO) or pathway identifiers) that can be visualized with coverage or gene count information using programs like KEGG Atlas[13], they do so using often incompatible formats. Such mapping is limited because there is no simple way to query, manipulate, or visualize the underlying implicit metabolic model directly. Moreover, prediction using KEGG results in amalgamated pathways with limited taxonomic resolution, impeding enrichment and association studies[11].

In responding to the deficiencies of existing tools, we recently developed a modular annotation and analysis pipeline enabling reproducible research[14] called MetaPathways, that guides construction of environmental Pathway Genome Databases (ePGDBs) from environmental sequence information[15] using Pathway Tools[16] and MetaCyc[17-19]. Pathway Tools is a production-quality software environment developed at SRI International that supports metabolic inference and flux balance analysis based on the MetaCyc database of metabolic pathways and enzymes representing all domains of life. Unlike KEGG, MetaCyc emphasizes smaller, evolutionarily conserved or co-regulated units of metabolism and contains the largest collection (over 2,400) of experimentally validated metabolic pathways[20]. Navigable and extensively commented pathway descriptions, literature citations, and enzyme properties combined within an ePGDB provide a coherent structure for exploring and interpreting predicted metabolic networks from the human microbiome across multiple levels of biological information (DNA, RNA, protein and metabolites). Over 9,800 PGDBs have been developed by researchers around the world, and thus ePGDBs represent a data format for metabolic reconstructions that exhibit a potential for reusability in further studies.

Here we present GutCyc, a compendium of over 418 ePGDBs constructed from public shotgun metagenome datasets generated by the HMP[5], the MetaHIT inflammatory bowel disease study[7], and the BGI diabetes study[6]. Relevant pipeline modules are summarized in Fig. 1. GutCyc provides consistent taxonomic and functional annotations, facilitates large-scale and reproducible comparisons between ePGDBs, and directly links into robust software and database resources for exploring and interpreting metabolic networks. This metabolic network reconstruction provides a multidimensional view of the microbiome that invites discovery and collaboration[21].

Figure 1

The GutCyc pipeline diagram.

The MetaPathways pipeline consists of five modular stages including (1) Quality control (QC) and open reading frame (ORF) prediction (2) Functional and taxonomic annotation, (3) Analysis (4) ePGDB construction, and (5) Pathway export. Inputs and programs are depicted on the left with corresponding output directories and exported files on the right.

Methods

Metagenomic data sources

We collected 418 assembled human gut shotgun metagenomes from public repositories and Supplementary Materials sourced from the HMP (American healthy subjects, n=148)[5], a MetaHIT (European inflammatory bowel disease subjects, n=125)[22], and a BGI (Chinese type 2 diabetes subjects, n=145) study[6]. See Supplementary Tables 1 and 2 for a detailed listing of accession numbers and file descriptors.

Data processing

Microbiome project sample metadata were manually curated to ensure compatibility with MetaPathways. ePGDBs were created for each sample by running the MetaPathways 2.5 pipeline and the Pathway Tools version 17.5, using the assembled metagenomes described above. The pipeline consists of five modular steps, including (1) quality control and ORF prediction, (2) homology-based functional and taxonomic annotation, (3) analyses consisting of tRNA and lowest common ancestor (LCA)[23] identification, (4) construction of ePGDBs using Pathway Tools and, finally, (5) pathway export[24,25] (see Fig. 1). The following paragraphs describe the individual processing steps required to construct an ePGDB for each sample, starting with assembled contigs in FASTA format.

Quality control

Contigs from each sample were collected from their respective repositories and curated locally. The MetaPathways pipeline performs a number of quality control steps. First, each contig was checked for the presence of ambiguous base pairs and homopolymer runs, splitting contigs into smaller sequences by removing such problematic regions. Next, the contigs were screened for duplicates. Finally, a length cutoff of 180 base pairs was applied to the remaining sequences to ensure that very short sequences were removed from downstream processing steps[26].

ORF prediction

Sequences passing quality control were scanned for ORFs using MetaProdigal[27], a robust ORF prediction tool for microbial metagenomes considered to be among the most accurate ORF predictors[28]. Resulting ORF sequences were translated to amino acid sequences using NCBI genetic code Table 11 for bacteria, archaea, and plant plastids[29]. Translated amino acid sequences shorter than 30 amino acids were removed as these sequences approached the so-called functional homology search ‘twilight zone’, where it becomes difficult to detect true homologs[30].

Functional annotation

The quality-controlled amino acid sequences were queried against a panel of functionally-annotated protein reference databases used in the validation of the MetaPathways pipeline[31]: KEGG[12] (downloaded 2011-06-18), COG[32] (downloaded 2013-12-27), MetaCyc[19] (downloaded 2011-07-03), RefSeq[33] (downloaded 2014-01-18), and SEED[34] (downloaded 2014-01-30). Protein sequence similarity searches were performed using the program FAST[35] with standard alignment result cutoffs (E-value less than 1×10−5, bit-score greater than 20, and alignment length greater than 40 amino acids[30]; and blast-score ratio (BSR) greater than 0.4 (ref. 36)). The choice of parameter thresholds were selected to maximize annotation accuracy, and were guided based on parameter choices used in previous studies[31,37,38].

Taxonomic annotation

Quality-controlled contigs were also searched against the SILVA[39] (version 115) and GreengeneS[40] (downloaded 2012-11-06) ribosomal RNA (rRNA) gene databases using BLAST version 2.2.25, with the same post-alignment thresholds applied as was previously described for the functional annotation. BLAST was applied for rRNA gene annotation because it has greater sensitivity for nucleotide-nucleotide searches than FAST, and the smaller reference database sizes make the relatively larger computational requirement justifiable.

Additionally, predicted ORFs were taxonomically annotated using the LCA algorithm[23] for taxonomic binning. In brief, the LCA in the NCBI Taxonomy Database[33] was selected based on the previously calculated FAST hits from the RefSeq database. This effectively sums the number of FAST hits at the lowest shared position of the NCBI Taxonomy Database. The RefSeq taxonomic names often contain multiple synonyms or alternative spellings. Therefore, names that conform to the NCBI Taxonomy Database were selected in preference over unknown synonyms.

tRNA scan

MetaPathways uses tRNAscan-SE version 1.4 (ref. 41) to identify relevant tRNAs from quality-controlled sequences. Resulting tRNA identifications are appended as additional functional annotations.

ePGDB creation

Functional annotations were parsed and separated into three files that serve as inputs to Pathway Tools, namely: (1) an annotation file containing gene product information (0.pf), (2) a catalog of contigs and scaffolds (genetic-elements.dat), and (3) a PGDB parameters file (organism-params.dat). The PathoLogic module[42,43] in the Pathway Tools software, was used to build the ePGDB and predict the presence of metabolic pathways based on functional annotations. Following ePGDB construction, the base pathways (i.e., pathways that are not contained within super-pathways) were extracted from ePGDBs to generate a summary table of predicted metabolic pathways for each sample.

Accessibility and flexibility

MetaPathways 2.5 generates data in a consistent file and directory structure. The output for each sample is contained within a single directory, which in turn is organized into sub-directories containing relevant files (see Fig. 1). The MetaPathways 2.5 graphical user interface (GUI) enables interactive exploration, visualization, and searches of individual sample results along with comparative queries of multiple samples, via a custom knowledge engine data structure. Input and output files are available for download from the GutCyc website (www.gutcyc.org) and may be readily explored in the MetaPathways GUI or Pathway Tools on Linux, Mac OS X and Windows machines.

Computational environment

Computational processing was performed using a local cluster of machines in the Hallam laboratory and on WestGrid's Bugaboo cluster part of Compute Canada's national platform of Advanced Research Computing resources. https://www.westgrid.ca/support/systems/bugaboo. The Hallam lab computers have a configuration profile of 2×2.4 GHz Quad-Core Intel Xeon processors with 64 GB 1,066 MHz DDR3 RAM. The Bugaboo cluster provides 4,584 cores with 2 GB of RAM per core on average. The average sample took 7–8 h to process on a single thread, and the span of the processing required to generate the GutCyc Collection was 135 days.

Software availability

MetaPathways 2.5, including integrated third party software, is available on GitHub (github.com/hallamlab/metapathways2), licensed under the GNU General Public License, version 3), along with a companion tutorial (github.com/hallamlab/mp_tutorial/) released under the Creative Commons Attribution License (allows reuse, distribution, and reproduction given proper citation). Pathway Tools is available under a free license for academic use, and may be commercially licensed (www.biocyc.org/download-bundle.shtml). MetaPathways outputs were processed using Pathway Tools version 17.5 under default settings except for disabling of the PathoLogic taxonomic pruning step (i.e., -no-taxonomic-pruning) as was described previously[31], and an additional refinement step of running the Transport Inference Parser[44] to predict transport reactions (i.e., -tip). FAST is freely available under the GNU General Public License, version 3 on our GitHub page (github.com/hallamlab/FAST).

Data Records

A list of each sample, its provenance, location and relevant data processing steps can be found in Supplementary Table 1. All records are available at the GutCyc project website (www.gutcyc.org), and at Figshare as described in (Data Citation 1). Each sample’s data records are contained within a single directory. Within this directory, sub-directories and files are located as depicted in Fig. 1. A summary of the data present in the GutCyc Collection is presented in Table 1. A full set of summary data for each ePGDB may be found in Supplementary Table 2.

Table 1

Summary statistics for the GutCyc Collection across 418 samples.

	Min	1st quartile	Median	3rd quartile	Max
The statistics for the number of bases processed is in units of Megabases. ‘Func. Annots.’ are functional annotations. ‘Trans. Reactions’ are transport reactions. ‘Compounds’ are small molecule metabolites. ‘Base Pathways’ include all pathways except complex pathways known as Super-Pathways.
Bases	0.98	54.75	81.35	113.75	370.51
Contigs	2,506	27,788	47,486.5	76,275.75	399,331
ORFs	2,448	61,703.5	95.531	139,690	550,312
Func. Annots.	2,176	57,102.25	86,054.5	123,747.25	425,033
Reactions	1,635	2,385.5	3,438	3,667.75	4,881
Trans. Reactions	12	26	31	34	46
Compounds	1,052	1,678	2,008.5	2,119.5	2,676
Base Pathways	257	350	616	654	832

preprocessed

For a sample with an identifier of , this directory contains two files: (1) .fasta, which contains the renamed, quality-controlled sequences, and (2) .mapping.txt, which maps the original sequence names to the new names assigned by MetaPathways. Sequences are renamed to _X where X is the zero-indexed contig number pertaining to the order in which the contig appears in the input file (e.g., a contig identified as DLF001_27 is interpreted as the 28th contig listed in the FASTA file for sample DLF001 ‘s assembly).

orf_prediction

This directory contains four files, (1) .fna which contains nucleic acid sequences of all predicted ORFs, (2) .faa which contains amino acid sequences of all predicted ORFs, (3) .qced.faa which contains amino acid sequences of all predicted ORFs meeting user defined quality thresholds (in this study, a minimum length of 60 amino acids), and (4) .gff, a general feature format (GFF) file[45] containing all quality-controlled sequences and information about the strand (− or +) on which the ORF was predicted. ORFs are named _X_Y, where X is the contig number pertaining to the order in which the contig appears and Y represents the order in which the ORFs were predicted on the contig.

results

This directory contains four sub-directories: (1) annotation_table, (2) rRNA, (3) tRNA, and (4) pgdb. The annotation_table sub-directory contains (1) statistics files for each functional database used to annotate the ORFs (._stats_.txt), (2) .functional_and_taxonomic_table.txt detailing the length, location, strand and annotation (functional and taxonomic) of each ORF, and (3) a file listing all ORFs and their functional annotations (.ORF_annotation_table.txt). The prokaryotic small subunit ribosomal RNA (SSU or 16S rRNA) gene is a standard marker gene used for measuring taxonomic diversity[46]. The rRNA sub-directory contains files detailing statistics for each taxonomic database used to annotate the ORFs (named as ..rRNA.stats.txt). The tRNA sub-directory contains (1) .trna.stats.txt, detailing the type, anticodon, location and strand of each predicted tRNA and (2) .tRNA.fasta containing all predicted tRNA sequences. The pgdb sub-directory contains a .pwy.txt file describing metabolic pathways predicted in the ePGDB, specifically, each predicted pathway, the ORF identities involved in each pathway, the enzyme abundance, and the pathway coverage in a tabular format navigable via the MetaPathways GUI.

genbank

This directory contains a file named .annotated.gff, a GFF file containing all quality-controlled sequences with their annotations.

ptools

This directory contains the three files necessary to build a ePGDB using Pathway Tools: (1) genetic-elements.dat, (2) organism-params.dat, and (3) 0.pf which contains all functional annotations to be processed by Pathway Tools. A sub-directory called flat-files contains flat files describing database objects such as compounds, reactions, pathways (each of which is described in more detail in[47]) for individual ePGDBs.

run_statistics

This directory contains three files: (1) .run.stats, the parameters used to process the sample; (2) .nuc.stats, the number and length of nucleic acid sequences before and after quality control filtering; and (3) .amino.stats, the number and length of amino acid sequences before and after quality control filtering.

Technical Validation

GutCyc was derived from third-party sequence data from three publicly-available human gut microbiome sampling projects with metagenomic assemblies, with published details on their own technical validation steps: the HMP[5], a MetaHIT study[22], and a BGI study[6]. The technical validation of third-party software used in MetaPathways may be found in the corresponding publications for MetaProdigal[27], BLAST[48], and tRNAscan-SE[41]. GutCyc functional sequence similarity was computed using FAST, an aligner based on a version of LAST[49], with multi-threading performance improvements and new support for generating BLAST-like E-values, with significant correlation with BLAST output (correlation of the log(E-value) outputs of BLAST and LAST: R2=0.887, P<0.01)[24]. The protocols undertaken in the MetaCyc project for the ongoing manual curation of new metabolic pathways, and its subsequent implications for accurate pathway prediciton, may be found in the following MetaCyc publications[17-19,50].

Validation of the overall MetaPathways pipeline may be found in previously published reports[31,51] with specific emphasis on how changes in taxonomic pruning, read length and metagenomic assembly coverage impact the accuracy and sensitivity of pathway recovery. In brief, pathway prediction is affected by taxonomic distance, sequence coverage and sample diversity, nearing an asymptote of maximum accuracy for metagenomes with increasing coverage. Additionally, like any alignment-based analysis, annotation quality is a function of both the level of errors in the input sequence data and the selection of reference databases. Summary data generated for each ePGDB as presented in Supplementary Table 2 was reviewed to detect samples with unusual statistics, such as a lack of 16S gene annotations. The metabolic reconstruction pathways were computationally predicted using the Pathway Tools PathoLogic module[52], which has an accuracy of 91% as evaluated using organism pathway databases with high levels of manual curation[43]). The performance of the Pathway Tools PathoLogic module has also been evaluated using datasets with different complexity and coding potential, including simulated metagenomes, a symbiotic system, and the Hawaii Ocean Time-series[31]. The authors provide detailed information about the effects of read length, coverage and sample diversity on pathway recovery but found that performance specificity was high (>85%) using all three datasets. The authors also provide a list of `prediction hazards’ such as the identification of dissimilatory nitrate reduction pathways of which the user should be aware and conclude that despite being imperfect, Pathway Tools provides a powerful means with which to predict metabolic interactions[31].

Usage Notes

Once a set of data such as GutCyc Collection has been crafted into a format that is both comprehensible to domain experts, and interpretable by machines, there are myriads of uses that can be explored. For example, comparing ePGDBs with sets of microbial PGDBs from the same environment can aid in identifying ‘distributed pathways’ present in the metagenome metabolic reconstruction, but absent from any individual genomic metabolic reconstruction[31]. Annotations from each of the protein reference databases can also be explored individually using MetaPathways. In addition, a file for each sample, located at sample/results/annotation_table/sample.2.txt provides a detailed overview of the annotation for each ORF in each database as well as score reflecting the confidence of the alignment and annotation. The predicted transport proteins can be used to predict trophism patterns within a community. Furthermore, the Pathway Tools software allows for sophisticated comparative analyses between ePGDBs, at the level of compounds, reactions, enzymes, and pathways[53]. The MetaFlux[54] module of Pathway Tools for performing flux balance analysis (FBA)[55] can be used with GutCyc ePGDBs to generate quantitative simulations of microbiome growth and pathway flux. A set of microbiome metabolic models also facilitates the exploration of the impact of xenobiotics[56], and provides a computational substrate for systems biology approaches to engineering the gut microbiome[57]. Figure 2 demonstrates the user interface for MetaPathways and Pathway Tools, along with example data analysis use cases.

Figure 2

GutCyc ePGDB use cases.

In the upper left and upper right insets, a GutCyc ePGDB is opened in MetaPathwaysIn the upper left, we display the Pipeline Execution step, and the Process Monitor interfaces. In the upper right, we display the Summary Table (with exportable sample statistics), and the Pathway Table (with exportable pathway abundances) interfaces. In the lower four inset images, a GutCyc ePGDB is opened in Pathway Tools. Clockwise from the upper left, we display the ePGDB summary statistics, interactive metabolic network visualization, the Pathway View, and the biochemical Reaction View.

In this section we motivate further two specific use cases for GutCyc. In the first case, we demonstrate how to use a GutCyc ePGDB to determine the metabolic path between two small molecules of interest. In the second case, we use GutCyc to visualize different levels of biological information, e.g., metabolomics data, in the context of a microbiome metabolic network.

Optimal metabolite tracing

The Pathway Tools software provides advanced biochemical querying capabilities for both PGDBs and ePGDBs. One such capability is energy-optimal metabolite tracing. To summarize, given both a starting and a terminal/target compound within an ePGDB, Pathway Tools is able to compute the shortest and most energetically-favorable route through the metabolic reaction network. While there is no guarantee that, in a complex milieu such as the gut microbiome, the syntrophic flux will necessarily follow a short and minimal energy path, these criteria allow us to narrow down a multiplicity of possible paths to a single parsimonious candidate path.

In a study by Koeth et al.[58], they demonstrated a causal connection between the intestinal gut microbiota’s metabolism of red meat and the promotion of atherosclerosis. In brief, the gut microbiome is capable of transforming excess L-carnitine into trimethylamine (TMA), which is further processed by the liver to form the cardiovascular disease-associated metabolite trimethylamine N-oxide (TMAO). Using this biotransformation as a motivating case, we queried an arbitrarily selected ePGDB from the GutCyc Collection, SRS015217Cyc, for the biochemical reaction path from L-carnitine to TMA, which is not provided explicitly by Koeth et al.[58] Utilizing the Pathway Tools Metabolic Route Search feature, we found an optimal path between L-carnitine to TMA for this ePGDB, using the MetaCyc carnitine degradation II pathway (PWY-3,602, expected in Proteobacteria) along with a betaine reductase reaction (EC 1.21.4.4; found in Clostridium sticklandii and Eubacterium acidaminophilum, both species affiliated with the order Clostridiales), minimizing the number of enzymes involved and chemical bond rearrangements. Pathway Tools found the optimal path in seconds.

The metabolic route identified may also help generate mechanistic hypotheses from microbiome study observations. L-carnitine and glycine betaine have known transporter families that facilitate their movement across the cell membrane[59], as do TMA and TMAO[60], and thus the metabolic route in this ePGDB may be a distributed pathway[31]. This demonstrates the power of ePGDBs in computing connections between nutritional or pharmaceutical inputs (such as L-carnitine) to identify potential interactions with known disease biomarkers (as TMAO is to cardiovascular disease).

High-throughput data visualization

Another capability of Pathway Tools is to visualize the results of high-throughput experiments mapped onto the Cellular, Genome, and Regulation Overviews, or as ‘Omics Pop-Ups’ when viewing a particular pathway[61]. Specifically, Pathway Tools provides support for the analysis of mass spectrometry data, by automatically mapping a list of monoisotopic masses to matching entries in MetaCyc, or in specific ePGDBs[62]. As a demonstration of this capability, we analyzed mass-spectrometry data from a metabolomic study of humanized mice microbiomes[63]. The dataset contained 867 unique masses, of which 453 masses were identified using MetaCyc by performing standard adduct monoisotopic mass manipulations[64], followed by monoisotopic mass search using Pathway Tools. We mapped the identified compounds on the Cellular Overview[65] of an arbitrarily-selected ePGDB from the GutCyc Collection for illustrative purposes, as seen in Fig. 3. This facilitates turning a massive table of data into a more intuitive construct based on the community metabolic interaction network, enabling more efficient pattern matching. For example, using the enrichment analysis tools in Pathway Tools[62], we identified the pathway class of ‘Secondary Metabolites Degradation’ as enriched for identified compounds (P=2.0×10−2, Fisher Exact Test with Benjamini-Hochberg multiple testing correction). By visually inspecting the pathways in the class, we can see that pathway P562-PWY, ‘myo-, chiro-, and scillo-inositol degradation pathway’, has four matched compounds from the metabolomics dataset.

Figure 3

The Cellular Overview for the SRS056259Cyc ePGDB, at three different zoom levels.

Compounds are highlighted in red if identified from a mass spectrometry analysis of the gut microbiome[63], and otherwise appear in grey. Reactions with enzyme data in SRS056259Cyc are drawn in blue. The top left inset shows a fraction of the full metabolic map. The middle inset shows a zoom-in of the ‘Secondary Metabolite Degradation’ pathway class. Bottom right inset shows zoom-in on Pathway P562-PWY, ‘myo-, chiro-, and scillo-inositol degradation pathway’, showing four mass-spectrometry identified compounds in red.

Additional Information

How to cite this article: Hahn, A. S. et al. A geographically-diverse collection of 418 human gut microbiome pathway genome databases. Sci. Data 4:170035 doi: 10.1038/sdata.2017.35 (2017).

Publisher’s note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.