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][18][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 .

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 METAPATH-WAYS. 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,  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 . . 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. 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 superpathways) 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. 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.
preprocessed For a sample with an identifier of o sample_ID>, this directory contains two files: (1) osample_ID>.fasta, which contains the renamed, quality-controlled sequences, and (2) osample_ID>.mapping.txt, which maps the original sequence names to the new names assigned by METAPATHWAYS. Sequences are renamed to o sample_ID>_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) osample_ID>.fna which contains nucleic acid sequences of all predicted ORFs, (2) osample_ID>.faa which contains amino acid sequences of all predicted ORFs, (3) osample_ID>.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) osample_ID>.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 osample_ID>_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.

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: R 2 = 0.887, P o0.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][18][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.
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 scilloinositol degradation pathway', has four matched compounds from the metabolomics dataset.   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.