Introduction

Bacteria exhibit metabolic diversity and can utilise a broad range of substrates for growth. It has become clear amongst pathogens that there is an intertwined relationship between metabolism and nutrient usage with virulence and antimicrobial resistance (1-7). Comparative analyses of metabolic profiles (e.g. substrate usage) are key to fully understanding these relationships. Traditionally, these profiles have been assessed via phenotypic growth on a limited number of substrates, such as those used to delineate between species (8-10) which form the basis of a number of commercial products for species identification. However, these methods are not sufficiently discriminatory for in-depth comparisons within species, and alternative approaches such as the Omnilog Phenotype MicroArray system (Biolog) are too expensive and/or labour intensive for application to large numbers of isolates. Similarly, probing of essential metabolism-associated genes via transposon mutant libraries (e.g. to identify novel virulence factors and therapeutic targets) (4, 11, 12) cannot be easily scaled across diverse bacterial populations.

Genome-scale metabolic models or metabolic reconstructions are a computational approach to analysing the metabolic potential of an organism, within which the entire biochemical network is represented as a stoichiometric matrix (13). Metabolic models are constructed programmatically, but typically informed and at least partially validated using phenotypic growth data (14-16). Once constructed, they can be run through simulations and analysed under various contexts, such as in silico growth experiments (Flux Balance Analysis [FBA]) to predict substrate usage profiles (17), evaluate the impact of single gene knockouts on growth (14, 18), and identify metabolic chokepoints for drug targets (19), among others. Traditionally, metabolic models are strain-specific (i.e. each model represents a unique individual http://bigg.ucsd.edu/models) and may not be applicable to other isolates due to unrepresented genetic diversity.

We recently described 37 curated strain-specific models for the Klebsiella pneumoniae Species Complex (KpSC) (14) comprised of K. pneumoniae and its close relatives (20). These organisms are a common cause of healthcare-associated infections world-wide, and among the World Health Organization’s priority antimicrobial resistant pathogens (21). KpSC are highly diverse and gene content can differ substantially between strains (22, 23). Accordingly, our models varied in terms of gene and reaction content, resulting in variable growth substrate usage profiles and metabolic redundancy (14). Similar variation has also been described in other key bacterial pathogens e.g. Escherichia coli (24), Salmonella enterica (25), Staphylococcus aureus (26) and Pseudomonas aeruginosa (27). This is highly relevant to the use of metabolic models for the exploration of virulence and antimicrobial resistance, and for the identification of novel drug targets. Therefore, such works should seek to include multiple strain-specific models, and in some cases 100s-1000s of models may be required to accurately represent population diversity (22, 28, 29).

There are several open source tools currently available that can rapidly produce strain-specific metabolic models, including CarveMe (30), gapseq (31), ModelSEED (32) and KBase (33) (see the recent review by Mendoza and colleagues for comparative descriptions (34)), as well as a recently published modelling and analysis pipeline, ChiMera, which leverages CarveMe for model construction (35). In their systematic analysis Mendoza et al. indicated CarveMe and ModelSEED to be of particular interest for large-scale studies due to their speed and model quality (34) (note that gapseq was published later and therefore not included in this review). KBase implements ModelSEED (32) as a web interface application, limiting its utility for high-throughput analysis of 100s – 1000s of bacterial genomes. CarveMe is a command line application; it is open source but is dependent on commercial solvers such as CPLEX (free for academic use). However, its use of a universal reference model may limit specificity of strain-specific models (36), and result in overestimation of model genes. These limitations can be overcome by manual curation of the output models, but such curation is highly labour intensive and not suitable for high-throughput analyses. Furthermore, the CarveMe database (BiGG universal_model) appears to be no longer actively maintained, meaning that there is no opportunity to integrate novel structural and/or biochemical data as these become available in the literature (as discussed in COBRA community forums). gapseq is a recently published command line tool which leverages an independent universal database (31). In their comparative analyses, the gapseq authors demonstrate superior accuracy to both CarveME and ModelSEED; however the concerns about the specificity of universal models remain and it has been reported that model construction takes considerable time (several hours in some cases (37)).

Here, we present Bactabolize (available at https://github.com/kelwyres/Bactabolize), an easy-to-use tool which allows scalable production of strain-specific draft metabolic models and prediction of growth phenotypes in under three minutes per input genome. Bactabolize builds upon the reference-based model reconstruction approach described by Norsigian et al. (36), leveraging the COBRApy framework (38) and BiGG nomenclature (39). We present a pan-metabolic reference model for the KpSC (derived from our 37 curated strain-specific models (14)), and describe an exemplar quality control framework for the application of Bactabolize to KpSC draft genome assemblies. We show that Bactabolize can rapidly produce strain-specific models from draft genomes with a high degree of completeness (as compared to models generated from completed genome assemblies), resulting in highly accurate growth predictions that match or exceed the accuracy of models from CarveMe and manual curation efforts.

Results

Description of Bactabolize

Bactabolize is written in Python 3 and utilises the metabolic modelling library COBRApy (38). Bactabolize has four main commands:

• i) Draft model generation (draft_model command), which generates a strain-specific draft metabolic reconstruction (‘model’) using the approach outlined previously (36), and uses gap-filling to identify any missing reactions required to simulate growth in the user-specified conditionsdraft_model command), which generates a strain-specific draft metabolic reconstruction (‘model’) using the approach outlined previously (36), and uses gap-filling to identify any missing reactions required to simulate growth in the user-specified conditions

• ii) Patching incomplete models (patch_model command) by the addition of missing reactions e.g. those identified by the automated gap-filling processpatch_model command) by the addition of missing reactions e.g. those identified by the automated gap-filling process

• iii) Substrate usage analysis via Flux Balance Analysis (FBA) (fba command) to predict growth outcomes for a specified range of substrates supported by the model(s) Figure 1).fba command) to predict growth outcomes for a specified range of substrates supported by the model(s) Figure 1).

Figure 1:

Additional processing scripts are provided alongside Bactabolize to improve model metadata annotation (improve_model_annotations.py), convert models generated using KBase and ModelSEED to Bactabolize/BiGG-compatible format (SEED_to_BiGG_model_convert.sh), generate network graph files from models (model_to_network_graph.py) and merging output FBA profiles (merge_fba_profiles_longtable.sh). Full documentation, including example code and test data are available at the Bactabolize code repository (https://github.com/kelwyres/Bactabolize).

For draft model construction, Bactabolize requires users to provide an input assembly (annotated or unannotated FASTA or Genbank format respectively), a reference model (JSON format) and the corresponding reference sequence data (gene and protein sequences in two separate multi-fasta files or a single Genbank annotation in a .gbk file) (Figure S1). For optimum results we suggest using a pan-model that captures as much diversity as possible for the target species or group of interest, because Bactabolize’s reconstruction method is reductive i.e. each output strain-specific model will include only genes, reactions and metabolites that are present in the reference or a subset thereof (although novel genes, reactions and metabolites can be added via manual curation).

If the input assembly is unannotated, Bactabolize will identify coding sequences using Prodigal (40) but will otherwise honour the existing coding sequence (CDS) notations and optionally use Prodigal to search for additional CDS. Draft genome-scale metabolic models are output in both SMBL v3.1 (41) and JSON formats (one pair of files for each independent strain-specific model), along with an optional MEMOTE quality report (42). Bactabolize will identify orthologs in the input genome(s) compared to the reference sequence data using Bi-directional BLAST (43) Best Hits (BBH) (44) using BLAST+ (36). Users can parameterise the ortholog finding settings (coverage and identity thresholds) for BBH. Alternatively, there is the option of using protein similarity to identify orthologs instead of identity.

Once a draft model has been constructed, it is validated via a simulated growth experiment on user-input choice of media and atmosphere (aerobic or anerobic). Predefined media include BG11 (Gibco), M9 + glucose (36), nutrient media (45), Luria-Bertani (LB) (45), Tryptic Soy (TSA) (45), TSA + sheep blood (45), LB as specified by the CarveMe developers (30), Chemically Defined Medium (CDM)-like (34), Plantarum Minimal Medium (PMM) PMM5-like (34) and PMM7-like (34). Users can also define custom media as Bactabolize supports several complex media ingredients, including peptone (peptic digest of bovine and porcine tissue) (46-48), tryptone (pancreatic digest of casein) (46, 48, 49), soy peptone/soytone (digest of soymeal) (46, 48, 50, 51), yeast extract (52-57) and beef extract (46, 48). If the model fails to simulate growth, gap-filling is performed to indicate missing reactions. Users can add these reactions to a patch JSON file and optionally use the patch_model command to correct the model (Figure S2). Bactabolize uses a conservative gap-filling approach that only adds the minimum number of reactions to enable growth under the chosen conditions. We recommend testing the models in minimal media and atmosphere expected to support growth for all isolates of the species of interest, unless the user has access to matched phenotypic data demonstrating growth for individual isolates in specific conditions. Aggressive gap-filling will effectively homogenise the models and should be avoided if the goal is to understand the underlying strain diversity.

Substrate usage analysis (the fba command) is performed iteratively for each possible carbon, nitrogen, sulfur and phosphor substrate supported by the model(s) (Figure S3), by replacing the default substrate in the user specified growth medium (specified in the fba_spec JSON file). For example, in M9 media the default substrates are glucose (carbon), ammonia (nitrogen), sulphate (sulfur) and phosphate (phosphor). Each substrate can be tested in aerobic and/or anaerobic conditions. Growth prediction output is recorded in a tab delimited file (one per strain). The merge_fba_profiles_longtable.sh helper script will combine the outputs for multiple strains into a single file for downstream analysis.

The growth impacts of single-gene knockout mutations can be simulated via the sgk command (Figure S4). Bactabolize will iterate through every gene in the model, temporarily removing it and its associated reactions (unless they are also associated with another gene) and running FBA to simulate growth in the user-specified conditions. The output is comparable to single-gene knockout studies such as transposon mutagenesis and can be used to probe gene essentiality.

KpSC pan-metabolic reference model

We constructed a species complex-specific pan-metabolic reference model by combining a collection of 37 manually curated models for which we have previously demonstrated high accuracy (range 88.3%–96.8% for prediction of 94 distinct growth phenotypes (14)). These models represent a diverse collection of KpSC (14) (including at least one each of the seven major taxa in the complex; K. pneumoniae, Klebsiella variicola subsp variicola, Klebsiella variicola subsp tropica, Klebsiella quasipneumoniae subsp quasipneumoniae, Klebsiella quasipneumoniae subsp similipneumoniae, Klebsiella quaisivariicola, Klebsiella africana). The combined pan-model, known as KpSC-pan v1, comprises a total of 1265 distinct genes, 2319 reactions and 1696 metabolites, and is available at github.com/kelwyres/KpSC-pan-metabolic-model.

Performance comparison

We compared the output and performance of Bactabolize to the two previously published tools that can support high-throughput analyses i.e. CarveMe (30) and gapseq (31). To aid interpretation in the context of community standard approaches, we also include a comparison to the popular web-based reconstruction tool, KBase (ModelSEED), and a manually curated metabolic reconstruction of K. pneumoniae strain KPPR1 (also known as VK055 and ATCC 43816, metabolic model named iKp1289) (15). This isolate was chosen as there is a completed genome sequence (Genbank accession: CP009208), single-source growth phenotype (15) and single-gene knockout growth essentiality data available (58). De novo draft models for strain KPPR1 were built using; i) Bactabolize with the KpSC pan v1 reference; ii) CarveMe, with its universal reference model (CarveMe universal); iii) CarveMe, with KpSC-pan v1 reference (CarveMe KpSC pan); iv) gapseq; and v) KBase (ModelSEED). Importantly, neither K. pneumoniae KPPR1 nor its genetic lineage (7 gene multi-locus sequence type, ST493), are represented in the KpSC pan-reference model, meaning these benchmarking comparisons were on equal footing. Subsequently, each model was used to predict growth phenotypes; i) in M9 minimal media with different sole sources of carbon, nitrogen, phosphorus and sulfur; and ii) for all possible single gene knockouts in LB under aerobic conditions. The predicted phenotypes were compared directly to the published phenotype data.

Among the high-throughput approaches, the Bactabolize draft model captured fewer genes and reactions (n = 1233 and 2307, respectively) than the gapseq (n = 1489 and 3186, respectively) and CarveMe universal models (n = 1960 and 2857, Figure 2A). The Bactabolize and CarveMe universal models contained similar numbers of metabolites (1696 vs 1737, respectively), while the gapseq model contained many more (n = 2519). Notably, the initial gapseq model was gap-filled by simulation of growth in M9 minimal media with glucose (consistent with Bactabolize and CarveMe media recipes), which resulted in the addition of 31 reactions to the draft model. This is in contrast to CarveMe and Bactabolize-generated models which produced biomass in M9 minimal media plus glucose without additional gap-filling. The CarveMe KpSC pan model captured considerably more genes than any of the other models (n = 2407), but these were associated with many fewer unique reactions and metabolites (1206 and 825, respectively). Upon further investigation we determined that this method resulted in the over prescription of gene reaction rules (GPRs) to multiple reactions (mean 2.2 GPRs per reaction when compared to Bactabolize using the same pan reference model: 1.94 GPRs per reaction.

Figure 2

Figures S5 and S6 show the overlaps of metabolites and reactions between the high-throughput reconstruction methods after processing with MetaNetX (59) to standardise the reaction and metabolite nomenclatures (excluding CarveMe pan for simplicity and given the likely problems of reaction oversubscription). The majority of the reactions included in the Bactabolize model were conserved in either the CarveMe universal model (n = 1225, 53.2%), gapseq model (n = 54, 2.3%) or both (n = 665, 28.9%). The reaction overlap was skewed to the CarveMe universal model which shared 1225 reactions that were conserved in the Bactabolize model but absent from the gapseq model. Notably, the gapseq model contained a large number (2200) of unique reactions (70.4% of those in the model). Similarly, the vast majority of metabolites in the Bactabolize model were conserved in one or both of the other models (n = 917, 85.6%). However, it is likely that true overlaps between methods are underrepresented due to the different reaction identifiers and chemical synonyms used within the BiGG (Bactabolize, CarveMe) vs ModelSEED nomenclatures (gapseq), which are difficult to harmonise in an automated manner even after the application of MetaNetX.

In comparison to the low-throughput reconstruction approaches, the Bactabolize model contained a similar number of genes, reactions and metabolites to the manually curated model (n = 1289, 2484 and 1827, respectively) and many more than the KBase model (72, 544 and 534, respectively). This is likely due to the low number of metabolism-associated genes identified by KBase, which has impacted the associated reaction and metabolite data.

MEMOTE scores, (produced by the MEMOTE report (42)) indicate the quality of the model metadata annotations, with the scores ranging between 0 – 100%. These provide a measure of model portability and the level of connected databases available to support the metabolite, reaction and genetic information represented in the model, but bear no reflection on model accuracy. Bactabolize performs on the lower end, with CarveMe universal and gapseq performing the best (Figure 2B). The KBase model appears to perform well in this regard, however this is due to the low number of genes, reactions and metabolites included in the model. Bactabolize using the KpSC-pan model outperforms the model propagation mode of CarveMe using the same reference model (Figure 2B). Work is ongoing to improve the metadata annotations in the KpSC-pan reference model, to support large-scale model propagation.

We assessed the performance of each model for in silico prediction of growth phenotypes compared to the previously published experimental data (15). Accuracy, sensitivity, specificity, precision and F1 scores were calculated (60). Note that the specific set of growth substrates and gene knockouts that can be simulated is determined by the sets of genes and metabolites captured by each model and is therefore model-dependent (Data S1 and S2). Among those with matched experimental phenotype data, the Bactabolize and CarveMe universal models were able to predict growth for a greater number of carbon, nitrogen, phosphorous and sulfur substrates than gapseq, CarveMe KpSC pan, KBase and iKp1289 models (Figure 2C, Data S1). While the CarveMe universal model had the highest number of true-positive growth predictions overall (n = 132 of 617 total predictions), it also had a comparably high number of false-positive predictions (n = 39 of 617 total predictions, Figure 2D). Similarly, the gapseq and iKp1289 models resulted in 31 (262 total predictions) and 50 (513 total predictions) false-positive predictions, respectively. In contrast, the Bactabolize model had fewer false-positive predictions (n = 21 of 505 total predictions) alongside a high number of true-positive predictions (n = 117 of 505 total predictions), resulting in the highest overall accuracy metrics (Figure 2E, Data S1). The KBase model was a notable outlier, associated with a high number of false-negative predictions (n = 31 of 103 total predictions) and low false-positive predictions (n = 3 of 103 total predictions), presumably resulting from the very low number of genes and reactions included in the model, driving low sensitivity and accuracy.

The gene essentiality results showed that gapseq produced the highest absolute number of true-positive gene essentiality predictions (n = 79 of 1403), followed by Bactabolize KpSC pan (n = 44 of 1220 total predictions), then CarveMe universal (n = 39 of 1951 total predictions). CarveMe universal had the largest number of true-negatives by a wide margin (n = 1599 of 1951 total predictions), followed by gapseq (n = 1085 of 1403 total predictions), then Bactabolize KpSC pan (n = 939 of 1220 total predictions), driving their high accuracies (83.96%, 82.96% and 80.57%, respectively). The Bactabolize model was associated with the greatest overall precision and specificity (Figures 2F & 2G) while the gapseq model resulted in the highest F1-score and sensitivity.

While model features and accuracy are essential metrics for comparison, computation time is also a key consideration for high-throughput analyses. We recorded the time required for each tool to build draft models for 10 of the completed KpSC genomes used in the quality control framework (see below) on a high-performance computing cluster (Intel Xeon Gold 6150 CPU @ 2.70GHz and 155 GB of requested memory on a CentOS Linux release 7.9.2009 environment. CarveMe KpSC pan was the fastest with a mean of 20.04 (range 19.90 - 20.18) seconds, followed by CarveMe universal at 30.28 (range 29.20 - 31.80) seconds, then Bactabolize KpSC pan at 98.05 (range 92.19 - 100.4) seconds. KBase took 183.50 (range 120.00 - 338.00) seconds per genome via batch analysis, including genome upload time and queuing. gapseq took 5.46 (range 4.55 - 6.28) hours to produce draft models (not including the required gap-filling), consistent with previous reports (37).

Quality control framework for input genome assemblies

There are now thousands of bacterial genomes available in public databases, the majority of which are in draft form, comprising 10s to 1000s of assembly contigs. This fragmentation of the genome is caused by repetitive sequences that cannot be resolved by the assembly algorithm and/or sequence drop-out. The latter can result in the loss of genetic information such that some portion of genes present in the underlying genome are lost from the genome assembly (either completely or partially). This in turn, poses a limitation for the reconstruction of metabolic models using these assemblies, since most published approaches use sequence searches to predict the presence/absence of genes and their associated enzymatic reactions. Therefore, if we are to use public genome data for high-throughput metabolic modelling studies, it is essential to evaluate the expected model accuracies and understand the minimum input genome quality requirements.

Here we performed a systematic analysis leveraging our published curated KpSC models (n=37, (14)), which were generated using completed genome sequences and were therefore considered to represent ‘complete’ models for which the underlying genome sequence contains all genes that are truly present in the genome (note the biological accuracy of these models was reported previously (14) and is not the subject of the current study). We randomly subsampled the corresponding Illumina read sets to various depths (10 – 100x, increments of 10) in triplicate and generated draft assemblies that were passed to Bactabolize for generation of draft metabolic models (Data S3). Due to low read depth (≤30x), two isolates were removed from this analysis. Additionally, ten 10x depth read samples failed to produce assemblies, leaving 1040 draft genomes for analysis. The resulting draft metabolic models were compared to the complete models to; i) determine the proportions of complete model genes and reactions captured in the draft models; and ii) compare 846 in silico aerobic growth predictions in M9 minimal media, where growth on 266 carbon, 153 nitrogen, 59 phosphorus and 25 sulfur sources were examined. Substrates containing multiple elements were tested as sole sources of each element independently and in combination, e.g. 1,5-Diaminopentane was tested as a sole carbon, sole nitrogen and sole carbon plus nitrogen source.

As expected, assembly quality generally increased with increasing sequencing depth i.e. assemblies generated from higher depth read sets were associated with higher N50 values, fewer contigs and fewer assembly graph dead-ends, although the rate of improvement drastically declined beyond 40-50x depth (Figure S7, Data S3). We noted that it was rare for draft models to capture 100% of model genes and reactions (just 420 of all 1040 draft assemblies were associated with models that captured 100% model genes) (Data S3, Figure S8), even when using the highest quality draft genomes. However, ≥99% of genes and reactions were commonly captured, which plateaued from 40x depth onwards (Figure S8). Therefore, we sought to evaluate whether ≥99% model capture would produce functionally accurate models.

We used FBA to simulate substrate growth profiles for the 40x depth assemblies, representing a sequencing depth that can be routinely achieved with standard Illumina library preparations. All but one assembly triplicate set (isolate SB4767 98% gene capture, 99% reaction capture) captured ≥99% but ≤100% model genes and/or reactions. The substrate growth profiles were then compared to those of the complete models. The vast majority of draft models produced accurate growth predictions; 102 of 108 models resulted in predictions with 100% concordance to those from the corresponding complete models. Three models for K. quasipneumoniae similipneumoniae isolate SB164 resulted in predictions with a mean of 99.8% concordance. The remaining three models were for isolate SB4767 and resulted in mean of 80.4% concordance. Notably, these models were those representing <99% gene capture. Together, these data suggest that draft models capturing ≥99% of the complete model genes/reactions generate highly accurate growth predictions and that these capture rates can be readily achieved from draft genome assemblies.

We investigated the relationships between assembly quality metrics and model gene/reaction capture in more detail. Variation in assembly graph dead-ends accounted for the greatest amount of variation in model capture, closely followed by raw contig counts (cubic polynomial fit, R² of ≥0.98 for graph dead-ends, R² of ≥0.9 for contig count). A segmented linear model was fitted to N50 length (R² ≥ 0.83), producing a breakpoint at 25153 bp (Figure 3).

Figure 3:

To further explore the optimum thresholds for assembly metrics, we tallied the number of draft assemblies resulting in ≥99% and <99% gene and reaction capture at increasing graph dead-end and contig count cut-offs, and decreasing N50 cut-offs. Draft models that captured ≥99% of the complete model genes/reactions were considered ‘good’ models, whereas draft models that captured <99% of complete model genes/reactions were considered ‘bad’ models. The optimum threshold for assembly graph dead end was determined to be ≤200. At this value, 94.44% of ‘good’ models were captured, and 0% ‘bad’ models. The optimum threshold for contig counts was determined as ≤130 contigs at which 67.92% of ‘good’ and 0% ‘bad’ models were captured (Figure 4). The optimum threshold for N50 was determined to be ≥65000, at which 94.97% of ‘good’ and 1.71% of ‘bad’ models were captured. The assembly graph dead-end threshold results in comparatively higher sensitivity (i.e. a higher proportion of ‘good’ models pass the threshold) than contig count and comparatively better specificity (i.e. lower proportion of ‘bad’ models pass the threshold) than N50, but the underlying metric information is not universally available because many isolate genomes are deposited in public databases only as assemblies without the associated assembly graph. We therefore recommend a three-tier approach, whereby the assembly graph dead-end criterion is preferenced if available, followed by N50 and then contig count.

Figure 4:

Impact of gap-filling models

Of the 901 draft genome assemblies which passed our QC criteria (≤200 assembly graph dead ends), 23 of the resulting draft models failed to simulate growth in M9 minimal media with glucose (despite capturing ≥99% of the genes and reactions in the corresponding complete models). It is expected that all KpSC models should be able to simulate growth on M9 media with glucose as a sole carbon source, as this central metabolism is universal amongst KpSC. To replace missing, critical reactions required for growth on M9 with glucose, we investigated model gap-filling using the patch_model command of Bactabolize. We then assessed the accuracy of the gap-filled models for prediction of growth on the full range of substrates, as compared to the predictions from the corresponding complete models.

Gap filling added 1 – 3 missing reactions to each model, with a median of one, fully restoring biomass production in M9 media with glucose in all but two of the 23 failed models. The missing reactions appeared to be random genes across these 23 genomes, likely due to missing information in these assemblies.

Substrate usage predictions from the 21 successfully gap-filled models were highly accurate, with 18/21 having a prediction concordance of ≥99% across all 846 growth conditions (12/21 had 100% concordance) (Figure S9). We therefore conclude that models generated for genome assemblies passing our QC criteria, which have been gap-filled to successfully simulate growth on minimal media plus glucose, are suitable for the prediction of growth across a range of substrates.

Predictive accuracy of draft models

We assessed the accuracy of Bactabolize for the construction of draft models for 10 novel KpSC clinical isolates, representing five of the major taxa in the complex. We included five isolates for which the associated STs were represented in the KpSC-pan v1 model and five isolates with STs that were not represented. Whole genome sequence data were generated on the Illumina platform and draft assemblies generated de novo. The resultant assemblies had 0-4 graph dead-ends, N50s of 151958-388486 bp and 83-187 contigs (Data S4), within the tiered threshold values.

FBA was performed, and the predicted growth profiles compared to matched phenotypic growth data for 16 carbon sources derived from Vitek GN ID cards. Though the number of tested carbon sources was limited, all were associated with high accuracy metrics (Figure 5, Data S4). As expected, models for isolates with STs represented in the KpSC-pan v1 reference performed slightly better (mean accuracy = 0.98) than those for non-represented STs (mean accuracy = 0.95).

Figure 5:

Discussion

In this work we described Bactabolize, a pipeline for rapid and scalable production of accurate bacterial strain-specific metabolic models and growth phenotype predictions. We describe a pan-reference model for the KpSC and demonstrate that a draft strain-specific model generated de novo via Bactabolize using the KpSC-pan v1 reference was highly accurate for growth phenotype prediction (85.79% accuracy for substrate usage across 190 substrates, and 80.57% for gene essentiality across 1220 genes). Importantly, we also described a quality control framework for the use of draft genome assemblies as input for metabolic reconstructions. We used a systematic analysis to; i) evaluate the proportion of gene and reaction capture compared to the corresponding ‘completed’ models; ii) define quality control thresholds for input assemblies (three tier approach for KpSC; ≤200 assembly graph dead ends, followed by ≥65000 N50, followed by ≤130 contigs); and iii) estimate the accuracy of the resultant growth predictions. While the quality control thresholds and accuracy estimates are specific to KpSC, the conceptual framework can be applied to any organism and is essential to support the confident application of metabolic modelling for large-scale genome datasets. We appreciate that assembly graphs may not be available for dead end count, e.g. for draft genome assemblies accessed via public repositories, however we encourage users to include this information in their quality control procedures wherever possible (e.g. using the recently published counter tool available at https://github.com/rrwick/GFA-dead-end-counter), because these counts represent a direct reflection of the completeness of the genome assembly. In contrast, contig counts and N50 are influenced by biological features such as repeat copy numbers as well as the underlying sequence data quality e.g. a bacterial genome harbouring many insertion sequence insertions will result in a draft assembly with a high number of contigs regardless of the sequence data quality and completeness.

Traditionally, genome-scale metabolic reconstruction approaches have relied upon significant manual curation efforts. While there will always remain a need for high quality curated models, such resource intensive approaches preclude their application at scale, and have therefore limited analyses to small numbers of individual strains (15, 16). However, automated reconstruction approaches can support the generation and comparison of multiple strain-specific draft models from which meaningful biological insights can be derived (61). Additionally, the quality of curated models is likely to vary depending on their age, level and type of curation, as well as the approach used for preliminary drafting. Indeed it is possible for automated approaches to outperform manually curated models; a draft model for K. pneumoniae KPPR1 generated using Bactabolize with the KpSC pan-v1 reference model outperformed the manually curated iKp1289 model representing the same strain (15). iKp1289 was published in 2017 (6 years prior to this study) and was initially drafted via the KBase pipeline (33), which uses RAST to annotate the sequences with Enzyme Commission numbers. It has been demonstrated several times that the Enzyme Commission scheme has systematic errors (62, 63), leading to a loss in accuracy when compared to the ortholog identification methods used by automated approaches. Consistent with this assertion, our draft KPPR1 model constructed with KBase (without manual curation) was an outlier in terms of the very low number of genes, reactions and metabolites that were included.

CarveMe with universal model (30) and gapseq (31) are the current gold standard automated approaches for model reconstruction, and we show that a draft KpSC model generated by Bactabolize with the KpSC pan v1 reference resulted in similar or better accuracy for phenotype prediction (Figure 2). Both the CarveMe universal and gapseq models resulted in high numbers of true-positive and true-negative growth predictions. However, these were also accompanied by comparatively higher numbers of false-positive predictions that resulted in a lower overall accuracy for substrate usage analysis compared to Bactabolize with the KpSC-pan v1 reference (Figure 2), and comparatively lower precision and specificity for the gene essentiality analysis. False-positive predictions may indicate that the relevant metabolic machinery are present in the cell but were not active during the growth experiments (e.g. due to lack of gene expression). In this regard, false-positives are not always a sign of model inaccuracy. However, false-positive predictions can also occur from incorrect gene annotations e.g. due to reduced specificity of ortholog assignment resulting from the use of the universal model without manual curation. Given a key objective here is to facilitate high-throughput analysis for large numbers of genomes, it is not feasible to expect that all models will be manually curated, and therefore we believe that identifying fewer genes with lower overall error rates provides greater confidence in the resulting draft models. We also note that the BiGG universal reference model which CarveMe leverages is no longer being actively maintained. In contrast, user defined reference models can be iteratively curated and updated to incorporate new knowledge and data as they become available.

Bactabolize’s reference-based reconstruction approach is reductive, meaning the resultant draft models will comprise only the genes, reactions and metabolites present in the reference, or a subset thereof, and will not include novel reactions unless they are manually identified and curated by the user. This is an important caveat that should be considered carefully for application of Bactabolize to large genome data sets, particularly for genetically diverse organisms such as those in the KpSC. For optimum results we suggest using a curated pan-model that captures as much diversity as possible for the target species or group of interest. While we acknowledge that a reasonable resource investment is required to generate a high-quality reference, we have shown that a pan-model derived from just 37 representative strains can be sufficient to support the generation of highly-accurate draft models (Figure 2 and 5). Additionally, we note that it is possible to use a single strain reference model, which should ideally represent the same or closely related species to that of the input genome assemblies, in order to facilitate accurate identification of gene orthologs. It is technically possible to use an unrelated reference model, but this is expected to result in inaccurate and/or incomplete outputs and has not been tested in this study. In circumstances were no high quality closely-related reference model is available, we recommend alternative reconstruction approaches that leverage universal databases e.g. CarveMe (30) or gapseq (31). However, gapseq’s long compute time makes it inappropriate for application to datasets comprising 100s-1000s of genomes (such as have become increasingly common in the bacterial population biology literature) (64, 65).

Bactabolize and the KpSC pan v1 model are freely available under open source licenses and satisfy the four features of the FAIR research principles (findability, accessibility, interoperability and reusability) (66). In addition to the KpSC pan-reference described here, a pan-reference model has been described previously for Salmonella enterica (representing 410 strains) (25), Bacillus subtilis (67) (representing 183 strains) and Escherichia coli (68) (representing 222 strains). We are actively working to expand and improve the KpSC pan reference model and welcome similar efforts to generate high quality references for other organisms. Together these resources will facilitate population wide metabolic analyses for global priority pathogens, which can be used to understand how they transmit, cause disease and evolve drug-resistance, and to identify novel therapeutic targets.

Methods

Bactabolize pipeline

Bactabolize utilises the existing metabolic modelling library COBRApy (38) and Python 3 (69). All code is freely available and open source at GitHub (www.github.com/kelwyres/Bactabolize) under a GNU General Public License v3.0. Users should additionally cite COBRApy (38) if Bactabolize is used.

Klebsiella pneumoniae Species Complex-pan metabolic model

The 37 metabolic models from a previous study (14) were combined with the iY1228 model using the create_master_model.py script (available at 10.6084/m9.figshare.21728717). Briefly, all GPRs from the iYL1228 model and the associated sequences were included, as well as new GPRs identified from the 36 additional strains by manual curation following comparison to the matched phenotype data (as described in (14)). Additionally, orthologous sequence variants with <75% nucleotide identity to gene sequences associated with these gene reaction rules (GPRs) were added if there was phenotype data supporting the reaction. The biomass reaction was updated, removing the metabolites udpgalur_c and udpgal_c as their production was strain-specific.

Metadata annotations were improved using the improve_model_annotations.py script (also available in the Bactabolize code repository) resulting in the KpSC_pan v1 used in this study, available at www.github.com/kelwyres/KpSC-pan-metabolic-model.

Draft model generation

The annotated and unannotated genome of K. pneumoniae KPPR1 were obtained from Genbank accession number: CP009208.

Bactabolize draft models were generated using the draft_model command in Bactabolize v1 with the KpSC-pan v1 model as a reference, the annotated K. pneumoniae KPPR1 as input and the following options:

--min_coverage 25 --min_pident 80 --media M9 --atmosphere aerobic

A draft model for K. pneumoniae KPPR1 was also generated via CarveMe version 1.5.1 using the universal reference, the annotated K. pneumoniae KPPR1 as input, with the following commands:

‘-g M9 -i M9’. Subsequently, the --universe-file mode was also used, so the KpSC-pan model could be used as a reference, with the previously described command.

A draft model was generated using gapseq version 1.2 with the ‘doall’ command using the unannotated genome (as gapseq does not take annotated input files). Gap-filling was subsequently performed using the ‘fill’ command and a custom M9 media file to match the nutrient list found in Bactabolize (https://github.com/kelwyres/Bactabolize/blob/main/data/media_definitions/m9_media.json).

Finally, a draft model was constructed using the annotated genbank K. pneumoniae KPPR1 file and the KBase narrative (15) (https://narrative.kbase.us/narrative/ws.14145.obj.1).

Speed calculations

Modelling methods were timed via a script using the date +%s.%N command run before and after command on the MASSIVE computing cluster (Intel Xeon Gold 6150 CPU @ 2.70GHz and 155 GB of memory, CentOS Linux release 7.9.2009 environment). 10 individual complete KpSC genomes used in the quality control framework were tested for each method and the mean and range reported using R version 4.0.3 (70).

Performance comparisons

The annotated and unannotated genome of K. pneumoniae KPPR1 were obtained from Genbank under the accession: CP009208, and draft metabolic models were generated using Bactabolize, CarveMe, ModelSEED and gapseq as described above. The previously described, manually curated model for KPPR1 (iKp1289) was also included for comparison (15). MEMOTE version 0.13.0 was used to collate basic model statistics. The following KPPR1 phenotype data were retrieved from published studies: BIOLOG Phenotypic Microarray data (15) and single gene knockout data inferred from the outputs of a TraDIS transposon mutagenesis library (58).

A list of BIOLOG growth substrates for plates PM1, PM2A, PM3B, and PM4A (71) were converted where possible to BiGG and SEED IDs by manual search of the BiGG (bigg.ucsd.edu) and SEED websites (https://modelseed.org/biochem/compounds). An updated BiGG to SEED dictionary can be found in Data S1. A total of 143 of 190 carbon, 82 of 95 nitrogen, 46 of 59 phosphor and 26 of 35 sulfur substrates were successfully matched to BiGG and SEED IDs (Data S1). These growth data were compared to in silico predictions generated via FBA using the fba command from Bactabolize to optimise the biomass objective function with the following options:

--fba_spec_name m9 --fba_open_value -20

Gene essentiality was inferred from single gene knockout growth predictions using the sgk command from Bactabolize with the following options to mirror the growth conditions of the TraDIS library (LB media grown aerobically):

--media_type lb --atmosphere aerobic

In all cases, an objective value cut-off of ≥10^-4 was used to indicate binarised growth as per previous studies (14, 72).

In silico predictions were compared to matched phenotype data and the following accuracy metrics were calculated:

Model metabolite and reaction IDs were harmonized for overlap comparisons using the ‘Import model’ function from MetaNetX.org (59) after import and export via the write_sbml_model function from COBRApy.

Quality control framework

Illumina read sets (250 bp paired end) and completed genome sequences for 37 KpSC isolates were described previously (14). Here we randomly subsampled the Illumina reads at various depths (10 – 100, by increments of 10) using rasusa version 0.3.0 (73) in technical triplicate. Reads were then trimmed using TrimGalore version 0.5.0 (74) and assembled de novo with Unicycler version 0.4.7 (75), default parameters. Assembly statistics and assembly graph dead ends were calculated using the GFA-dead-end-counter version 1.0.0 (https://github.com/rrwick/GFA-dead-end-counter) (76). Draft metabolic models were generated with Bactabolize using the KpSC-pan v1 reference, and growth substrate profiles were predicted as described above. We compared the outputs from models generated for draft genome assemblies to those generated for the corresponding completed genomes. Where necessary models were gap-filled via the patch_model command.

Predictive accuracy of draft models

Novel growth phenotype data were generated for 10 KpSC clinical isolates from our in house collection using the VITEK 2 GN ID card system as described previously (14). Briefly, isolates were grown on Tryptic Soy (OXOID) agar plates overnight at 37°C, then analysed using VITEK 2 GN ID cards (bioMérieux) and read on the VITEK 2 Compact (bioMérieux) as per manufacturer’s instructions using software version 8.0. DNA was extracted for whole-genome sequencing via Genfind v3 extraction kit, library preparation performed using Nextera Flex (Illumina) using ¼ reagents. Paired-end read data (300 bp) were generated on an Illumina NovaSeq6000 SP v1.0 and have been deposited in the European Nucleotide Archive under Bioproject PRJNA777643 (individual read accession numbers are given in Data S4). Draft genome assemblies were generated with Unicycler, and draft metabolic models and growth predictions were generated with Bactabolize as described above.

Statistics and visualisation

Statistical analysis and graphical visualisation were performed using R version 4.0.3 (70), RStudio version 1.3.1093 (77), with the following software packages: tidyverse version 1.3.1 (78), viridis version 0.5.1 (79), RColorBrewer version 1.1-2 (80), ggpubr version 0.4.0 (81) ggpmisc version 0.4.4 (82), aplot version 0.1.6 (83), colorspace version 2.0-2 (84), ggpattern version 0.4.3-3 (85), ggtext version 0.1.1 (86) and glue version 1.4.2 (87).

Linear regression analysis was performed in R using the lm function and a third degree polynomial model was fitted to plots with the following equation: y ∼ poly(x, 3, raw = TRUE). The segmented linear model was fitted using segmented version 1.6-2 (88).

All code used to generate results can be found as supplemental material, https://github.com/kelwyres/Bactabolize, https://github.com/kelwyres/KpSC-pan-metabolic-model and on Figshare (10.6084/m9.figshare.21728717).

Logo

The Bactabolize logo was constructed in Inkscape version 1.0.1 (89). The font used is Proportional TFB (90) and Element (91).

Author contributions

Conceptualization: BV, JH, JMM, KEH, KLW

Methodology: SCW, BV, LMJ, JH, JMM, KLW

Software: SCW, BV

Validation: BV, HBC, KLW

Formal analysis: BV, KLW

Investigation: BV, KLW

Resources: AJ

Writing - Original Draft: BV, KLW

Writing - Review & Editing: All authors

Visualization: BV, KLW

Supervision: JMM, KEH, KLW

Project administration: KLW

Funding acquisition: JMM, KEH, KLW

All authors contributed to and approve of the manuscript in its current form.

Acknowledgements

We thank Sylvain Brisse for the isolates used in this study and review of the manuscript.

Funding

This work was funded by the Australian Research Council Discovery Project DP200103364.

KLW is supported by the National Health and Medical Research Council of Australia (Investigator Grant APP1176192).

Data availability

All genomes used are publicly available and generated by referenced studies. All Bactabolize code and data generated for this study can be found at https://github.com/kelwyres/Bactabolize. The KpSC-pan v1 reference model is freely available at https://github.com/kelwyres/KpSC-pan-metabolic-model.

Conflict of Interest

The authors declare that they have no competing or conflicts of interest.

Supplementary figures

Figure S1

Figure S2

Figure S3

Figure S4

Figure S5

Figure S6

Figure S7

Figure S8

Figure S9