Given the ever-increasing number of microbial genomes from microbiome studies, we developed METABOLIC to enable the metabolic pathway analysis and the visualization of biogeochemical cycles and community-scale functional networks. METABOLIC has an improved methodology to get fast, accurate, and robust annotation results, and it integrates a variety of visualization functions for better interpreting the community-level functional interactions and microbial contributions. While METABOLIC relies on microbial genomes and metagenomic reads for underpinning its analyses for community-level functional interactions, it can easily integrate transcriptomic datasets to provide an activity-based measure of community networks. The scalable capacity, wide utility, and compatibility for analyzing datasets from various environments make it a well-tailored tool for metabolic profiling of large sets of genomes. In the following sections, the microbial community consisting of 98 MAGs from a deep-sea hydrothermal vent was used as the input dataset if not mentioned otherwise.
Workflow to determine the presence of metabolic pathways
METABOLIC is written in Perl and R and is expected to run on Unix, Linux, or macOS. The prerequisites are described on METABOLIC’s GitHub wiki pages (https://github.com/AnantharamanLab/METABOLIC/wiki). The input folder requires microbial genome sequences in FASTA format and an optional set of genomic/metagenomic reads which were used to reconstruct those genomes (Figure 1). The annotated proteins from input genomic sequences are queried against HMM databases (KEGG KOfam, Pfam, TIGRfam, and custom HMMs) using hmmsearch implemented within HMMER [41] which applies methods to detect remote homologs as sensitively and efficiently as possible. After the hmmsearch step, METABOLIC subsequently validates the primary outputs by a motif-checking step for a subset of protein families; only those protein hits which successfully pass this step are regarded as positive hits.
METABOLIC relies on matches to the above databases to infer the presence of specific metabolic pathways in microbial genomes. Individual KEGG annotations are inferred in the context of KEGG modules for a better interpretation of metabolic pathways. A KEGG module is comprised of multiple steps with each step representing a distinct metabolic function. We parsed the KEGG module database [61] to link the existing relationship of KO identifiers to KEGG module identifiers to project our KEGG annotation result into the interactive network which was constructed by individual building blocks – modules – for better representation of metabolic blueprints of input genomes. In most cases, we used KOfam HMM profiles for KEGG module assignments. For a specific set of important metabolic marker proteins and commonly misannotated proteins, we also applied the TIGRfam/Pfam/custom HMM profiles and motif-validation steps. The software has customizable settings for increasing or decreasing the priority of specific databases, primarily meant to increase annotation confidence by preferentially using custom HMM databases over KEGG KOfam when both targeting the same set of proteins.
Since individual genomes from metagenomes and single-cell genomes can often have incomplete metabolic pathways due to their low completeness compared to isolate genomes, we provide an option to determine the completeness of a metabolic pathway (or a module here). A user-defined cutoff is used to set the threshold of completeness for a given module to be assigned as present (the default cutoff is the presence of 75% of metabolic steps/genes within a given module), which is then used to produce a KEGG module presence/absence table. All modules exceeding the cutoff value are determined to be present. Meanwhile, the presence/absence information for each module step is also summarized in an overall output table to facilitate further detailed investigations.
Outputs consist of six different results that are reported in an Excel spreadsheet (Additional file 1: Figure S1). These contain details of protein hits (Additional file 1: Figure S1A) which include both presence/absence and protein names, presence/absence of functional traits (Additional file 1: Figure S1B), presence/absence of KEGG modules (Additional file 1: Figure S1C), presence/absence of KEGG module steps (Additional file 1: Figure S1D), carbohydrate-active enzyme (CAZyme) hits (Additional file 1: Figure S1E) and peptidase/inhibitor hits (Additional file 1: Figure S1F). For each HMM profile, the protein hits from all input genomes can be used to construct phylogenetic trees or further be combined with reference protein collections for detailed evolutionary analyses.
Quantitative visualization of biogeochemical cycles and sequential reactions
After METABOLIC generates protein and pathway annotation results, the software further identifies and highlights specific pathways of importance in microbiomes associated with energy metabolism and biogeochemistry. To visualize pathways of biogeochemical importance, it generates schematic profiles for nitrogen, carbon, sulfur, and other elemental cycles for each genome. The set of genomes used as input is considered the “community”, and each genome within is considered an “organism”. A summary schematic diagram at the community level integrates results from all individual genomes within a given dataset (Figure 2) and includes computed abundances for each step in a biogeochemical cycle if the genomic/metagenomic read datasets are provided. The genome number labeled in the figure indicates the number/quantity of genomes that contain the specific gene components of a biogeochemical cycling step (Figure 2) [2]. In other words, it represents the number of organisms within a given community inferred to be able to perform a given metabolic or biogeochemical transformation. The abundance percentage indicates the relative abundance of microbial genomes that contain the specific gene components of a biogeochemical cycling step among all microbial genomes in a given community (Figure 2) [2].
Microorganisms in nature often do not encode pathways for the complete transformation of compounds. For example, microorganisms possess partial pathways for denitrification that can release intermediate compounds like nitrite, nitric oxide, and nitrous oxide in lieu of nitrogen gas which is produced by complete denitrification [62]. A greater energy yield could be achieved if one microorganism conducts all steps associated with a pathway (such as denitrification) [2] since it could fully use all available energy from the reaction. However, in reality, few organisms in microbial communities carry out multiple steps in complex pathways; organisms commonly rely on other members of microbial communities to conduct sequential reactions in pathways [2, 63, 64]. Thus, to study this metabolic scenario in microbial communities, METABOLIC summarizes and enables visualization of the genome number and coverage (relative abundance) of microorganisms that are putatively involved in the sequential transformation of both important inorganic and organic compounds (Figure 3). This provides a quantitative calculation of microbial interactions and connections using shared metabolites associated with inorganic and organic transformations. Additionally, it shows the intuitive pattern of quantity and abundance of microorganisms that are able to conduct partial or all steps for a given pathway, which potentially reflects the degree of resilience of a microbial community.
Calculation and visualization of functional networks, metabolic weight scores (MW-scores), and microbial contribution to metabolic reactions
Given the microbial pathway abundance information generated by METABOLIC, we identified co-existing metabolisms in microbial genomes as a measure of connections between different metabolic functions and biogeochemical steps. In the context of biogeochemistry, this approach allows the evaluation of relatedness among biogeochemical steps and the connection contribution by microorganisms. This is enabled at the resolution of individual microbial groups based on the phylogenetic classification (Figure 4) assigned by GTDB-Tk [50]. As an example, we have demonstrated this approach on a microbial community inhabiting deep-sea hydrothermal vents. We divided the microbial community of deep-sea hydrothermal vents into 18 phylum-level groups (except for Proteobacteria which were divided into their subordinate classes). The functional network diagrams were depicted at the resolution of both individual phyla and the entire community level (Additional file 10: Dataset S3). Figure 4 demonstrates metabolic connections that were represented with individual metabolic/biogeochemical cycling steps depicted as nodes, and the connections between two given nodes depicted as edges. The size of a given node is proportional to the gene coverage associated with the metabolic/biogeochemical cycling step. The thickness of a given edge was depicted based on the average of gene coverage values of two biogeochemical cycling steps (the connected nodes). More edges connecting two nodes represent more connections between these two steps. The color of the edge corresponds to the taxonomic group. At the whole community level, more abundant microbial groups were more represented in the diagram (Figure 4). Overall, METABOLIC provides a comprehensive approach to construct and visualize functional networks associated with important pathways of energy metabolism and biogeochemical cycles in microbial communities and ecosystems.
To address the lack of quantitative and reproducible measures to represent potential metabolic interactions in microbial communities, we developed a new metric that we termed MW-score (metabolic weight scores) (Equations 1 and 2). MW-scores quantitatively measure “function weights” within a microbial community as reflected by the metabolic profile and gene coverage. As metabolic potential for the whole community was profiled into individual functions that either mediated specific pathways or transformed certain substrates into products, a function weight that reflects the abundance fraction for each function can be used to represent the overall metabolic potential of the community. MW-scores resolved the functional capacity and abundance in the co-sharing functional networks as studied and visualized in the above section. More frequently shared functions and their higher abundances lead to higher MW-scores, which quantitatively reflects the function weights in functional networks (Figure 5). MW-score reflects the same functional networking pattern as the above description on the edges (networking lines) connecting the nodes (metabolic steps) that – more edges connecting two nodes indicates two steps are more co-shared, thicker edges indicate higher gene abundance for the metabolic steps. The MW-scores can integratively represent these two networking patterns and serve as metrics to measure these function weights. At the same time, we also calculated each microbial group’s (phylum in this case) contribution to the MW-score of a specific function within the community (Figure 5). A higher microbial group contribution percentage value indicates that one function is more represented by the microbial group (for both gene presence and abundance) in the functional networks. MW-scores provide a quantitative measure of comparing function weights and microbial group contributions within functional networks.
To understand the contributions of microbial groups associated with specific metabolic and biogeochemical transformations, we developed an approach to visualize the connections among specific taxonomic groups, metabolic reactions, and entire biogeochemical cycles such as carbon, nitrogen, and sulfur cycles. Our approach involves the use of Sankey diagrams (also called ‘Alluvial’ plots) to represent the fractions of metabolic functions that are contributed by various microbial groups in a given community (Figure 6). It allows visualization of metabolic reactions as the link between microbial contributors clustered as taxonomic groups and biogeochemical cycles at a community level (Figure 6 and Additional file 10: Dataset S3). The function fraction was calculated by accumulating the genome coverage values of genomes from a specific microbial group that possesses a given functional trait. The width of curved lines from a specific microbial group to a given functional trait indicates their corresponding proportional contribution to a specific metabolism (Figure 6). Alternatively, the genomic/metagenomic datasets which are used in constructing the above two diagrams: functional network diagram (Figure 4) and metabolic Sankey diagram (Figure 6), can be replaced by transcriptomic/metatranscriptomic datasets, and correspondingly, the gene coverage values will be replaced by gene expression values, and therefore, diagrams will represent the transcriptional activity patterns of functional network and microbial contribution to metabolic reactions (Additional file 2, 3, 4, and 5: Figure S2, S3, S4, and S5).
To demonstrate this part of the workflow in reality, the microbial community consisting of 98 MAGs from a deep-sea hydrothermal vent was used as a test dataset. After running the bioinformatic analyses described above, resulting tables and diagrams were compiled and visualized accordingly (Figure 4, 5, 6 and Additional file 10: Dataset S3). Results for functional networks and MW-scores of the deep-sea hydrothermal vent environment indicate that the microbial community depends on mixotrophy and sulfur oxidation for energy conservation and involves arsenate reduction potentially responsible for detoxification/arsenate resistance [65]. MW-scores indicate that amino acid utilization, complex carbon degradation, acetate oxidation, and fermentation are the major heterotrophic metabolisms for this environment; CO2-fixation and sulfur oxidation also occupy a considerable functional fraction, which indicates heterotrophy and autotrophy both contribute to energy conservation (Figure 5). As represented by both MW-scores and metabolic Sankey diagram, Gammaproteobacteria are the most numerically abundant group in the community and they occupy significant functional fractions amongst both heterotrophic and autotrophic metabolisms (MW-score contribution ranging from 59-100%) (Figure 5, 6), which is consistent with previous findings in the Guaymas Basin hydrothermal environment [52, 66]. Meanwhile, MW-scores also explicitly reflect the involvement of other minor electron donors in energy conservation which are mainly contributed by Gammaproteobacteria, such as hydrogen and methane (Figure 5). This is also consistent with previous findings [52, 66] and indicates the accuracy and sensitivity of MW-scores to reflect metabolic potentials.
METABOLIC performance demonstration
To test METABOLIC’s performance on speed, we applied the software (METABOLIC-C mode) to analyze the metagenomic dataset which includes 98 MAGs from a deep-sea hydrothermal vent, and two sets of metagenomic reads (that are subsets of original reads with 10 million reads for each pair comprising ~10% of the total reads). The total running time was ~3 hours using 40 CPU threads in a Linux version 4.15.0-48-generic server (Ubuntu v5.4.0). The most compute-demanding step is hmmsearch, which took ~45 mins. When tested on another dataset comprising ~3600 microbial genomes (data not shown), METABOLIC could complete hmmsearch in ~5 hours by using 40 CPU threads, indicating its scalable capability on analyzing thousands of genomes.
In order to test the accuracy of the results predicted by METABOLIC, we picked 15 bacterial and archaeal genomes from Chloroflexi, Thaumarchaeota, and Crenarchaeota which are reported to have 3 hydroxypropionate cycle (3HP) and/or 3-hydroxypropionate/4-hydroxybutyrate cycle (3HP/4HB) for carbon fixation. METABOLIC predicted results in line with annotations from the KEGG genome database which can be visualized in KEGG Mapper (Table 1). Our predictions are also in accord with biochemical evidence of the existence of corresponding carbon fixation pathways in each microbial group: 1) 3 out of 5 Chloroflexi genomes are predicted by both METABOLIC and KEGG to possess the 3HP pathway and none of all these Chloroflexi genomes are predicted to possess the 3HP/4HB pathway. This is consistent with current reports based on biochemical and molecular experiments that only organisms from the phylum Chloroflexi are known to possess the 3HP pathway [67] (Table 1). 2) All 5 Thaumarchaeota genomes and 2 out of 5 Crenarchaeota genomes are predicted by both METABOLIC and KEGG to possess the 3HP/4HB pathway and none of these Thaumarchaeota and Crenarchaeota genomes are predicted to possess the 3HP pathway. This is consistent with current reports that only the 3HP/4HB pathway could be detected in Crenarchaeota and Thaumarchaeota [68, 69] (Table 1). We have also applied METABOLIC on a large well-studied dataset comprising 2545 metagenome-assembled genomes from terrestrial subsurface sediments and groundwater [2]. The annotation results of METABOLIC are consistent with previously described reports (Additional file 6, 10: Figure S6, Dataset S3). These results suggest that METABOLIC can provide accurate annotations and perform well as a functional predictor for microbial genomes and communities.
Table 1
The carbon fixation metabolic traits of 15 tested bacterial and archaeal genomes predicted by both METABOLIC and KEGG genome database
| METABOLIC result | KEGG genome pathway |
Carbon fixation | Carbon fixation |
Accession ID | Organism | KEGG Organism Code | Group | 3HP cycle | 3HP/4HB cycle | 3HP cycle | 3HP/4HB cycle |
GCA_000011905.1 | Dehalococcoides mccartyi 195 | det | Chloroflexi | Absent | Absent | Absent | Absent |
GCA_000017805.1 | Roseiflexus castenholzii DSM 13941 | rca | Chloroflexi | Present | Absent | Present | Absent |
GCA_000018865.1 | Chloroflexus aurantiacus J-10-fl | cau | Chloroflexi | Present | Absent | Present | Absent |
GCA_000021685.1 | Thermomicrobium roseum DSM 5159 | tro | Chloroflexi | Absent | Absent | Absent | Absent |
GCA_000021945.1 | Chloroflexus aggregans DSM 9485 | cag | Chloroflexi | Present | Absent | Present | Absent |
GCA_000299395.1 | Nitrosopumilus sediminis AR2 | nir | Thaumarchaeota | Absent | Present | Absent | Present |
GCA_000698785.1 | Nitrososphaera viennensis EN76 | nvn | Thaumarchaeota | Absent | Present | Absent | Present |
GCA_000875775.1 | Nitrosopumilus piranensis D3C | nid | Thaumarchaeota | Absent | Present | Absent | Present |
GCA_000812185.1 | Nitrosopelagicus brevis CN25 | nbv | Thaumarchaeota | Absent | Present | Absent | Present |
GCA_900696045.1 | Nitrosocosmicus franklandus NFRAN1 | nfn | Thaumarchaeota | Absent | Present | Absent | Present |
GCA_000015145.1 | Hyperthermus butylicus DSM 5456 | hbu | Crenarchaeota | Absent | Absent | Absent | Absent |
GCA_000017945.1 | Caldisphaera lagunensis DSM 15908 | clg | Crenarchaeota | Absent | Present | Absent | Present |
GCA_000148385.1 | Vulcanisaeta distributa DSM 14429 | vdi | Crenarchaeota | Absent | Absent | Absent | Absent |
GCA_000193375.1 | Thermoproteus uzoniensis 768-20 | tuz | Crenarchaeota | Absent | Present | Absent | Present |
GCA_003431325.1 | Acidilobus sp. 7A | acia | Crenarchaeota | Absent | Absent | Absent | Absent |
Currently, several software packages and online servers are available for genome annotation and metabolic profiling. Comparing to other software/online servers including GhostKOALA [70], BlastKOALA [70], KAAS [71], RAST/SEED [33], and eggNOG-mapper [72], METABOLIC is unique in its ability to integrate multi-omic information towards elucidating and visualizing community-level functional connections and the contribution of microorganisms to biogeochemical cycles (Figure 7A). Additionally, in order to compare the prediction performance of METABOLIC to others, we conducted parallel in silico experiments (Figure 7B). We used two representative bacterial genomes as the test datasets. We randomly picked 100 protein sequences from individual genomes and submitted them to annotation by these six software/online servers. Predicted protein annotations by individual software and online servers were compared to their original annotations that were provided by the NCBI database (Additional file 11, 12: Dataset S4, S5). According to statistical methods of evaluating binary classification [73], the following parameters were used to make the comparison: 1) recall (also referred to as the sensitivity) as the true positive rate, 2) precision (also referred to as the positive predictive value) which indicates the reproducibility and repeatability of a measurement system, 3) accuracy which indicates the closeness of measurements to their true values, and 4) F1 value which is the harmonic mean of precision and recall, and reflects both these two parameters. Among the tested software/online servers, the performance parameters of METABOLIC consistently placed it as the top 3 and top 2 software for recall and F1 and the top 1 and top 2 software for precision and accuracy. These results demonstrate that METABOLIC (Figure 7B) provides robust performance and consistent metabolic prediction that facilitate accurate and reliable applicability for downstream data visualization and community-level analyses.
To demonstrate the application and performance of METABOLIC in different samples, we tested eight distinct environments (marine subsurface, terrestrial subsurface, deep-sea hydrothermal vent, freshwater lake, gut microbiome from patients with colorectal cancer, gut microbiome from healthy control, meadow soil, wastewater treatment facility). Overall, we found METABOLIC to perform well across all the environments to profile microbial genomes with functional traits and biogeochemical cycles (Additional file 10: Dataset S3). Among these tested environments, we also performed community-scale metabolic comparisons based on the MW-score (Figure 8). MW-score fraction at the community scale reflects the overall metabolic profile distribution pattern. Specifically, we compared samples from terrestrial and marine subsurface and samples from hydrothermal vent and freshwater lake. We observed that terrestrial subsurface contains more abundant metabolic functions related to nitrogen cycling compared to the marine subsurface (Figure 8A), consistent with the previous characterization of these two environments [2, 74]. Deep-sea hydrothermal vent samples had a considerably high concentration of methane and hydrogen [52] as compared to Lake Tanganyika (freshwater lake). Consistent with this phenomenon, the deep-sea hydrothermal vent microbial community has more abundant metabolic functions associated with methanotrophy and hydrogen oxidation (Figure 8B). In order to focus on a specific biogeochemical cycle, we applied METABOLIC to compare sulfur-related metabolisms at the community scale for these two environment pairs (Additional file 7: Figure S7). Terrestrial subsurface contains genomes covering more sulfur cycling steps compared to marine subsurface (7 steps vs 3 steps) (Additional file 7: Figure S7A). Freshwater lake contains genomes involving almost all the sulfur cycling steps except for sulfur reduction, while deep-sea hydrothermal vent contains less sulfur cycling steps (8 steps vs 6 steps) (Additional file 7: Figure S7B). Nevertheless, deep-sea hydrothermal vent has a higher fraction of genomes (59/98) and a higher relative abundance (73%) of these genomes involving sulfur oxidation compared to the freshwater lake (Additional file 7: Figure S7B). This indicates that the deep-sea hydrothermal vent microbial community has a more biased sulfur metabolism towards sulfur oxidation, which is consistent with previous metabolic characterization on the dependency of elemental sulfur in this environment [52, 75–77]. Collectively, by characterizing community-scale metabolism, METABOLIC can facilitate the comparison of overall functional profiles as well as for a particular elemental cycle.
METABOLIC enables accurate reconstruction of cell metabolism
To demonstrate applications of reconstructing and depicting cell metabolism based on METABOLIC results, two microbial genomes were used as an example (Figure 9). As illustrated in Figure 9A, Hadesarchaea archaeon 1244-C3-H4-B1 has no TCA cycling gene components, which is consistent with previous findings in archaea within this class [78]. Gluconeogenesis/glycolysis pathways are also lacking in the genome; since gluconeogenesis is the central carbon metabolism responsible for generating sugar monomers which will be further biosynthesized to polysaccharides as important cell structural components [79], the lack of this pathway could be due to genome incompleteness. As an enigmatic archaeal class newly discovered in the recent decade, Hadesarchaea have distinctive metabolisms that separate them from conventional euryarchaeotal groups. They almost lost all TCA cycle gene components for the production of acetyl-CoA; while they could metabolize amino acids in a heterotrophic lifestyle [78]. It is posited that the Hadesarchaea genome has been subjected to a streamlining process possibly due to nutrient limitations in their surrounding environments [78]. Due to their metabolic novelty and limited available genomes at the current time, there are still uncertainties on unknown/hypothetical genes and pathways and unclassified metabolic potential across the whole class. The previous metabolic characterization on four Hadesarchaea genomes indicates Hadesarchaea members could anaerobically oxidize CO, and H2 was produced as the side product [78]. In the Hadesarchaea archaeon 1244-C3-H4-B1 genome, METABOLIC results indicate the loss of all anaerobic carbon-monoxide dehydrogenase gene components, which suggests the distinctive metabolism of this Hadesarchaea archaeon from others and highlights the accuracy of METABOLIC in reflecting functional details.
We also reconstructed the metabolism for Nitrospirae bacteria M_DeepCast_50m_m2_151, a member of the Nitrospirae phylum reconstructed from Lake Tanganyika [55] (Figure 9B). It contains the full pathway for the TCA cycle and gluconeogenesis/glycolysis. Furthermore, it also has the full set of oxidative phosphorylation complexes for energy conservation and functional genes for nitrite oxidation to nitrate. Other nitrogen cycling metabolisms identified in this genome include ammonium oxidation, urea utilization, and nitrite reduction to nitric oxide. The reverse TCA cycle pathway was identified for carbon fixation. The metabolic profiling result is in accord with the fact that Nitrospirae is a well-known nitrifying bacterial class capable of nitrite oxidation and living an autotrophic lifestyle [79]. Additionally, their more abundant distribution in nature compared to other nitrite-oxidizing bacteria such as Nitrobacter indicates their significant contribution to nitrogen cycling in the environment [79]. This highlights the ability of METABOLIC in reflecting functional details of more common and prevalent microorganisms compared to the Hadesarchaea archaeon. Notably as discovered from METABOLIC analyses, this bacterial genome also contains a wide range of transporter enzymes on the cell membrane, including mineral and organic ion transporters, sugar and lipid transporters, phosphate and amino acid transporters, heme and urea transporters, lipopolysaccharide and lipoprotein releasing system, bacterial secretion system, etc., which indicates its metabolic versatility and potential interactive activities with other organisms and the ambient environment. Collectively, METABOLIC result of functional profiling provides an intuitively-represented summary of a single microbial genome which enables depicting cell metabolism for better visualizing the functional capacity.
METABOLIC accurately represents metabolism in the human microbiome
In addition to resolving microbial metabolism and biogeochemistry in environmental microbiomes, METABOLIC also accurately identifies metabolic traits associated with human microbiomes. The implications of microbial metabolism on human health largely remain a black box, much like microbial contributions to biogeochemical cycling. We demonstrate the utility of METABOLIC in human microbiomes using publicly available data from stool samples collected from patients with colorectal cancer and healthy individuals. From this study, we selected stool metagenomes from one colorectal cancer (CRC) and an age and sex-matched healthy control to conduct the comparison. The heatmap indicates the human microbiome functional profiles of both samples based on the marker gene presence/absence patterns (Figure 10). As an example of METABOLIC’s application, we demonstrate that there were 28 markers with variations > 10% in terms of the marker-containing genome fractions between these two samples (Figure 10, Additional file 13: Dataset S6). These 28 markers involved all the ten metabolic categories except for lipid metabolism and translation, suggesting the broad functional differences between these two samples. Recently, METABOLIC was applied to stool metagenomic samples from 667 individuals who either were healthy or had adenomas or carcinomas of the colon, to profile organic/inorganic sulfate reduction and sulfide production [80]. In addition to analyzing human microbiome specific functional markers, METABOLIC can be used to visualize elemental nutrient cycling and analyze metabolic interactions in human microbiomes. Overall it enables systematic characterization of the composition, structure, function, and dynamics of microbial metabolisms in the human microbiome and facilitates omics-based studies of microbial community on human health [59].