国家海洋药物和生物制品产业联盟

2.1. Data curation

This study utilized two approaches to collect the data. The first approach involved extracting compound information from journal articles identified through the Scopus and Web of Science databases and published within the last 5 years. The search keywords and literature-screening procedure are provided in the Supplementary Material. The second approach was to manually curate the compound list from the data repository that had been identified from seaweed https://np-mrd.org/natural_products/ [33]. The keyword used in the search was obtained by inputting the seaweed genus into the web repository, and the compound's name was retrieved. The seaweed genus used in the search was derived from the algae database (algaebase.org) [34].

2.2. Data filtering and classification

During the compound listing, this study excluded the primary metabolites and macromolecules such as polysaccharides, proteins, and lipids. This study focused on small molecules and other secondary metabolites (Table S1). From the journal references [26], [35], [36], [37], [38], [39], [40], [41], [42], [43], [44], [45], [46], [47], [48], [49], [50], [51], [52], [53], [54], [55], [56], [57], [58], [59], [60], [61], [62], [63], [64], [65], [66], [67], [68], [69], [70], [71], [72], [73], [74], [75], [76], [77], [78], [79], [80], [81], [82], [83], [84], [85], [86], [87], [88], [89], [90], [91], [92], [93], [94], [95], [96], [97], [98], [99], [100], [101], [102], [103], [104], [105], [106], [107], [108], [109], [110], [111], [112], [113], [114], [115], [116], [117], [118], [119], [120], [121], [122], [123], [124], this study grouped the collected compounds into four categories, i.e. 1) compounds that were isolated and characterized as a single pure compound with anticancer properties (Table S2, 64 compounds); 2) compounds that were isolated and characterized as a single pure compound without any anticancer assay (Table S3, 208 compounds); 3) compounds identified in the seaweed extract with anticancer properties (Table S4, 151 compounds); and 4) compounds identified in the seaweed extract without any anticancer assay (Table S5, 307 compounds). Thus, this study retrieved 730 compounds from journal references. These compounds were collected from different seaweed species across phyla.

Meanwhile, web-repository exploration was carried out to enrich the compound dataset from seaweed. From the https://np-mrd.org/natural_products/, this study retrieved 2446 compounds from 132 red seaweed species, 1686 compounds from 100 brown seaweed species, and 381 compounds from 32 green seaweed species (Table S6). Then, the compound-species duplications retrieved from the journal references and the web repository were removed to avoid redundancy in further analysis. Following duplicate checking, structure standardization, and compound–species mapping, 5203 valid occurrence records were retained for the final screening dataset. These comprised 2785 records from Rhodophyta, 1965 from Ochrophyta, and 453 from Chlorophyta (Fig. 1). At the structure level, these records represented 4083 unique standardized structures. The collected compounds were classified based on the Natural Products Classifier according to https://npclassifier.ucsd.edu/ web server into their respective pathway, superclass, and class [125]. Then, the natural products (np) class richness was calculated as Eq. (1), following the species richness measurement by Domingo-Fernández et al. [126].(1) $$R_j=\underset{c=1}{\overset{K}{\sum }}I\left(n_{\mathit{jc}}>0\right)$$where, $R_j$ is the natural products (np) class richness of species $j$, $K$ is the total number of np classes, $n_{\mathit{jc}}$ is the number of compounds belonging to class $c$ in species $j$, and $I\left(\bullet \right)$ is an indicator function that equals 1 when $n_{\mathit{jc}}>0$ ad 0 otherwise. Missing or unclassified np-class assignments were excluded.

1. Download: Download high-res image (1MB)
2. Download: Download full-size image

Fig. 1. The general schematic protocol from data curation to model development in this study.

2.3. Protein target selection

2.4. Training-set leakage and chemical-space overlap assessment

A structural overlap analysis was performed between the seaweed screening library and the final model-development dataset of each target. This analysis was conducted independently for BCL2, CASP3, CASP8, CASP9, CDK4, mTOR, and NFkB. The seaweed library and training compounds were processed using the same molecular standardization procedure. Duplicate seaweed records were collapsed by standardized structure for this analysis, resulting in 4083 unique seaweed structures. The final target-specific training datasets contained 350 compounds for BCL2, 1751 for CASP3, 252 for CASP8, 25 for CASP9, 226 for CDK4, 4285 for mTOR, and 294 for NFkB.

Exact overlap was identified by matching the standardized structure keys between the seaweed library and the corresponding training dataset. Structural similarity was then evaluated using a 2048-bit Extended-Connectivity Fingerprint, up to a maximum diameter of 4 bonds (ECFP4) generated with the Morgan algorithm at radius 2, followed by pairwise Tanimoto similarity calculation as Eq. (2) [136]. For each seaweed structure, the maximum similarity to any compound in the corresponding training set was recorded. Exact matches were separated before analysis of non-identical analogues. For descriptive analysis, non-exact structures with a Tanimoto similarity ≥0.80 were classified as highly similar analogues [137].(2) $$T\left(A,B\right)=\frac{c}{a+b-c}$$where, $T\left(A,B\right)$ is the Tanimoto similarity between compounds $A$ and $B$, $a$ is the number of fingerprint bits present in compound $A$, $b$ is the number of fingerprint bits present in compound $B$, and $c$ is the number of fingerprint bits shared by both compounds.

Scaffold-level overlap was examined using Bemis–Murcko molecular approach as Eq. (3) [138]. A seaweed structure was recorded as having scaffold overlap when its molecular framework was also present in the corresponding target-specific training dataset. Exact structure matches, maximum non-exact Tanimoto similarity, and shared-scaffold occurrence were summarized separately for each target. The compound-level results were retained to allow individual screening compounds with direct training-set overlap to be identified and flagged.(3) $$P_{\mathit{scaffold}}=\frac{N_{\mathit{shared\; scaffold}}}{N_{\mathit{seaweed}}}\times 100$$

where, $P_{\mathit{scaffold}}$ is the percentage of novel seaweed structures sharing a training-set scaffold, $N_{\mathit{shared\; scaffold}}$ is the number of novel seaweed structures with a scaffold also present in the corresponding target training set, and $N_{\mathit{seaweed}}$ is the number of seaweed structures after exact training-set overlaps were excluded.

2.5. Machine-learning model development and molecular representation

The following hyperparameter space was used. SVR used a radial basis function kernel with $C$ values of 1, 10, and 100, epsilon values of 0.05 and 0.10, and gamma = scale as Eqs. (4) and (5). RF and ETR consisted of 500 trees, with max_features evaluated at sqrt and 0.5 and min_samples_leaf at 1 and 2 as Eq. (6) [147], [149], [150]. XGB used 500 estimators, subsample = 0.8, colsample_bytree = 0.8, and evaluated maximum depths of 4 and 6 and learning rates of 0.03 and 0.05 as Eqs. (7)–(10). All ML models were executed using high-performance computing (HPC) Mahameru at the National Research and Innovation Agency of Indonesia (BRIN). The ML script was provided in the supplementary material.

2.6. Uniform true nested cross-validation

The same nested cross-validation procedure was applied to each of the seven targets by adopting the protocol of previous studies [151], [152]. The outer loop consisted of fivefold cross-validation repeated three times, producing 15 held-out outer evaluations for each target. Within each outer training partition, fourfold cross-validation with shuffling was used as the inner loop. RF, ETR, XGB, and SVR were optimized separately using only the compounds in the outer training partition. For each algorithm, the hyperparameter combination with the highest mean inner-cross-validation predictive squared correlation (Q²) was selected [141], [143], [147], [148], [150]. When two configurations showed similar Q² values, the lower inner-cross-validation root mean squared error (RMSE) was used as the secondary selection criterion. The tested hyperparameter ranges and selected settings for each algorithm are provided in Table S7 and S8 (Supplementary).

For each outer iteration, the best-performing algorithm and hyperparameter combination identified in the inner loop were refitted using the complete outer training partition. The resulting model was then evaluated only on the corresponding held-out outer test fold. Thus, compounds in the outer test fold were not used for algorithm selection or hyperparameter optimization. The predictive performance for each held-out outer fold was tested using predictive squared correlation (Q²) as Eq. (11) [153], root mean squared error (RMSE) as Eq. (12) [154], mean absolute error (MAE) as Eq. (13) [155], and Lin's concordance correlation coefficient (CCC) as Eq. (14) [156]. The mean (Eq. 15) and standard deviation (Eq. 16) for each of the 15 outside evaluations were used to summarize performance. Repeated outer forecasts for each compound were averaged to generate aggregated out-of-fold predictions and pooled out-of-fold Q², RMSE, MAE, and CCC values. Following the completion of nested performance estimation, model selection was performed on the whole curated dataset for each target using fivefold internal cross-validation.(11) $${\mathrm{Q}}^2=1-\frac{{\sum }_{i=1}^n{\left(y_i-{\mathrm{ŷ}}_i\right)}^2}{{\sum }_{i=1}^n{\left(y_i-\bar{y}\right)}^2}$$(12) $$\text{RMSE}=\sqrt{\frac{1}{n}\sum {\left(y_i-{\mathrm{ŷ}}_i\right)}^2}$$(13) $$\mathrm{MAE}=\frac{1}{n}\sum \left|y_i-{\mathrm{ŷ}}_i\right|$$(14) $$\mathrm{CCC}=\frac{2\mathit{cov}\left(y,\stackrel{̂}{y}\right)}{{\sigma }_y^2+{\sigma }_{\stackrel{̂}{y}}^2+{\left({\mu }_y-{\mu }_{\stackrel{̂}{y}}\right)}^2}$$(15) $$\bar{M}=\frac{1}{K}\underset{k=1}{\overset{K}{\sum }}{M}_k$$(16) $${\mathit{SD}}_M=\sqrt{\frac{{\sum }_{k=1}^K{\left(M_k-\bar{M}\right)}^2}{K-1}}$$where, $y_i$ is observed pIC₅₀ of compound i, ${ŷ}_i$ is predicted pIC₅₀ of compound i, $\bar{y}$ is mean pIC₅₀ of observed pIC₅₀ values, $n$ is number of compounds in the evaluated set, ${\mu }_y$ is mean of the observed values, ${\mu }_{\stackrel{̂}{y}}$ is mean of the predicted values, ${\sigma }_y^2$ is variance of the observed values, ${\mu }_{\stackrel{̂}{y}}$ is variance of the predicted values, $\mathit{cov}\left(y,\stackrel{̂}{y}\right)$ is covariance between observed and predicted values, $\bar{M}$ is the mean across the 15 outer evaluations, $M_k$ is the metric obtained in outer evaluation k, ${\mathit{SD}}_M$ is standard deviation, and K = 15.

2.7. Applicability-domain analysis

An applicability-domain (AD) assessment was performed for each target-specific machine-learning model after model development and validation and before the interpretation and multi-target prioritization of the seaweed metabolites following previous studies [157], [158]. Molecular frameworks were standardized with RDKit and represented as ECFP4 fingerprints generated using the Morgan algorithm with radius 2 and 2048 bits without chirality. Structural similarity was calculated using the Tanimoto coefficient as Eq. (2).

For each target, the structural AD threshold was determined independently from the empirical distribution of these leave-one-out (LOO) nearest-neighbour similarities as Eq. (17) [159].(17) $$s_{i,t}^{\mathit{LOO}}=\underset{\begin{array}{c}j=1,\dots ,n_t\\ j\ne i\end{array}}{\max }T\left(F_i,F_j\right)$$where, $s_{i,t}^{\mathit{LOO}}$ is the highest similarity of the training compound $i$ to another compound in the same target-specific training set, $n_t$ is the training set for target $t$, $F_i$ and $F_j$ are the ECFP4 fingerprints of compounds $i$ and $j$, respectively, and $T\left(F_i,F_j\right)$ is their Tanimoto similarity.

The fifth percentile (Q05) of this distribution was used as the primary criterion as Eq. (18).(18) $${\theta }_t=Q_{05}\left(s_{1,t}^{\mathit{LOO}},s_{2,t}^{\mathit{LOO}},\dots ,s_{n_t,t}^{\mathit{LOO}}\right)$$where, ${\theta }_t$ is the structural AD threshold for target $t$, $Q_{05}$ denotes the fifth percentile, and $n$ is the number of compounds in the target-specific training set. This approach allowed the AD boundary to follow the structural diversity of each training dataset instead of applying a single Tanimoto cutoff to all seven models. For each seaweed compound $x$, its maximum similarity to the training set of target $t$ was calculated as Eq. (19) [159], [160].(19) $$s_t\left(x\right)=\underset{j=1,\dots ,n_t}{\max }T\left(F_x,F_j\right)$$where $s_t\left(x\right)$ is the maximum similarity of seaweed compound $x$ to the target-specific training set, $F_x$ is its ECFP4 fingerprint, and $F_j$ represents the fingerprint of training compound $j$.

The compounds were then assigned to the structural AD based on the highest resemblance to the associated target-specific training set as Eq. (20).(20) $${\mathit{AD}}_{\mathit{struct},t}\left(x\right)=\left\{\begin{array}{c}1,s_t\left(x\right)\ge {\theta }_t\\ 0,s_t\left(x\right)<{\theta }_t\end{array}\right.$$where, ${\mathit{AD}}_{\mathit{struct},t}\left(x\right)=1$ means compound $x$ is inside the structural AD for target $t$, ${\mathit{AD}}_{\mathit{struct},t}\left(x\right)=0$ means compound $x$ is outside the structural AD, $s_t\left(x\right)$ is maximum similarity of compound $x$ to the target-specific training set, ${\theta }_t$ is target-specific AD threshold.

By repeating the classification using the first, fifth, and tenth percentiles (Q01, Q05, and Q10) of the training leave-one-out similarity distribution, the robustness of this definition was investigated. The similarity categories of 0.60 ≤ Tanimoto < 0.80 and Tanimoto ≥ 0.80 were not employed as universal AD criteria; rather, they were primarily utilized for descriptive chemical-space summaries.

A response-domain criterion was also evaluated for each target. A prediction was considered within the response domain when its predicted pIC₅₀ value remained within the experimental pIC₅₀ range represented by the corresponding training dataset as Eq. (21) [161].(21) $$y_{t,\mathit{min}}\le {\stackrel{̂}{y}}_t\left(x\right)\le y_{t,\mathit{max}}$$where ${\stackrel{̂}{y}}_t\left(x\right)$ is the predicted pIC₅₀ of compound $x$ for target $t$, and $y_{t,\mathit{min}}$ and $y_{t,\mathit{max}}$ are the minimum and maximum experimental pIC₅₀ values in the target-specific training set, respectively. A novel compound was considered within the combined AD when it satisfied both the structural and response-domain criteria. Exact standardized-structure matches between the seaweed library and the training data were identified separately and were flagged as training-set overlaps rather than novel external predictions.

2.8. Multi-objective prioritization and uncertainty analysis

Following the research of Fromer et al. and Angelo et al. [162], [163], a multi-objective Pareto approach was used because the seven target-specific models differed in training-set size, chemical space, predictive performance, error, and response distribution. All experimental IC₅₀ values had been standardized to molar concentration and transformed to pIC₅₀ before model development. Experimental IC₅₀ values were converted to pIC₅₀ using $p{\mathit{IC}}_{50}=-{\log }_{10}\left({\mathit{IC}}_{50}\left[\text{in}\mathrm{M}\right]\right)$; equivalently, $p{\mathit{IC}}_{50}=9-{\log }_{10}\left({\mathit{IC}}_{50}\left[\text{in}\mathrm{nM}\right]\right)$. All seven models therefore generated predictions on the same logarithmic molar pIC₅₀ scale. Pairwise pIC₅₀ differences were retained as the primary screening contrasts because they preserve both the direction and magnitude of the predicted potency difference. However, the target-specific models were trained independently using datasets with different sizes, chemical spaces, prediction errors, and response distributions. Thus, ΔpIC₅₀ was interpreted as a model-based screening contrast rather than an experimentally calibrated cross-target selectivity value. Target-specific prediction error was considered through uncertainty propagation. Sensitivity to differences in prediction distribution was examined separately using percentile normalization.

BCL2, CDK4, mTOR, and NFkB were the four oncogenic targets that were considered as inhibition targets, with larger anticipated pIC₅₀ being advantageous. CASP3, CASP8, and CASP9 were addressed differently since significant suppression of these apoptosis-related proteins was undesirable in the proposed priority scheme. In this study, “caspase sparing” refers only to lower predicted inhibitory potency toward CASP3, CASP8, and CASP9. It does not indicate predicted or experimentally demonstrated caspase activation or apoptosis induction.

2.9. Chemical-class enrichment analysis

Chemical-class enrichment was investigated to see if certain chemical classes were more common among prioritized compounds than in the entire dataset. The analysis was performed separately for chemical class and chemical superclass. Three groups of prioritized compounds were evaluated, i.e., the top 5%, the top 10%, and the final Robust Priority Front 1 [171]. The complete set of 5203 compounds was used as the background. The enrichment factor ( ${\mathit{EF}}_{\mathit{class}}$) was calculated as Eq. (31).(31) $${\mathit{EF}}_{\mathit{class}}=\frac{a/\left(a+b\right)}{\left(a+c\right)/\left(a+b+c+d\right)}$$

The numerator represents the frequency of the tested chemical class in the prioritized group. The denominator represents its frequency in the complete dataset. Where, $a$ is number of compounds from the tested chemical class in the prioritized subset, $b$ is number of compounds from all other chemical classes in the prioritized subset, $c$ is the number of compounds from the tested chemical class outside the prioritized subset, and $d$ is the number of compounds from all other chemical classes outside the prioritized. An ( $\mathit{EF}$ > 1) indicates enrichment. An ( $\mathit{EF}$ = 1) indicates the same frequency as the background. An ( $\mathit{EF}$ < 1) indicates depletion.

The odds ratio ( $\mathit{OR}$) was calculated as Eq. (32). An $\mathit{OR}$ > 1 indicates that the tested chemical class occurred more frequently among prioritized compounds relative to the background, whereas $\mathit{OR}$ < 1 indicates under-representation [172].(32) $$\mathit{OR}=\frac{\mathit{ad}}{\mathit{bc}}$$

Statistical significance was tested using a one-sided Fisher's exact test as Eq. (33). For the enrichment test, the probabilities of the observed table and all tables showing stronger enrichment were summed to obtain the one-sided ( $P$) value. Fisher's exact test was selected because some chemical classes contained small numbers of compounds [173].(33) $$P_{\mathit{enrich}}=\underset{k=a}{\overset{\mathit{min}\left(a+b,a+c\right)}{\sum }}\frac{\left(\begin{array}{c}a+b\\ a\end{array}\right)\left(\begin{array}{c}c+d\\ c\end{array}\right)}{\left(\begin{array}{c}a+b+c+d\\ a+c\end{array}\right)}$$

Because several chemical classes were tested, multiple comparisons were controlled using the Benjamini-Hochberg false discovery rate method, Eq. (34). The ( $P$) values were ordered from the smallest to the largest [172], [173], [174].(34) $$p_{\left(i\right)}\le \frac{i}{m}\alpha$$where ( $i$) is the rank of the ( $p$) value, ( $m$) is the total number of tests, and $\alpha$ = 0.05.

To support reproducibility, fixed random seeds were used for stochastic procedures, including data splitting, cross-validation, model fitting, hyperparameter optimization, and random sampling, where applicable. The data sets used for model development, train/test split information, selected model parameters and hyperparameters, and scripts used to generate the main results are provided in the Supplementary material.

2.10. Molecular docking and retrospective active–decoy benchmarking

The docking analysis was performed after completion of the multi-objective prioritization. The final Robust Priority Front 1 contained 16 occurrence records corresponding to 12 unique standardized structures. Because molecular docking was performed at the structure level, repeated occurrence records sharing the same standardized structure were represented once in the docking set. The 12 unique structures were each docked against CASP3 (PDB ID: 2XYG) [175], CASP8 (PDB ID: 1QTN) [176], CASP9 (PDB ID: 1JXQ) [177], BCL2 (PDB ID: 6O0K) [178], CDK4 (PDB ID: 7SJ3) [179], mTOR (PDB ID: 4JT6) [180], [181], and NFkB (PDB ID: 1IKN) [182]. Thus, the final candidate docking analysis contained 84 structure–target combinations (12 structures × 7 targets). All docking and validation were executed using high-performance computing (HPC) Mahameru at the National Research and Innovation Agency of Indonesia (BRIN). The protein crystallization model was obtained from the Research Collaboratory for Structural Bioinformatics (RCSB) database (https://www.rcsb.org/) [183]. AutoDock Vina was used for molecular docking, and more negative Vina scores (kcal/mol) were interpreted as more favorable docking scores within the same protein target [184]. Docking scores were interpreted separately for each protein target. Because Vina scores from different receptor systems are not assumed to be directly comparable on an absolute scale, cross-target mean scores were not used for compound ranking or as estimates of overall binding affinity.

Three-dimensional ligand structures were generated from standardized Simplified Molecular Input Line Entry System (SMILES). Hydrogen atoms were added and initial conformations were generated using the Experimental Torsion Knowledge Distance Geometry (ETKDG) approach implemented in RDKit [185]. Geometry optimization was performed using Merck Molecular Force Field (MMFF) when the required parameters were available and Universal Force Field (UFF) otherwise. The structures were subsequently converted to the Protein Data Bank + Partial Charge + Atom Type (PDBQT) format for AutoDock Vina calculations. The same receptor preparation, docking box, and Vina parameters were used for prioritized metabolites and reference compounds to allow direct target-wise comparison (Table S10, supplementary).

Native-ligand redocking was retained as a test of pose reproduction (RMSD < 2) [186]. However, redocking Root-Mean-Square Deviation (RMSD) was not considered evidence of virtual-screening enrichment. Therefore, a separate active–decoy benchmark was performed for each target [187]. Twenty known active compounds and 100 decoys were used for each protein. Database of Useful Decoys - Enhanced (DUD-E) reference decoys were used where suitable target-specific sets were available [188], while LIDeB's Useful Decoys (LUDe)-generated decoys were used for CDK4 and NFkB [189]. These approaches provide property-matched decoys for retrospective evaluation of virtual-screening methods. All reference molecules were processed and docked using the same workflow as the seaweed compounds.

Docking discrimination was evaluated using the Receiver Operating Characteristic–Area Under the Curve (ROC-AUC) and enrichment factors (EF) at the top 1%, 5%, and 10% of each target-specific ranked list. ROC-AUC describes the ability of a ranking method to separate known actives from decoys over the complete ranked set, while the enrichment factor measures the concentration of actives near the top of the ranked list [190]. The true-positive and false-positive rates were defined as Eqs. (35)–(37).(35) $$\mathrm{ROC}-\mathrm{AUC}={\int }_0^1\mathrm{TPR}\left(\mathrm{FPR}\right)d\left(\mathrm{FPR}\right)$$(36) $$\mathrm{TPR}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}$$(37) $$\mathrm{FPR}=\frac{\mathrm{FP}}{\mathrm{FP}+\mathrm{TN}}$$where TP and FN are active compounds classified above and below a given ranking threshold, respectively, and FP and TN are the corresponding decoys.

For the top fraction (x) of the ranked library, the enrichment factor ( ${\mathit{EF}}_{x\%}$) was calculated as Eq. (38).(38) $${\mathit{EF}}_{x\%}=\frac{a_{x}N}{n_{x}A}$$

where a_x is the number of known actives found within the top x fraction, n_x is the number of compounds in that fraction, A is the total number of actives, and N is the total number of benchmark compounds. The retrospective active–decoy analysis was used to characterize docking performance separately for each target. Because discrimination differed among targets, docking scores were interpreted according to the corresponding target-specific benchmarking performance. Docking was used as complementary structural support and not as validation of the ML-predicted activities.

3.1. Compounds' diversity

The distribution of natural products was uneven across the investigated seaweed species and phyla (Fig. 2a). Based on the compounds collected from literature and databases, red seaweeds (Rhodophyta) contributed the largest number of compounds in the curated dataset. In comparison, brown seaweeds (Ochrophyta) and green seaweeds (Chlorophyta) contributed fewer compounds. The large number of compounds reported from red seaweeds reflects both their rich secondary metabolism and the extensive investigation of this group in marine natural-product research.

1. Download: Download high-res image (5MB)
2. Download: Download full-size image

Fig. 2. The compound diversity across seaweed species and phylum: (a) The diversity of compound superclass and class based on the seaweed species retrieved in this study, (b) the top 30 species with the richest compounds' class, and (c) the top 30 compound class distribution across three seaweed phyla.

The NP class richness analysis showed considerable variation among species (Fig. 2b). Class richness was used as an indicator of chemical diversity by counting the number of distinct natural product classes reported for each species. Laurencia dendroidea exhibited the highest richness score, followed by Dictyota dichotoma and Laurencia obtusa. Species belonging to the genus Laurencia dominated the top-ranked species, with several representatives appearing among the 30 most chemically diverse seaweeds. Other red seaweed genera, including Gracilaria and Plocamium, were also frequently represented. Among brown seaweeds, Dictyota, Undaria, Saccharina, Fucus, and Sargassum appeared within the top 30 species. This trend indicated substantial chemical diversity in selected members of Ochrophyta. Green seaweeds were less represented, although species from the genera Ulva and Halimeda were also present among the most compound-rich taxa.

The distribution of natural product classes revealed clear chemotaxonomic differences between seaweed phyla (Fig. 2c). Red seaweeds were strongly enriched in sesquiterpene-derived metabolites, particularly chamigrane sesquiterpenes, laurene sesquiterpenes, and bisabolene sesquiterpenes. These classes were abundant in Rhodophyta, although they occurred in limited numbers in brown seaweeds. In contrast, brown seaweeds were characterized by prenyl quinone meroterpenoids, dolabellane diterpenoids, and tetraketide meroterpenoids, which were much less common in red seaweeds. The heatmap also showed that several classes were almost exclusively associated with one phylum, suggesting distinct biosynthetic capabilities among macroalgal lineages. Green seaweeds generally contained fewer compounds and fewer dominant classes, although some metabolite groups were shared with both Rhodophyta and Ochrophyta.

3.2. Target-specific training datasets and activity distributions

The final training datasets differed considerably in size among the seven targets. The largest dataset was obtained for mTOR, with 4285 compounds, followed by CASP3 with 1751 compounds. The remaining datasets contained 350 compounds for BCL2, 294 for NFkB, 252 for CASP8, 226 for CDK4, and only 25 for CASP9. Thus, the amount of experimental activity data available for model development was strongly target-dependent. The experimental pIC₅₀ distributions also differed among the targets (Fig. 3). The pIC₅₀ ranges were 2.37–9.09 for BCL2, 3.05–10.00 for CASP3, 4.00–9.92 for CASP8, 4.03–8.70 for CASP9, 3.63–9.15 for CDK4, 2.00–10.15 for mTOR, and 4.01–12.00 for NFkB. Median experimental pIC₅₀ values also varied among the datasets, with the lowest median observed for CASP3 (5.29) and the highest for CDK4 (7.49).

1. Download: Download high-res image (592KB)
2. Download: Download full-size image

Fig. 3. Final sizes and experimental pIC₅₀ distributions of the seven target-specific ChEMBL datasets after data audit and curation. The dashed vertical lines indicate the median experimental pIC₅₀ for each target.

3.3. Training-set overlap and chemical-space comparison

3.4. Model comparison and selection under nested cross-validation

3.5. Applicability domain of seaweed metabolites

The structural applicability domain differed strongly among the seven target-specific models (Fig. 4, Table 5). The analysis was performed for 4083 unique standardized seaweed structures, with exact target-specific training-set matches excluded from the novel-compound coverage calculation. The primary structural AD was defined using the fifth percentile (Q05) of the leave-one-out nearest-neighbour Tanimoto similarity distribution of each training set. The Q05 threshold ranged from 0.207 for CASP9 to 0.658 for mTOR. NFkB showed the largest coverage of the seaweed chemical space, with 526 structures (12.88%) inside the Q05 structural AD. Coverage was much lower for CASP9 (102 structures (2.50%)), CASP3 (84 (2.06%)), and BCL2 (65 (1.59%)). No novel seaweed structure met the Q05 threshold for CASP8, CDK4, or mTOR. Coverage changed when Q01 or Q10 was used, confirming that the AD classification was sensitive to the selected boundary. No novel seaweed structure was inside the Q05 structural AD of all seven targets. Thus, most large-scale predictions represented extrapolation relative to the corresponding target-specific training chemical space and were interpreted with greater caution.

1. Download: Download high-res image (375KB)
2. Download: Download full-size image

Fig. 4. Applicability-domain assessment of the seven target-specific machine-learning models: (a) Percentage of seaweed structures inside the primary structural applicability domain defined at the Q05 threshold, (b) Target-specific Tanimoto similarity thresholds derived from the Q01, Q05, and Q10 percentiles of leave-one-out nearest-neighbour similarities in each training set, (c) Sensitivity of structural applicability-domain coverage to the selected threshold.

Table 5. Target-specific structural applicability-domain coverage of the seaweed screening library.

Target	Q05 Tanimoto threshold	Novel structures evaluated, n	Inside structural AD, n (%)	Outside structural AD, n (%)	Exact training matches, n
BCL2	0.2765	4083	65 (1.59%)	4018 (98.41%)	0
CASP3	0.3226	4083	84 (2.06%)	3999 (97.94%)	0
CASP8	0.5801	4083	0 (0.00%)	4083 (100.00%)	0
CASP9	0.2075	4083	102 (2.50%)	3981 (97.50%)	0
CDK4	0.3764	4083	0 (0.00%)	4083 (100.00%)	0
mTOR	0.6575	4082	0 (0.00%)	4082 (100.00%)	1
NFkB	0.2753	4083	526 (12.88%)	3557 (87.12%)	0

3.6. Robust multi-target prioritization under prediction uncertainty

The target-distribution sensitivity analysis showed a Spearman correlation of ρ = 0.601 between the primary pIC₅₀-based and percentile-normalized rankings across the 5203 records. Twelve of the 16 primary Robust Priority Front-1 records (75.0%) were also included in the percentile-normalized Front 1. This result indicates that most final candidates remained supported after removal of target-specific prediction-scale differences, although the overall ranking was moderately sensitive to normalization. CDK4 was the strongest predicted oncogenic target for 5169 compounds (99.35%). CASP3 had the lowest predicted inhibitory potency among the three caspases for 5182 compounds (99.60%) and was described as the most-spared caspase in the computational screening. Here, sparing refers to weaker predicted inhibition and not to caspase activation. The largest median oncogenic–caspase contrasts involved CDK4-CASP3 and CDK4-CASP8, followed by mTOR-CASP3.

Propagation of target-level model error changed the stability of several contrasts. Under the RMSE-based Gaussian perturbation model, CDK4-CASP3 remained the most stable oncogenic–caspase comparison, while CDK4-CASP8 also provided substantial robust-pair support. In contrast, none of the four oncogenic-CASP9 pairs reached the predefined pair-stability threshold of ≥0.95. Thus, CASP9 did not contribute a robust positive pair to the criterion for any compound. Its influence on the final ranking was mainly through the lower pair-stability values and potentially the minimum pair-stability criterion. CASP9 acted as an uncertainty-related constraint rather than as positive evidence supporting final candidate selection.

The final ranking was based on the joint trade-off among Monte Carlo oncogenic Pareto Front-1 stability, the number of robust positive oncogenic–caspase pairs, and minimum pair stability (Fig. 5). Sixteen occurrence records, corresponding to 12 unique standardized structures, formed the final Robust Priority Front 1, including 5-OHKF (SW5194), hydroxypheophytin a (SW0481), pheophytin a (SW0482), aplysiatoxin (SW2423), 6-hydroxystigmasta-4,22-dien-3-one (SW0494), pheophorbide A (SW0476), pheophorbide B (SW0325), and clerosterol 3-glucoside (SW4827) (Fig. 6, Fig. 7). Since some compounds have International Union of Pure and Applied Chemistry (IUPAC) names, and those are long names, we put the complete list of Robust Priority Front 1 into Table S11 (Supplementary).

1. Download: Download high-res image (557KB)
2. Download: Download full-size image

Fig. 5. Monte Carlo Pareto Front-1 stability frequency plotted against the minimum oncogenic–caspase pair stability for 5203 compounds.

1. Download: Download high-res image (464KB)
2. Download: Download full-size image

Fig. 6. Predicted pIC₅₀ profiles of the 16 Robust Priority Front-1 occurrence records across the seven targets.

1. Download: Download high-res image (540KB)
2. Download: Download full-size image

Fig. 7. Pairwise Monte Carlo stability of the 16 Robust Priority Front-1 occurrence records. Stability values represent the fraction of RMSE-based Gaussian perturbations in which each of the 12 oncogenic–caspase contrasts remained positive. These values are model-based stability measures and not experimentally calibrated probabilities.

3.7. Chemical-class enrichment across compound-priority thresholds

Chemical-class enrichment was evaluated at the Top 5%, Top 10%, and final Robust Priority Front-1 thresholds (Fig. 8). Here, enrichment refers to statistical overrepresentation of a chemical class among computationally prioritized compounds relative to its background frequency in the complete dataset. A high EF does not indicate greater biological potency or stronger in vitro or in vivo anticancer activity of compounds belonging to that class. At the Top 5% threshold, several chemical classes remained significantly enriched after FDR correction. Related enrichment patterns were also observed in the Top 10% set, although the magnitude and statistical support varied among classes and priority thresholds.

1. Download: Download high-res image (406KB)
2. Download: Download full-size image

Fig. 8. Chemical-class enrichment across the compound-priority thresholds: (a) top 5%, (b) top 10%, and (c) robust priority Front-1. Red bars indicate classes significant after Benjamini–Hochberg correction (FDR < 0.05), whereas blue bars indicate non-significant enrichment. The dashed line indicates EF = 1, *FDR < 0.05. only classes with ≥10 compounds in the complete dataset are shown. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

The enrichment pattern changed in the final 16-record Robust Priority Front 1. Four occurrence records were classified as stigmastane steroids, while two records were classified as prenyl quinone meroterpenoids. Aminosugars, cyclic peptides/depsipeptides/lipopeptides, aplysiatoxins, and carboline alkaloids were each represented by one Front-1 record. Several of these classes showed high enrichment factors because their background frequencies in the complete dataset were low. However, high EF values based on only one or a few Front-1 records should be interpreted cautiously, particularly when statistical significance was not retained after multiple-testing correction.

3.8. Molecular docking investigation

The final Robust Priority Front 1 contained 16 occurrence records corresponding to 12 unique standardized structures. Repeated occurrence records sharing the same standardized structure were represented once in the docking analysis. These 12 structures were docked against all seven targets, giving 84 structure–target combinations (Fig. 9). Vina scores ranged from −10.637 to −4.223 kcal/mol. The most negative scores were observed for SW0500 against mTOR (−10.637 kcal/mol) and for SW0494 against CDK4 (−10.507 kcal/mol). SW0500 also showed a low score against NFkB (−9.247 kcal/mol), while SW0494 gave the lowest score for BCL2 (−8.788 kcal/mol) and CASP8 (−8.004 kcal/mol). Within the other caspases, the lowest scores were obtained for SW0500 against CASP3 (−7.767 kcal/mol) and SW4832 against CASP9 (−7.818 kcal/mol).

1. Download: Download high-res image (769KB)
2. Download: Download full-size image

Fig. 9. Molecular docking-score profiles of the final Robust Priority Front-1 structures across seven targets.

Agreement between the ML predictions and docking rankings was target-dependent and generally limited (Fig. 10). Spearman correlations ranged from −0.816 for CASP3 to 0.116 for mTOR. Negative correlations were observed for CASP3 (ρ = −0.816), CASP8 (ρ = −0.599), BCL2 (ρ = −0.508), CASP9 (ρ = −0.403), and CDK4 (ρ = −0.144), while weak positive correlations were obtained for mTOR (ρ = 0.116) and NFkB (ρ = 0.053). The inverse CASP3 correlation was statistically significant and remained significant after Benjamini–Hochberg correction (p = 0.0012; FDR = 0.0084). CASP8 showed a nominal association (p = 0.0396), yet it did not remain significant after multiple-testing correction (FDR = 0.139). The remaining target-wise correlations were not statistically significant after correction. The significant inverse CASP3 association indicates disagreement between the ML and docking rankings rather than independent validation of the ML predictions. Docking was retained as a complementary structural assessment and was not used to rerank or confirm the ML predictions.

1. Download: Download high-res image (344KB)
2. Download: Download full-size image

Fig. 10. Target-wise ML–docking rank concordance for the 12 final unique structures. Points show Spearman ρ between predicted pIC₅₀ and −Vina score, with 95% confidence intervals. FDR values were obtained using Benjamini–Hochberg correction.

3.9. Taxonomic distribution and species-normalized prioritization

The 16 occurrence records in the final Robust Priority Front 1 were reported from 11 seaweed species. Dictyota dichotoma contributed four Front-1 records, while D. fasciola and Codium decorticatum each contributed two records; the remaining species contributed one record each. Normalization by the number of screened records changed the apparent species priority. Among species represented by at least 10 records, C. decorticatum showed the highest Front-1 hit rate, with 2 of 13 records (15.38%). This was followed by Bryopsis pennata (1/16 (6.25%)), Halimeda gracilis (1/22 (4.55%)), Gracilaria coronopifolia (1/24 (4.17%)), D. fasciola (2/57 (3.51%)), and Caulerpa racemosa (1/50 (2.00%)). In contrast, the more extensively represented D. dichotoma, Sargassum siliquastrum, and Laurencia dendroidea showed lower normalized hit rates of 1.33%, 0.76%, and 0.27%, respectively. At genus level, Gracilaria contained 2 Front-1 records among 147 screened records (1.36%), compared with 6/657 (0.91%) for Dictyota and 1/1684 (0.059%) for Laurencia (Fig. 11). These normalized hit rates are descriptive proportions and should not be interpreted as statistical evidence of species-level chemical enrichment. In particular, high rates based on small denominators can be unstable even when a minimum threshold of 10 screened records is applied.

1. Download: Download high-res image (367KB)
2. Download: Download full-size image

Fig. 11. Species-normalized representation of the final Robust Priority Front 1. Labels indicate the number of Front-1 records relative to the total number of screened records for each species. The dashed line indicates the overall Front-1 rate across the complete screening dataset.