Transporter genes in biosynthetic gene clusters predict metabolite characteristics and siderophore activity

Genome mining of transporters associated with biosynthetic gene clusters

Using two compiled sets of transporter-specific HMMs (Pfam and CATH) (https://www.cathdb.info/), we cataloged all classes of transporters across the MIBiG 2.0 database of characterized and experimentally validated biosynthetic gene clusters. We found that 56% of the bacterial BGCs in MIBiG contained at least one Pfam transporter hit and an additional 6% contained a CATH transporter hit without a Pfam domain (Fig. 1C). These percentages increased among BGCs that produce antibiotics (71%) and siderophores (78%), indicating that BGCs with these activities are more likely to contain at least one transporter. BGCs with transporters contained 2.5 transporter-associated domains across transport-annotated genes on average (Fig. 1D), which is expected because many ATP-dependent transporter systems have at least two domain complexes.

However, some BGCs contain considerably more transporter ORFs and domains, indicating that sometimes multiple transport systems can be associated with one BGC, although the number of protein domains that function as one transport system can often vary. The number of transporters in a BGC had no association with the number of metabolite structures reported for that BGC. The ATP-binding ABC transporter (ABC_tran) domain and the Major Facilitator Superfamily 1 (MFS_1) domain were the two most common transporter domains found in BGCs (Fig. 1B). A variety of proteins had nucleotide binding domains along with several different transmembrane domains—ABC_membrane and ABC2_membrane domains were most common but ABC2_membrane_2, -_3, and -_4 domains were also represented.

Examining domains specific for export, the ABC_membrane domain is often characteristic of exporters (e.g., Sav1866) (Velamakanni et al. 2008), but recently has been reported in the genes for siderophore uptake (YbtPQ) in Yersinia (Wang et al. 2020) and is therefore not necessarily indicative of export or import alone. The second most common transmembrane domain, ABC2_membrane, has been observed in the O-antigen polysaccharide exporter (Bi et al. 2018). The MacB-FtsX tripartite efflux pump was found in 60 BGCs, whereas the RND (ACR_tran) efflux pump was less common and found in only 10 BGCs. Other known efflux systems, such as SMR, MatE, and the MFS families 2 through 5 were comparatively rare across BGCs.

We next calculated co-occurrence correlations between all transporter protein families across BGCs and observed a strong negative correlation between MFS transporters and ATP-dependent transporters relying on the nucleotide binding domain, and a weaker negative correlation of MatE domains from the ATP-dependent NBD (Fig. 2A). This points toward a dichotomous choice between ATP-dependent and ATP-independent transport associated with a BGC.

Figure 2.

Presence of importer-specific domains and co-occurrence between transporters across BGCs. (A) Spearman's correlations between commonly occurring Pfam transporter domains across MIBiG BGCs—only correlations with P < 0.001 are shown. (B) Counts of Pfam transporter substrate binding domain families and corresponding substrate binding protein clusters described by Berntsson et al. (2010). (C) Counts of importer-specific CATH domains across MIBiG BGCs. The CATH functional family “Iron ABC permease” is essentially synonymous with the FecCD Pfam.

Open in new tabDownload PowerPointLink to figure

Multiple lines of annotation evidence indicated that many of the transporter genes associated with BGCs were likely to be importers. Importers can be involved in the uptake or re-uptake of molecules like siderophores and may also play roles in importing precursor metabolites for a BGC. The membrane domains specific to Type I importers (BPD_transp_1) and Type II importers (FecCD) were the most often observed ATP-dependent transmembrane domains besides ABC_membrane and ABC2_membrane (Fig. 1B). CATH Protein Structure Classification database HMMs (CATHs) that were specific for importer families in the Transporter Classification Database and found in BGCs included permeases for sugars, oligopeptides, and iron siderophores (Fig. 2C). Eleven percent of BGCs with a transporter also contained a substrate binding protein. Among the substrate binding proteins we searched for, the most common contained domain was Peripla_BP_2, also found in the E. coli B12 importer complex BtuCDF, and variants of this SBP (cluster A-II) are specific for siderophores and cobalamin (Berntsson et al. 2010). We also observed many substrate binding proteins with specificities predicted to include carbohydrates, oligopeptides, and peptide uptake (Fig. 2B; Berntsson et al. 2010). For example, the gene family SBP_bac_1 (SBP cluster D-I) (Berntsson et al. 2010) is specific for uptake of sugars and was found in the BGCs for the glycopeptide Mannopeptimycin and the aminoglycoside spectinomycin and may play a role in sugar precursor uptake (Supplemental Fig. S1). It was also found in the acarbose and acarviostatin BGCs (Supplemental Fig. S1), consistent with their putative roles as carbophors (Guo et al. 2012).

In the co-occurrence data, we observed pairing of different substrate binding proteins with different transmembrane domains. Peripla_BP_2 positively correlated strongly with FecCD and TonB_dep_Rec, genes known to be involved in siderophore uptake. BPD_transp_1 co-occurred with either SBP_bac_5 (SBP cluster C) or SBP_bac_1, whereas BPD_transp_2 co-occurred with Peripla_BP_4 (SBP cluster B-I). Taken together, these results show a logical organization of importer-specific transporter domains within BGCs that may be involved in either siderophore uptake, precursor uptake, or other roles. Regardless of the substrate specificity of these proteins, care must be taken when assuming that a transporter in a BGC is definitively for export of the matured product.

Prediction of siderophore and antibacterial activity from biosynthetic transporters

Because transporters are required for the ecological functions of biosynthesized specialized metabolites, we used machine learning to test if transporter classes were predictive of BGC-synthesized metabolite structures and functions. We noticed that metabolite activity labels in MIBiG were strongly associated with phylogeny: 83% of antibiotic BGCs were from Gram-positive bacteria, but only 40% of siderophore BGCs were from Gram-positive bacteria in the data set of curated MIBiG BGCs. To reduce the impact of this potential bias, we created separate training and testing data sets for activity prediction for Gram-positive and Gram-negative organisms.

Using our curated set of BGCs with transporters from the MIBiG 2.0 database, we tested for associations between BGC transport genes and metabolite function. We generated two activity classification tasks: (1) distinguishing siderophores (including known ionophores) (n = 16 Gram-positive, n = 24 Gram-negative) from non-siderophores (n = 142 Gram-positive, n = 42 Gram-negative); and (2) distinguishing antibiotics and antifungals (n = 131 Gram-positive, n = 27 Gram-negative) from non-antibiotics (n = 57 Gram-positive, n = 37 Gram-negative). We observed several statistically significant (Fisher's exact test; Q < 0.05) associations in the distribution of transporter types between both the siderophore and other activity classes. Among Gram-positive bacteria, 60% of siderophore BGCs contained the SBP Peripla_BP_2 and 55% contained the FecCD importer, but no BGCs with other activities had either (Fig. 3A; Supplemental Fig. S2). The situation was similar for Gram-negative siderophore BGCs. The TonB-dependent receptor (completely absent from Gram-positive bacteria) was the strongest signal, found in almost 80% of Gram-negative siderophore BGCs with a transporter and never in BGCs with other activities.

Figure 3.

Transporter domains are predictive of siderophore BGCs. (A) The frequencies of common transporter Pfam domains across siderophore BGCs and BGCs of other known activities in Gram-positive and Gram-negative bacteria. Bars in green were significantly different in frequency between the two classes (Fisher's exact test; Q < 0.05) (B) Precision-recall curves for two-layer decision trees classifying siderophore BGCs using Pfam transporter, CATH transporter, and Pfam biosynthetic gene features in Gram-negative and Gram-positive bacteria. (C) Examples of three siderophore BGCs without activity labels in MIBiG 2.0, which could be identified using transporter frequencies. Transporter genes are blue, core biosynthetic genes (NRPS and PKS) are dark red, accessory biosynthetic genes are light red, and regulatory genes are green.

Open in new tabDownload PowerPointLink to figure

To assess siderophore predictability from BGC gene content, we used decision trees with only two layers applied to different feature sets of protein domain annotations, transport-affiliated Pfams, transport-affiliated CATH HMMs, and biosynthetic Pfams. To avoid issues with class imbalance, we report precision and recall on the siderophore class, because siderophore prediction requires searching for a minority class (siderophores) within a background of mostly non-siderophores. With the transport-only features, we found that just with two gene decisions, it is possible to achieve 100% precision with <80% recall using either Pfam or CATH transporter annotations for Gram-negative siderophores, and 100% precision with <80% recall using CATH transporter annotations for Gram-positive siderophores (Fig. 3B; Supplemental Table S5). On the other hand, when using all biosynthetic annotations within BGCs, we found that two-layer decision trees trained on biosynthetic genes performed substantially worse at predicting siderophore activity than those trained on transporter genes (Fig. 3B; Supplemental Table S5). We further validated our results by training LASSO linearized regression models, which do not model interactions between features. These models obtained a slightly improved area under the precision-recall curve (AUPRC), indicating that very simple transporter patterns are highly predictive of whether a BGC is siderophore producing or not in our data set (Fig. 3B). Transporter features predictive of siderophores were consistently selected by LASSO across stratified cross-validation repeats, giving evidence that these patterns are robust (Supplemental Fig. S3). The top predictive biosynthetic features of siderophores were the IucA/IucC protein family (used by antiSMASH to label siderophores and which is known to be involved in aerobactin biosynthesis) and condensation domains (likely to capture nonribosomal peptide siderophores), but predictive effect sizes were smaller than those for transporters (Supplemental Fig. S3).

Using siderophore-specific transporter genes, we attempted to predict siderophore classes for any remaining gene clusters that have no annotated function in MIBiG2 (that we had not already hand curated). We searched for gene clusters containing the siderophore-predictive genes FecCD and Peripla_BP_2 and found six additional BGCs with no annotated activity, three of which were experimentally validated by the literature to be siderophores (Fig. 3C; Carran et al. 2001; Matsuo et al. 2011; Chen et al. 2013). The remaining three identified BGCs were false positives. One of them, the BGC for the antibiotic Ficellomycin, only contained these transport genes in flanking regions not known to be involved in biosynthesis (Liu et al. 2017), whereas the herbimycin A BGC contains the transporters genes in the reverse reading frame from the BGC, separated by an unusual 16-kb intergenic region and the genes appear to be fragmented (Rascher et al. 2005). The other false positive was Lividomycin, an antibiotic that does seem to be a rare non-siderophore BGC with FecCD and Peripla_BP_2 transporters in the MIBiG2.0 database.

To understand the extent of transporter specificity for particular classes of BGCs, we expanded our analysis to 95,293 BGCs in the antiSMASH database, which is a set of predicted BGCs in microbial genomes from the RefSeq database. Frequencies of general transporters and transporter genes that were specific for siderophore biosynthesis in MIBiG (FecCD and the TonB-dependent Receptor) were calculated across all antiSMASH BGCs by bacterial genus (Fig. 4; Supplemental Fig. S4). The general ABC transporter ATP-binding domain and MFS superfamily transporter genes varied in their frequencies across different classes of BGCs, but there were some consistent patterns. Thiopeptides, NRPS-independent siderophores, and arylpolyene BGCs tended to have an MFS exporter, but ribosomally synthesized products (e.g., Lasso peptides and Lantipeptides) consistently had ATP-dependent transport mechanisms (Fig. 4A).

Figure 4.

Presence of general and siderophore-specific transporters by biosynthetic class and bacterial genus across the antiSMASH database. (A) ATP-dependent and ATP-independent (MFS) transporters are commonly associated with a variety of BGCs in the antiSMASH database across a wide range of genera. Each point is the percentage of a BGC class with a transporter within a particular genus. Each genus is colored by its Gram status, and genera with fewer than 20 BGCs of a particular class are excluded. (B) Siderophore-specific transporters are associated with few BGC classes in the antiSMASH database. Each point is the percentage of a BGC class with a transporter within a particular genus. Each genus is colored by its Gram status, and genera with fewer than 20 BGCs of a particular class are excluded.

Open in new tabDownload PowerPointLink to figure

Alternatively, the TonB-dependent receptor and the FecCD Type II importer were highly restricted to specific classes of BGCs. They only consistently appeared in NRPS-independent siderophore and NRPS clusters and were nearly absent from entire other classes of BGCs (Fig. 4B). This is consistent with the known NRPS-independent and NRPS biosynthetic pathways for siderophores. The TonB-dependent receptor was also found associated with Lasso peptides (possibly functioning as a resistance gene) (Mathavan et al. 2014), but FecCD was not. Twenty-one percent of NRPS clusters contained a FecCD gene, and 28% contained a TonB-dependent receptor, possibly indicating that at least one in five uncharacterized NRPSs may function as siderophores. Thus, even in an uncurated data set of thousands of BGCs, it appears as though siderophore-specific transporters only rarely occur in biosynthetic gene clusters that are unlikely to have siderophore functions.

Classifying antibiotics and antifungals from either transporters or biosynthetic genes proved more challenging than classifying siderophores (Supplemental Table S5; Supplemental Fig. S5). In Gram-negative bacteria, we observed positive associations between the MacB-FtsX tripartite efflux pump with antibacterial or antifungal activity (Fisher's exact test; Q < 0.05). MacB and associated components (FtsX and OEP) were positively associated with antibiotic activity by LASSO logistic regression. This result was stable across both cross-validation folds and repeats of cross-validation (Supplemental Fig. S5). We found that seven of 27 Gram-negative antibacterial BGCs contained a MacB, and no BGCs in our classes of other activities contained MacB. Although MacB is involved in export of the siderophore Pyoverdine (for which there is not an accurate BGC in MIBiG) in Pseudomonas (Greene et al. 2018a), in general, MacB may be a strong indicator of antibacterial activity for a BGC. Previously, we identified a number of MacB-FtsX exporters in BGCs from novel Acidobacteria (Crits-Christoph et al. 2018), possibly indicating a role in antibacterial activity for these BGCs. LASSO effect sizes for individual biosynthetic genes were substantially lower (Supplemental Fig. S5).

Association of biosynthetic transporter classes with the molecular weights and lipophilicity of their putative substrates

We next hypothesized that transporter classes could be predictive of other molecular features beyond functional activity. There was no strong correlation between the molecular size of the metabolite produced and Gram status of the bacteria encoding each corresponding BGC in the MIBiG data set. Thus, we tested for differences in transporter classes in BGCs producing metabolites that were (1) less than and (2) greater than 1000 Da in size across all bacteria. There was a significant difference in the frequencies of some transporters between BGCs with different metabolite molecular weights (Fig. 5A). The strongest difference was in the distribution of MFS transporters—found in 57% of BGCs with products <1000 Da, but only 14% of those >1000 Da (Fisher's exact test; Q < 0.05). The 95th percentile of metabolite molecular weights for clusters with MFS_1 and without ABC_tran was 1082 Da, which may be approaching a biological limit for the molecular weights of substrates for these transporters. Conversely, the ATP-dependent ABC_tran domain was found in 89% of BGCs producing high molecular weight compounds but only 42% of those producing low molecular weight compounds (Fisher's exact test; Q < 0.05). MacB/FtsX and the rarer transmembrane domains were also associated with higher molecular weight compounds.

Figure 5.

Transporter domains associated with molecular size and partition coefficient. (A) The frequencies of common transporter Pfam domains in BGCs that synthesize metabolites >1000 Da (left) and <1000 Da (right). Bars in green were significantly different in frequency between the two classes (Fisher's exact test; Q < 0.05). (B) Precision-recall curves for two-layer decision trees and LASSO logistic regression models classifying BGCs producing metabolites >1000 Da using Pfam transporter, CATH transporter, and Pfam biosynthetic gene features. (C) The distribution of metabolite molecular weights synthesized by BGCs with at least one NBD-binding ABC transporter domain, at least one MFS domain, and the ABC2_membrane_3 transmembrane domain. (D) Predicted partition coefficients (log P) for metabolites synthesized by BGCs that contain at least one variant of two different ABC transporter transmembrane domains.

Open in new tabDownload PowerPointLink to figure

We tested for an association within the two largest chemical classes in the data set, PKs and NRPs, and found that also within the PKS and NRPS biosynthetic classes ATP-dependent transporters were associated with larger metabolite molecular weights than the MFS family (Fig. 5C). We also observed that the ABC2_membrane_4 was associated almost exclusively with large RiPPs, with almost all of the BGCs in which it is found in producing compounds >1500 Da in size. After training both LASSO logistic regression and two-layer decision tree models to classify whether produced molecules are >1000 Da, we found that transporter genes were able to distinguish large from small metabolites with moderate precision and recall (Supplemental Table S5) and an AUPRC of up to 42%, and biosynthetic genes performed similarly (Fig. 5B). Top transporter features consistently had larger effect sizes than top biosynthetic features, indicating that transporter-based features provided clearer signals. This result was stable across both cross-validation folds and repeats (Supplemental Fig. S6). The biosynthetic protein family most associated with high molecular weight metabolites was Glycos_transf_2, likely owing to the addition of sugar groups to metabolites by these enzymes in BGCs.

It has previously been reported that transporters can be specific for compounds with similar hydrophilicity (Rempel et al. 2020). The lipophilicity of a metabolite is often considered critical for its success in clinical development for human therapeutics (Arnott and Planey 2012). With LASSO logistic regression, we predicted partition coefficients (log P), a measure of lipophilicity, for all of the metabolites in MIBiG, and tested how well metabolite partition coefficients could be predicted by gene content. We found that the presence of five transporter classes was significantly associated with increased lipophilicity (Fisher's exact test; Q < 0.05). In particular, we observed an association between varying ATP-dependent transmembrane domains and log P, with ABC2_ membrane_3 domain co-occurring with BGC-metabolites with a high log P (median 4.4) (Fig. 5D) and ABC2_membrane_4 domain co-occurring with BGC-metabolites with a low log P (median −6.2). Although the ABC2_membrane_4 association is likely a result of its exclusive association with large RiPP products, the ABC2_membrane_3 domain occurred in multiple BGC classes, mostly polyketides, and was still associated with a decrease in log P just within the polyketide class. LASSO logistic regression distinguished log P > 0 from log P < 0 with 77% AUPRC (Supplemental Table S5; Supplemental Fig. S7). On this task biosynthetic genes were distinctly superior over transporters at prediction of log P, obtaining a 83% AUPRC with a LASSO logistic regression trained on biosynthetic genes (Supplemental Table S5).

Identifying novel siderophore-like biosynthetic gene clusters in the human microbiome

To show the predictive utility of BGC-associated transporters, we mined BGCs with siderophore-specific transporters in metagenomic genomes (dereplicated per species) from (1) the gut microbiomes of neonatal infants in the intensive care unit (Olm et al. 2019), and (2) a cross-study collation of genomes assembled from multiple human gut studies (Nayfach et al. 2019). Identified “siderophore-like” BGCs putatively produce siderophores, as they contained the transporter classes that can achieve near 100% siderophore specificity in the MIBiG database: (1) Peripla_BP_2 and FecCD in Gram-positive bacteria, and (2) Peripla_BP_2, FecCD, and TonB_dep_Rec in Gram-negative bacteria. We identified 1442 BGCs with siderophore-like transporter classes (Fig. 6A; Supplemental Table S6) and then grouped them into novel gene cluster families using BiG-SCAPE, resulting in 75 siderophore-like gene cluster families (Fig. 6B).

Figure 6.

BGCs with siderophore-like transporters from human gut microbiomes. (A) Concatenated ribosomal protein tree (collapsed to the genus level) for high-quality genomes from the infant and adult gut microbiomes that encode siderophore-like BGCs. On the right are counts of siderophore BGCs from infant gut genomes (blue) and adult gut genomes (red). (B) Gene Cluster Families of BGCs containing known siderophores (bright red) and human microbiome–derived BGCs with siderophore-like transporters. BGCs are connected by similarity to other BGCs in the same gene cluster family, calculated using BiG-SCAPE. (C) Families of siderophore-like BGCs without any similarity to existing known BGCs. BGCs in the network are colored by the taxonomy of the genome of origin and are grouped and labeled by the antiSMASH reported biosynthetic class: for NRPS gene clusters, the adenylation domain specificities are reported.

Open in new tabDownload PowerPointLink to figure

Most siderophore-like BGCs were in large gene cluster families with other BGCs that also contained the same set of transporter hits. Twenty-three percent of microbial genomes from the neonatal infant gut microbiomes had siderophore-like BGCs, but only 3% of those assembled from adult gut microbiomes did. Siderophore-like BGCs were identified across a range of bacterial genera, but the majority were from the Staphylococcus or Enterobacteriaceae, which are known to be in high abundance in the neonatal gut microbiome. The genera with the most siderophore-like BGCs were Staphylococcus and Klebsiella, common hospital-acquired pathogens of neonates. Five of the 75 identified siderophore-like BGC families included a known representative gene cluster in MIBiG, and four of the known representatives were siderophores, again pointing to the specificity of these transporter classes.

The rest of the siderophore-like gene cluster families that were identified had no closely characterized representative in MIBiG, indicating that there is likely capacity for production of multiple novel siderophores in the human gut microbiome (Fig. 6C). Almost all siderophore-like BGCs were either NRPS or NRPS-independent siderophore classes, the latter of which is based on the presence of the IuA/IuC gene family, as in aerobactin biosynthesis. Of the NRPS siderophore-like gene clusters, the majority had adenylation domain specificities for serine (ser) and 2,3-dihydroxybenzoate (dhb), indicating similar catechol-containing nonribosomal biosynthetic pathways to siderophores like enterobactin and salmochelin (Supplemental Table S7). We observed substantial genetic diversity between gene cluster families containing similar NRPS domains, which may indicate the existence of possible unknown derivatives of these siderophores in the human microbiome. In adult gut microbiome samples, one large novel NRPS siderophore-like gene cluster with unknown adenylation specificity was identified in Coprococcus, members of the Lachnospiraceae, often considered to be important commensals in the human gut (Chen et al. 2017; Duvallet et al. 2017). Thus, although a majority of siderophore-like BGCs in the human microbiome contained core enzymes similar to known siderophore biosynthetic pathways, there was substantial genetic diversity that could indicate further unexplored structural variation.

Curation and selection of BGCs and transporter annotations

We parsed the MIBiG 2.0 database of biosynthetic gene clusters metadata and extracted information including host genus, compound count, chemical structures of the metabolite product, known metabolite activities, and the number of open reading frames for each BGC. Using Entrez and NCBI, we assigned the expected Gram status for each BGC based on phylum, coded 0 = Gram-negative, 1 = Gram-positive, 2 = Fungal, 3 = other. For the purpose of this manuscript, we only analyzed BGCs from Gram-positive and Gram-negative bacteria. We noticed that the activity labels in MIBiG 2.0 were often incomplete, and manually added a set of antibacterial, siderophore, and antifungal labels derived from the literature (Supplemental Table S1). We found that 28 BGCs (1.8%) in MIBiG were unusually large in length, and comparisons to published papers on these BGCs showed that their MIBiG counterparts were overextended in comparison to the validated BGC. For this reason, we eliminated the 28 BGCs over 60 ORFs in length. Using Python and RDKit, we calculated molecular weights and partition coefficients (log P) using the algorithm described in Wildman and Crippen (1999) for all 1042 MIBiG BGCs with a single associated compound structure. Two hundred thirty-eight BGCs have more than one associated metabolite structure, and these multistructure BGCs were not used in our structural association analyses. We annotated biosynthetic genes in MIBiG with HMMER hmmsearch (Eddy 1998) and cath-resolve-hits (Lewis et al. 2019) on a set of the 99 most commonly represented biosynthetic Pfams in antiSMASH BGCs obtained from Cimermancic et al. (2014) to generate a counts table of biosynthetic protein families for each BGC (Supplemental Table S4).

To obtain a comprehensive overview of the distribution of transport-associated protein domains in biosynthetic gene clusters, we generated two separate feature tables: (1) using CATH (Sillitoe et al. 2019) HMMs, and (2) using Pfam (El-Gebali et al. 2019) HMMs. For the first, we downloaded all proteins in the Transporter Classification Database (TCDB) (Saier et al. 2016) and annotated them with hmmsearch and cath-resolve-hits, using all CATH Functional Family HMMs. We then selected all CATH HMMs that were represented at least five times. We then manually curated this list down to 180 final CATH HMMs that were transport specific. We then calculated the specificities of each CATH HMM for TCDB families; 80 were specific for exactly one TCDB family.

For the second set of features, we took all protein sequences with an annotation including “Transport” in the antiSMASH database v2 (Blin et al. 2019a) and annotated these proteins with the Pfam-A set of HMMs using hmmsearch and the option‐‐cut_ga. We then selected highly represented HMMs and manually curated this list to be transporter specific and representative of the major transporter classes in TCDB. We then also selected the Pfam Substrate Binding Protein and Periplasmic Binding Protein HMMs that were represented more than five times in MIBiG. When comparing both the Pfam and CATH set of HMMs, we found substantial overlap, but the CATH set is composed of 166 domain features, whereas the Pfam set only contains 18 (Supplemental Tables S2, S3).

Machine learning to predict metabolite structural and functional characteristics

To identify associations between metabolite functional classes and structural properties with BGC gene content, we used traditional statistical tests and different machine learning models. Our classification tasks were (1) siderophore (n = 16 Gram-positive; n = 24 Gram-negative) versus other activity (n = 142 Gram-positive; n = 42 Gram-negative), (2) antibiotic and antifungal (n = 131 Gram-positive; n = 27 Gram-negative) versus other activity (n = 57 Gram-positive; n = 37 Gram-negative), (3) metabolite molecular weight > 1000 Daltons (n = 149) versus metabolite molecular weight < 1000 Daltons (n = 421), and (4) predicted partition coefficient log P < 0 (n = 220) versus log P ≥ 0 (n = 350). For the functional classification tasks, we noticed a strong class imbalance with Gram status, so we performed functional classification separately for BGCs from Gram-positive and Gram-negative bacteria. We first tested for univariate differences in proportional representation of each transporter BGC between classes for classification tests using Fisher's exact test in Python (Q < 0.05, Benjamini–Hochberg correction) (Hochberg and Benjamini 1990).

We then assessed the predictive power of the three sets of features for each BGC: transporter Pfam HMMs, transporter CATH HMMs, and biosynthetic Pfam HMMs. Features were counts of protein families, which were standardized using the StandardScaler function in the scikit-learn package. Given the nature of our study, we used simple models to ensure reliability and interpretability of our results. We fit two classes of machine learning models: (1) LASSO-penalized logistic regression, which fits a linear model with a sparsity penalty on weights; and (2) shallow decision trees (of depth one or two), which can classify based on splitting at most two features (Franklin 2005). All models were trained using the Python package scikit-learn. Because of data size and class imbalance, we fit models using repeated, stratified k-fold cross-validation (“RepeatedStratifiedKFold” in scikit-learn) with five repeats and five folds. On each cross-validation split of our data, we computed the area under the precision-recall curve (AUPRC) to evaluate performance for our class-imbalanced tasks. Thus for the final output of this procedure, we reported the mean of accuracy, precision, recall, and AUPRC each generated from five repeats of different random fivefold partitions of the data. We further use these repeated cross-validation splits to see which features are consistently used for classification across repeats.

Annotating metagenomic siderophore-like BGCs from the human microbiome

To assess the distribution of siderophore-like BGCs in the human microbiome, we downloaded two sets of genomes assembled from metagenomes obtained from the human gut microbiome: 2425 genomes from a neonatal intensive care unit (NICU) premature infant microbiome (Olm et al. 2019) and 24,345 genomes from a diverse set of mostly adult human cohorts (Nayfach et al. 2019). We ran antiSMASH 5.0 on these genomes and then scanned predicted BGCs for at least two of the Pfams that were found to be specific for siderophore BGCs (FecCD, Peripla_BP_2, and Ton_dep_Rec). We dereplicated these BGCs and compared them to known MIBiG siderophore BGCs using the software BiG-SCAPE (Navarro-Muñoz et al. 2020) run with default settings. We then considered BGCs with either set of hits and reported their genomic taxonomic distribution based on the closest BLAST hit representatives of genomic ribosomal proteins to taxonomic genera defined by GTDB (minimum percent identity of hits >80%) (Parks et al. 2018).