Predicting precursors of plant specialized metabolites using DeepMol automated machine learning

2.1 Data

This study employed a manually curated dataset consisting of 566 alkaloids and their biosynthetic precursors, sourced from scientific literature and the Kyoto Encyclopedia of Genes and Genomes (KEGG), as described by Eguchi et al. [13]. The dataset features 17 precursors, including 10 amino acids, one aromatic acid (anthranilate), five terpenoids, and indole-3-glycerol phosphate (IGP).

To further expand this dataset, we extracted the biosynthetic pathways for terpenoids, phenylpropanoids, glucosinolates, and benzoxazinoids from KEGG. These pathways were converted into directed acyclic graphs, with compounds represented as nodes and reactions as edges, and we employed depth-first search to identify precursors.

For terpenoids, the identified precursors included GPP (precursor of monoterpenoids [18]), GGPP (precursors of diterpenoids [18] and carotenoids [19]), FPP (precursors of sesquiterpenoids and triterpenoids [18]), and campesterol (precursor of brassinosteroids [20]).

The precursors for phenylpropanoids included L-phenylalanine, the universal starter for all phenylpropanoids, p-Coumaroyl-CoA (precursors of flavonoids and stilbenoids [8]), naringenin and liquiritigenin (isoflavonoids), as well as apigenin and kaempferol (flavones and flavonols [21]), and pelargonin, cyanidin, and delphinidin (anthocyanins [21]).

The precursors for glucosinolates are amino acids [22], which are already included in the dataset, while the precursor for benzoxazinoids is IGP [23]. Figure 1 illustrates the number of compounds associated with each precursor.

Figure 1:

Number of compounds associated with each precursor. Note that one compound can have more than one precursor.

The dataset was divided into training (70 %), validation (20 %), and test (10 %) sets. Figure 2 shows clusters of compounds grouped by their chemical similarities using t-distributed Stochastic Neighbor Embedding (t-SNE) from a similarity matrix (Supplementary Material), positioning compounds with greater chemical similarity closer together, while those with significant differences are placed further apart. Figure 2A and B show the distribution of compounds by chemical class and dataset, respectively.

Figure 2:

Chemical space of the dataset visualised using t-SNE. (A) Natural product classes. (B) Training, validation and test sets. (C) A training and test set with compounds from distant similarity clusters. (D) A training and test set where compounds from different natural product classes overlap in similarity.

Additionally, the dataset was split into two alternative ways to evaluate model robustness under more challenging training-test conditions. First, we selected more distant compound clusters (Figure 2C), ensuring at least minimal precursor representation by including a few compounds in the training set; otherwise, the model could not predict that specific precursor. Second, compounds that are highly similar but belong to different NP classes were excluded from the training set and included in the test set, as illustrated in Figure 2D. This configuration creates a challenging prediction task because, despite the molecules’ high similarity, they originate from different precursors.

2.2 Metrics

The metrics used include the F1 score, macro F1 score (mF1), recall, macro recall (mRecall), precision, and macro precision (mPrecision), each offering insights into different aspects of model performance.

The mF1 score is calculated as the harmonic mean of precision and recall for each label individually, and these values are then averaged to produce a single score for the entire dataset. A similar approach is used for mRecall and mPrecision, where the averages of precision and recall across labels are computed. The exact formulations can be found below.

Where N represents the total number of labels, with Precision _i and Recall _i indicating the precision and recall for label i, respectively. True Positives _i are the true positive predictions for label i, and False Negatives _i are the missed predictions for label i. Finally, the False Positives _i represent the incorrect positive predictions.

2.3 Automated machine learning

Our experimental framework utilized the DeepMol AutoML engine (version 1.1.5) [14] to train all classifiers for multi-label classification available in scikit-learn version 1.2.0, employing molecular fingerprints (FPs) and descriptors as inputs on the training set, and evaluating their performance on the validation set. The models implemented include ridge classifiers, decision trees, random forests, and k-nearest neighbours, among others.

AutoML streamlines ML tasks by automating the search for the optimal pipeline, which integrates a feature set, feature selection, and classifier. The main goal during the optimization process was to maximize the mF1 score on the validation set. We employed the Tree-structured Parzen Estimator (TPE) algorithm over 500 trials to determine the most appropriate hyperparameters and methods. Upon completing these trials, we identified the best pipeline and subsequently retrained it using the combined training and validation data to assess its performance on the test set.

2.4 Fingerprints for comparison

To benchmark our method, we compared it against several established molecular fingerprints (FPs) that have shown promise for NP classification and bioactivity prediction.

The following fingerprints were considered:

1. MinHashed Fingerprints (MHFP) is a string-based FP algorithm that operates similarly to circular fingerprints, but instead of using atom identifiers, it uses SMILES substrings as fragment identifiers. These identifiers are then hashed and stored in a fixed-size vector using the MinHash algorithm, as shown by Probst et al. [24]. A recent study also highlighted MHFP’s effectiveness in NP bioactivity prediction [17].

2. Neural Natural Product Fingerprints (NeuralNPFP) is a neural network-based molecular representation specifically designed for NPs. This FP encodes implicit NP-relevant information while maintaining details about their chemical space by training neural networks to differentiate between natural and synthetic compounds [25].

3. NPClassifier Fingerprints are a variation of Morgan FPs designed for NP classification, where non-negative integers represent the count of repeated atomic substructures (or molecular fragments) rather than binary values. Morgan fingerprints represent molecular structures by encoding the substructures within a circular region around each atom into a fixed-length bit vector. These substructures are hashed into unique identifiers, with corresponding bits in the vector set to indicate the presence of specific atomic environments. The fingerprint’s size and detail are controlled by adjusting the radius, which defines the neighbourhood size around each atom. NPClassifier FPs combine three Morgan FPs, each with different radii (0, 1, 2), to capture structural information at varying distances [16].

This set of FPs highlights the diversity of approaches in molecular representation, from neural network-based embeddings to substructure counting, which provides a strong foundation for assessing the strengths of our method.

To ensure a fair comparison, we implemented each FP in DeepMol and used its optimization framework to identify the best combination of methods and models, while keeping the same fingerprint consistent across all trials. For each FP, we conducted 50 optimization trials using the TPE algorithm, following the standard AutoML approach in DeepMol.

2.5 Statistical methods

We employed the Wilcoxon Signed-Rank (WSR) Test to evaluate differences in the metric values between the models for all the test sets. In this context, the samples consist of the metric values obtained by two models for each label within the multi-label classification framework. The objective is to determine the significance of any difference in the median value of a metric across all labels. This hypothesis test is formulated as follows. Given metric values for two models across n labels, m _1i and m _2i , calculate the differences d _i = m _1i − m _2i for each label i. For these differences, ignore d _i = 0 and rank the absolute differences |d _i |, assign ranks R _i and compute W + = ∑ d i > 0 R i and W − = ∑ d i < 0 R i , the test statistic W is defined as W = min(W ⁺, W ⁻). A p-value is calculated as the probability of observing a value of W _ref , determined by a reference distribution under the null hypothesis, as extreme as or more extreme than the observed value (W). The null hypothesis is that there are no significant differences between the median metric values between the two models. We considered a p-value lower than 0.05 sufficient to reject the null hypothesis.

Additionally, we applied the same test to each pair of models for the experimental setup outlined in [13], where cross-validation is conducted solely for predicting the precursors of alkaloids. This hypothesis test is formulated as follows. Given two models evaluated across n labels and r folds, resulting in performance metrics m _Aij and m _Bij for models A and B respectively, for each label i and fold j, we perform the following steps: calculate the differences d _ij = m _Aij − m _Bij , rank the absolute differences |d _ij |, and apply the Wilcoxon Signed-Rank Test as explained above.

2.6 Precursor predictions for LOTUS-DB plant compounds

After selecting the best model, we retrained it with the collected data (train, validation and test sets). Next, we downloaded all the data from LOTUS-DB and filtered the compounds annotated as plant compounds. We then predicted the precursors for all these compounds using the retrained model and analyzed the main patterns of predictions by class of compounds, afterwards.

3.1 AutoML consistently outperform MGCNN and NP-focused fingerprints

Figure 3 provides a comparative analysis of several classification models across three metrics: macro F1-score (mF1), macro Precision (mPrecision), and macro Recall (mRecall).

Figure 3:

Results of each model evaluated on four different settings (see methods). (A) The test set originated from the stratified split. (B) Cross-validation datasets for predicting the precursors of alkaloids from the Eguchi et al. 2019 paper. (C) The distant similarity clustered compounds dataset. (D) the overlapping natural product classes dataset.

The models were selected based on the outcomes of an AutoML search, where Ridge classifiers using Morgan FPs and Layered FPs emerged as the top performers in terms of model and feature extraction methods, respectively. Among the models, the standout performer was a Ridge classifier trained on Layered FPs. The Layered FPs are molecular FPs that capture multiple structural attributes of a molecule using a multi-layered approach, where each layer captures different molecular features, including topology, bond order, atom types, ring presence and sizes, and aromaticity. The construction process involves generating a multi-layered hash for each molecular subgraph to set bits in the fingerprint, allowing for flexible representation. On top of this, the AutoML engine from DeepMol improved this model’s validation results by applying feature selection, specifically filtering out features with a variance below 0.0628, which was identified as the optimal threshold by the AutoML. The approach that combines Layered FPs, a variance threshold filter, and a Ridge classifier is called the Layered FP+Ridge Classifier. In addition, we evaluated a Ridge classifier trained with Morgan FPs, known for their straightforwardness in mapping molecular substructures directly to specific FP bits. Similar performance patterns were observed for NPClassifier FP and MHFP, where Ridge classifiers achieved the best results. However, for the Neural NPFP, a classifier chain model leveraging a Multi-Layer Perceptron (MLP) as the base estimator delivered the top performance.

Layered FP+Ridge Classifier delivered the best performance on the test dataset (Figure 3A), achieving an mF1 score of 0.898 and an mPrecision of 0.975. This classifier demonstrated a well-balanced trade-off between mPrecision and mRecall, with a slight bias towards mPrecision. Notably, this model exhibited significant improvements over the MGCNN model [13]. It achieved markedly higher F1 and precision scores than MGCNN, with statistical significance (p − value _F1 < 0.05, p − value _precision < 0.01, in a WSR test, see Supplementary Material and Methods). While other models, such as NPClassifierFP+Ridge Classifier, MHFP+Ridge Classifier, and Morgan FP+Ridge Classifier, also outperformed MGCNN in terms of mF1 score, the Layered FP+Ridge Classifier was the only one to achieve a statistically significant improvement. As depicted in Figure 3A, MGCNN yielded a marginally higher mRecall (0.871) than the other models, though this difference was statistically significant only in comparison to MHFP+Ridge Classifier (p − value _recall < 0.05, in a WSR test) and NeuralNPFP+MLPs (p − value _recall < 0.01, in a WSR test). The pipeline utilizing NeuralNPFP underperformed across all metrics when compared to the other approaches.

The Layered FP+Ridge Classifier with the Ridge classifier outperformed MGCNN in predicting nearly all precursors, except cholesterol, campesterol, and FPP. This indicates that MGCNN exhibits slightly stronger predictive capabilities for steroid, sesquiterpenoid, and triterpenoid precursors (for further details, refer to the Supplementary Material). When focusing specifically on the alkaloids dataset from Eguchi et al. [13] (Figure 3B), Layered FP+Ridge Classifier consistently outperformed the state-of-the-art MGCNN across all metrics, with statistically significant improvements in mF1 and mPrecision. Although the MHFP+Ridge Classifier achieved a higher mPrecision than the Layered FP+Ridge Classifier, this difference was not statistically significant, and the MHFP model underperformed in terms of mRecall and mF1. Furthermore, the Layered FP+Ridge Classifier was the only approach that surpassed MGCNN across all metrics.

These performance trends persisted with the clustered dataset (Figure 3C), where the Ridge classifiers consistently outperformed MGCNN across all metrics. However, the differences were not statistically significant, likely due to the small number of labels evaluated (13) and the high variability in label scores across different methods (see Supplementary Material). The only statistically significant differences were observed when comparing all Ridge classifiers to the NeuralNPFP+MLPs pipeline.

In the overlapping dataset (Figure 3D), the Morgan FP+Ridge Classifier and MGCNN models showed the best performance, although their advantage was slight and not statistically significant, except in precision scores, where the Morgan FP+Ridge Classifier significantly outperformed the Layered FP+Ridge Classifier. All models showed a substantial improvement over the NeuralNPFP+MLPs in terms of F1 scores. While other models exhibited lower precision scores, these differences were not statistically significant due to the large variability across labels (more details can be found in the Supplementary Material).

The ML pipelines discovered by DeepMol AutoML demonstrated clear advantages over models utilizing Natural Product-focused fingerprints and the state-of-the-art tool MGCNN. While the Natural Product-focused representations such as NPClassifierFP and MHFP had isolated strengths in certain metrics, they generally underperformed compared to the Ridge classifiers using Layered and Morgan FPs. These findings highlight the importance of feature selection and optimization, as facilitated by DeepMol’s AutoML framework, in achieving state-of-the-art performance in molecular classification tasks.

3.2 Predictions for LOTUS-DB compounds show consistent patterns

Next, we predicted the precursors for all plant compounds in LOTUS-DB to annotate them and provide comprehensive insights. We also conducted a large-scale analysis of the patterns in our predictions across all the 190,273 compounds.

The pie charts in Figure 4 offer insights into the predictions made by the Layered FP+Ridge Classifier for plant compounds in LOTUS-DB [15], a comprehensive database of NPs with taxonomic details and chemical class ontologies. This database enables the filtering of compounds based on their chemical class and the organisms responsible for their production, facilitating our further analysis.

Figure 4:

Predictions for LOTUS-DB using Layered FP+Ridge Classifier. The top pie chart illustrates the percentage of annotated and unannotated plant compounds in LOTUS-DB. A majority (77.6 %) of the compounds are annotated, while 22.4 % remain unannotated. The bottom four pie charts show the percentage of predictions for each precursor. Each chart represents a different subset of compound classes.

The top chart (Figure 4A) illustrates that our tool predicts the precursors of 77.6 % of plant compounds in LOTUS-DB, leaving 22.4 % unannotated. These unannotated compounds likely stem from primary metabolism or have precursors not included in the current dataset. While compounds from primary metabolism probably represent only a small portion of this group, the majority might consist of metabolites with precursors outside our current scope. For example, polyketides are synthesized from precursors like acetyl-CoA, malonyl-CoA, and propionyl-CoA, none of which were included in our dataset, as plant polyketide production was not a focus of this study. Expanding the dataset to include additional precursors, especially those involved in the biosynthesis of polyketides and other important plant metabolite classes, will be an important next step.

The bottom charts display the percentage of predicted precursors involved in the biosynthesis of alkaloids, terpenoids, phenylpropanoids, and glucosinolates. Benzoxazinoids were excluded from the analysis as this class of compounds was not in the database.

In the first chart (Figure 4B), which focuses on predictions for alkaloids, tryptophan is identified as the most prevalent precursor, accounting for 20.9 % of the predictions, followed by L-tyrosine at 17.7 %, and secologanin at 13.8 %. Other precursors contribute smaller percentages. The prevalence of amino acids, some terpenoids such as secologanin, and the aromatic acid anthranilate align well with expectations, as these are known precursors of alkaloids. However, the dominance of tryptophan, L-tyrosine, and secologanin could indicate a bias in the predictions, potentially caused by label imbalance, where overrepresented labels skew predictions in their favour (see Methods, Figure 1). Alternatively, it may reflect the actual distribution of compounds in the database, where a significantly larger proportion of compounds are produced from these most prevalent metabolites compared to others. Interestingly, despite the label imbalance, anthranilate did not appear as the 5th most common label as expected, even though it holds that position (as 5th most prevalent label for alkaloids) in our dataset (see Methods, Figure 1). Although GGPP was not included in the dataset as an alkaloid precursor, alkaloids are known to be hybrid structures that can incorporate terpenoids into their carbon skeletons. This observation suggests that the model has learned chemical patterns from terpenoids, enabling it to generalize and apply this knowledge to predict for other compound classes, even when specific precursors, like GGPP, are not explicitly present in the alkaloid dataset.

In the second chart (Figure 4C), which focuses on the precursors of terpenoids, GGPP is the most prevalent label, accounting for 49.1 % of the predictions, while FPP contributes 31.9 %. The contributions from other precursors are significantly smaller, with the third precursor (GPP) contributing 6.8 %. Another precursor, campesterol, had a lower prevalence compared to L-phenylalanine and cholesterol. The predictions involving L-phenylalanine are questionable, as terpenoids are not synthesized from this amino acid. Cholesterol, on the other hand, is a terpenoid, and despite not being explicitly included in the dataset as a terpenoid precursor, it is understandable that the model would make this prediction, especially when steroids are being analyzed. Once again, the overrepresented labels also appeared most frequently in the predictions.

The third chart (Figure 4D) focuses on the prediction of phenylpropanoid precursors. Unsurprisingly, the dominant precursor, responsible for 59.1 % of the predictions, was L-phenylalanine, followed by p-Coumaroyl-CoA at 10.2 % and Kaempferol at 7.9 %. Most of the other precursors were expected as they were assigned as phenylpropanoid precursors. FPP and GGPP were also predicted as precursors, which is questionable, as they are not known to be involved in phenylpropanoid biosynthesis. However, the fact that L-phenylalanine was predicted as a precursor for some terpenoids, and FPP and GGPP as phenylpropanoid precursors, suggests that these metabolites could be part of complex natural products, such as hybrid molecules that incorporate elements from both pathways. Once again, the most frequent predictions were associated with the overrepresented labels.

Finally, in the fourth chart (Figure 4E), focused on the production of glucosinolates, L-methionine accounts for 64.8 % of the predictions, with tryptophan contributing 20.8 % and L-tyrosine 12.6 %. These results are in accordance with the precursors defined for glucosinolates.

The predictions made for all plant compounds in LOTUS-DB highlight the effectiveness of our model in annotating a dataset of 190,273 compounds and revealing insights into the underlying biosynthetic pathways. The detailed breakdown of precursor predictions across various compound classes reveals clear patterns and interesting challenges. The model correctly identified common precursors such as tryptophan, L-phenylalanine, and GGPP, but also seems to have faced issues with overrepresented labels, which may have skewed certain predictions. For instance, in alkaloid biosynthesis, the high prevalence of tryptophan and L-tyrosine may reflect a bias in the data, while in terpenoid and phenylpropanoid classes, a few questionable predictions might indicate the complexity of natural product biosynthesis and the possibility of hybrid structures.

Overall, these results demonstrate the robustness of our predictive model, while also highlighting areas for future improvement, such as refining the handling of label imbalances and expanding the dataset to cover other biosynthetic pathways.

3.3 Ridge classifier can make and explain relevant predictions

Next, we aimed to evaluate the predictive capabilities of the Morgan FP+Ridge Classifier pipeline on several complex cases and analyze the significance assigned to specific features. This pipeline was selected for its ability to trace FP bits back to specific molecular structures, a feature that is not available with Layered FPs in the Layered FP+Ridge Classifier. We selected four compounds for this analysis: two monoterpenoid indole alkaloids recognized for their cryptic biosynthesis and two others with partially elucidated biosynthetic pathways.

Figure 5 illustrates the relationship between the model coefficients of the Morgan FP+Ridge Classifier and the substructures of the compounds. These coefficients represent the weights assigned to each feature, defining the hyperplane that separates the labels. Each coefficient corresponds to a specific feature (or molecular substructure) and indicates the extent and direction (positive or negative) of its influence on the label prediction.

Figure 5:

Predictions made by the Morgan FP+Ridge Classifier model and the importance given by the model to structural features in the compounds. On the left, bar plots show the top Ridge classifier coefficients, which correspond to bits in the FP. Red bars represent inactive negative coefficients, green bars show inactive positive bits and dark blue bars indicate active positive coefficients. The molecules subjected to the predictions are presented on the right, and the most important bits of the FP (represented in dark blue in each respective bar plot) are highlighted, where the extension of the bits (atoms and bonds within the radius defined by MorganFPs) is shaded by faded blue.

The carbon skeletons of the selected monoterpenoid indole alkaloids, quinine and ajmaline, are derived from both tryptophan and secologanin. Quinine, an alkaloid obtained from the Cinchona trees of the Rubiaceae family, has a long history as an antimalarial treatment [26]. Ajmaline, on the other hand, is an antiarrhythmic drug used in diagnosing Brugada syndrome, a rare type of arrhythmia [27].

The biosynthesis of these compounds begins with the condensation of tryptamine (derived from tryptophan) and secologanin, catalyzed by strictosidine synthase. Despite this, the substructures originating from secologanin largely disappear in the final structures, especially in the case of ajmaline. This suggests that predicting these molecules as derivatives of secologanin is challenging, given that the distinct features of secologanin may not be evident in the final chemical structure.

An experiment was conducted in which we identified the key substructures of secologanin and selected 18 molecules where these substructures largely disappear in the final structures. These molecules are still related to secologanin but with few of its substructures retained during biosynthesis, making them hard to annotate and link to the precursor. They were excluded from the training set. Despite this, secologanin was still predicted as a major precursor for these molecules, achieving a recall of 0.89, precision of 1.00, and an F1 score of 0.94, showing the model’s strong predictive ability.

Figure 5 presents the interpretability of the predictions for quinine. The most important features for predicting quinine primarily include the quinoline heterocyclic ring and a propene fragment. The high coefficient for the FP bit corresponding to the quinoline ring system may suggest an unknown step in the biosynthesis of quinine. It is believed that the quinoline structure found in cinchoninone arises from the rearrangement of cinchonaminal, which is thought to be involved in quinine’s biosynthesis [26]. Interestingly, the strong signal for the quinoline ring cannot be attributed to any significant presence of compounds with similar structures in the training set, as it did not contain such compounds; instead, these were all reserved for the test set. In contrast, the high coefficient for the propene fragment’s FP bit aligns with chemical intuition, as this fragment is the only part of the secologanin carbon skeleton retained in quinine [26].

For ajmaline, no substructure directly associated with secologanin was identified. However, rather than a propene fragment, a propane fragment exhibited a strong signal and played a crucial role in model prediction. This observation could be linked to a reduction event during the biosynthesis of ajmaline, where a propene fragment derived from the secoiridoid moiety of secologanin undergoes reduction by the 1,2-dihydrovomilenine 19,20(S)-reductase enzyme [27]. This reduction results in the propane moiety, which is highlighted in dark blue in Figure 5, further supporting the model’s predictive power in capturing biosynthetic subtleties.

The third compound featured in Figure 5, tangutorid E, is a β-carboline alkaloid found naturally in tomatoes and recently identified as a biomarker for human consumption of tomato juices [28]. β-carboline alkaloids are produced through a Pictet-Spengler reaction between indolethylamines, such as tryptophan-derived compounds, and aldehydes or α-ketoacids [29].

The model predicts that tangutorid E originates from tryptophan, with its β-carboline moiety playing a crucial role in this prediction. The presence of the indole structure and the α-amino group in tryptophan are key features influencing this outcome. This highlights the model’s focus on critical structural elements related to the biosynthesis of β-carboline alkaloids, emphasizing its ability to capture essential chemical patterns associated with these complex molecules.

The fourth compound is a phenylpropanoid produced by Vitis vinifera, Petunidin-3-(6-O-acetyl)glucoside, an anthocyanin discovered in the grapevine berry, which contributes to its blue colour. It is produced through the methylation of delphinidin glycosides [30], 31].

The model identifies delphinidin as its precursor and focuses on the structural features of anthocyanin precursors – pelargonidin, cyanidin, and delphinidin – that differ primarily by the number of hydroxyl groups attached to the benzene ring (atoms 22–27 in Figure 5). Delphinidin, for instance, contains three hydroxyl groups, while pelargonidin and cyanidin have one and two, respectively. The model assigns particular importance to the bits corresponding to these hydroxyl fragments to differentiate between the anthocyanin precursors.

Additionally, the model highlights the central pyran ring in its cationic form, which is a distinguishing feature for anthocyanin precursors, setting them apart from other phenylpropanoid precursor groups such as flavones, flavanols, and isoflavonoids. This emphasizes the model’s sensitivity to the key structural elements that define the biosynthetic pathways of anthocyanin compounds.