
==== Front
Brief Bioinform
Brief Bioinform
bib
Briefings in Bioinformatics
1467-5463
1477-4054
Oxford University Press

10.1093/bib/bbae046
bbae046
Problem Solving Protocol
AcademicSubjects/SCI01060
Multiple and Optimal Screening Subset: a method selecting global characteristic congeners for robust foodomics analysis
Xu Rui Human Nutrition Program, Department of Human Sciences, The Ohio State University, Columbus, Ohio, USA  43210
Comprehensive Cancer Center, The Ohio State University, Columbus, Ohio, USA  43210

Zhang Huan Human Nutrition Program, Department of Human Sciences, The Ohio State University, Columbus, Ohio, USA  43210
Comprehensive Cancer Center, The Ohio State University, Columbus, Ohio, USA  43210

Crowder Michael W Department of Chemistry and Biochemistry, Miami University, Oxford, Ohio, USA  45056

https://orcid.org/0000-0002-4548-8949
Zhu Jiangjiang Human Nutrition Program, Department of Human Sciences, The Ohio State University, Columbus, Ohio, USA  43210
Comprehensive Cancer Center, The Ohio State University, Columbus, Ohio, USA  43210

Corresponding authors: Michael W. Crowder, Department of Chemistry and Biochemistry, Miami University, Oxford, OH 45056, USA. Tel.: +513-529-3754; Fax: 614-292-7969; E-mail: crowdemw@miamioh.edu; Jiangjiang Zhu, Comprehensive Cancer Center and Department of Human Sciences, The Ohio State University, 400 W 12th Ave, Columbus, OH 43210, USA. Tel.: +614-685-2226; Fax: 614-292-7969; E-mail: zhu.2484@osu.edu
3 2024
21 2 2024
21 2 2024
25 2 bbae04628 8 2023
4 1 2024
26 1 2024
© The Author(s) 2024. Published by Oxford University Press.
2024
https://creativecommons.org/licenses/by-nc/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (https://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com

Abstract

Metabolomics and foodomics shed light on the molecular processes within living organisms and the complex food composition by leveraging sophisticated analytical techniques to systematically analyze the vast array of molecular features. The traditional feature-picking method often results in arbitrary selections of the model, feature ranking, and cut-off, which may lead to suboptimal results. Thus, a Multiple and Optimal Screening Subset (MOSS) approach was developed in this study to achieve a balance between a minimal number of predictors and high predictive accuracy during statistical model setup. The MOSS approach compares five commonly used models in the context of food matrix analysis, specifically bourbons. These models include Student’s t-test, receiver operating characteristic curve, partial least squares-discriminant analysis (PLS-DA), random forests, and support vector machines. The approach employs cross-validation to identify promising subset feature candidates that contribute to food characteristic classification. It then determines the optimal subset size by comparing it to the corresponding top-ranked features. Finally, it selects the optimal feature subset by traversing all possible feature candidate combinations. By utilizing MOSS approach to analyze 1406 mass spectral features from a collection of 122 bourbon samples, we were able to generate a subset of features for bourbon age prediction with 88% accuracy. Additionally, MOSS increased the area under the curve performance of sweetness prediction to 0.898 with only four predictors compared with the top-ranked four features at 0.681 based on the PLS-DA model. Overall, we demonstrated that MOSS provides an efficient and effective approach for selecting optimal features compared with other frequently utilized methods.

Graphical Abstract

Graphical Abstract

molecular feature selection
optimal screening
mass spectrometry
foodomics
metabolomics
National Institute of General Medical Sciences of the National Institutes of Health R35GM133510
==== Body
pmcINTRODUCTION

The field of foodomics, metabolomics, and other omics analyses often faces certain challenges that can introduce elements of arbitrariness throughout its entire process. For instance, in the preliminary stage of data preprocessing, a commonly used software such as MetaboAnalyst suggests eliminating the lowest 40% of features based on their average abundance. However, there are also studies that used mean-difference plots to offer a more feasible and arbitrary truncation point for users by aligning biological samples with blanks. This allows for a better equilibrium between eliminating low-quality features and retaining high-quality features [1]. Meanwhile, during the compound identification process, it was recognized as a common mistake to use a top or an arbitrary match for subsequent pathway analysis, because the top match is usually only one of many possible candidates for the true features [2]. To solve the arbitrary top match issue in compound identification, the mummichog algorithm proposed a strategy by accommodating the ambiguity of metabolite identification in the statistical framework of pathway analysis. The algorithm searches the path using a combination of all candidate paths, containing most of the incorrectly predicted candidates and combinations, and computes the enrichment P-value for each of them, outputting the few results that will fall on the true metabolite and the correct pathway [3].

In our previous foodomics/metabolomics data analysis, between data preprocessing and feature identification, we utilized statistical models for feature selection. In this step, we typically reduce the feature dimensions by eliminating features that contribute minimally to the classification models, aiming to achieve a balance between high predictive accuracy and a minimal set of predictors. From a practical application standpoint, minimizing the complexity of input variables while enhancing model predictive capability is crucial for deploying models in real-world settings, thereby aiding disease diagnosis, efficacy assessment, product classification, and other similar applications. To enable feasible feature selection, a model-based approach utilizing a variable importance score is often employed, which ranks features based on their contribution to a classification problem. In a critical review of metabolomics data analysis, univariate analysis, unsupervised analysis, supervised analysis, and receiver operating characteristic (ROC) curve analysis were identified as essential criteria for ensuring comprehensive data analysis [4]. Many studies have suggested that supervised analysis using partial least squares–discriminant analysis (PLS-DA)/orthogonal partial least squares–discriminant analysis (OPLS-DA) is the most frequently employed method, followed by random forests (RFs) and support vector machines (SVMs) [5, 6]. Although reduced dimensions were found to account for the resulting variables, variable-picking models typically provide only a measure of the importance of the variables and do variable selection unnaturally. The list generated through this process offers the analyst a ranking of the most crucial metabolites/molecules or mass spectral features, but there is no good way to select a threshold to determine which features are considered significant. Conventionally, a commonly used threshold, such as P-value < 0.05 or 0.01, was abused in this study [7]. Nonetheless, it is worth noting that these thresholds may not universally apply to all biomarker selection scenarios, which could result in arbitrary feature selection that does not align with the specific context [8]. In essence, a definitive method for identifying the metabolites of utmost significance in predicting outcomes remains elusive, potentially resulting in suboptimal results.

To simultaneously achieve non-arbitrary feature selection in contextualized models, high prediction accuracy, and reduced computational complexity, we developed a three-step method, named Multiple and Optimal Screening Subset (MOSS). First, an initial candidate list of data features with ambiguous significance is generated. Then, an optimal subset size is determined by comparing the multiple-candidate combination to the corresponding top-ranked features. Finally, the optimal feature subset for classification is determined by traversing all possible feature combinations with the optimal size from feature candidates. We utilized an example foodomics analyses of bourbon whiskey and applied our newly developed strategy for the search of molecular biomarkers that can differentiate these commercial bourbons based on their age, proof, and taste. Our study is innovative because, through the three-step process, we not only avoid the need to determine an arbitrary threshold but also prevent the large computational costs. We demonstrated that the new method we developed can achieve a balance between a minimal set of predictors and high predictive accuracy.

METHODS

Materials

All the solvents including LC/MS-grade methanol, acetonitrile, water, ammonium acetate, formic acid, and acetic acid were purchased from Fisher Scientific (Pittsburgh, PA, USA), and 121 commercial bourbon products were collected from local retailers and stored in clear glass vials at 4°C. Detailed information regarding the bourbon samples is available in Table S1 (see Supplementary Data available online at http://bib.oxfordjournals.org/).

Sample preparation

Each bourbon sample was diluted 1:10 in a mixture of methanol and water (1,1) by adding 100 μl, in accordance with a published method [9]. To prepare the quality control (QC) samples, aliquots of each bourbon sample (0.5 ml) were combined, treated as described above, and injected after every 10 samples.

Sample analysis

A Thermo Vanquish UPLC system (ThermoFisher, Waltham, MA, USA) coupled with a Hybrid Quadrupole Orbitrap Q Exactive™ mass spectrometer (ThermoFisher) was used for small molecule analysis of study samples. Bourbon molecules were separated using an Xbridge BEH Amide (2.5 μm, 2.1 × 150 mm, Waters, Milford, MA, USA) column. The mobile phases consisted of solvents A and B, where A contained 5 mM ammonium acetate, 0.1% acetic acid in 90% H2O/10% acetonitrile (ACN), and B contained 5 mM ammonium acetate, 0.1% acetic acid in 10% H2O/90% ACN. A gradient of mobile phase injection was used for 20 min at a flow rate of 0.3 ml/min, starting with 30% of mobile phase A and increasing to 70% at 5 min. The proportion of mobile phase A was maintained at 70% until 9 min and then started to decrease to 30% until 11 min. The percentage of solvent A was then kept at 70% until 20 min, after which it gradually returned to 30% to prepare for the next injection. To observe and prove the consistency of the instrument during testing and to enable data normalization, a pooled QC was added to the sample tray after every 10 sample injections in both sequences. A blank sample was added after each QC sample.

Data preprocessing

To obtain a mass spectral feature table for all the samples, raw data obtained by LC–MS were transformed into mzXML files using ProteoWizard. MS-DIAL was then used to conduct spectral deconvolution and peak integration. Data collection in positive ion mode (ESI+) was kept for further analysis because more spectral information was contained than that from negative ion mode (ESI−).

The coefficient of variation (CV) was calculated for the pooled QC samples. Peaks meeting the criteria of CV less than 0.3 and an average intensity greater than 1e5 were identified. Only those peaks detected consistently were retained. The signal-to-noise (S/N) ratios of each feature were calculated by comparing the average intensity from samples with that of blanks. Features with a S/N ratio greater than 3 were retained as true signals. The data matrix that was filtered as described previously was utilized for subsequent analysis. Log transformation and auto-scaling were employed for feature normalization, while zero values were utilized for missing value imputation.

MOSS model development

In MOSS workflow (Figure 1), three-step feature selection procedures (from left to right) were implemented, involving different datasets, models, and feature sets, while maintaining the use of a consistent ROC validation metric for each stage, resulting in the respective outputs (from top to bottom). For the input dataset and features, a matrix with binary categorized samples (0 or 1) and multiple features (N) is required to go through a three-step feature selection. N is a predetermined constant representing the total feature number, while n denotes any variable within the range of 2 ~ N. In the first step, the categorized sample matrix is divided proportionally into a training set and a testing set in a ratio of 7:3 for 1000 permutations. For models, the training set is utilized as input for five different analytical models. We draw upon frequently utilized models from the literature to select and reconstruct a comprehensive and robust model for our analysis. This selection encompasses both univariate models, such as the Student t-test, and ROC analysis, and multivariate models. In particular, we emphasize the use of supervised models such as PLS-DA, RF, and SVM. All the features are ranked based on their importance indicators that corresponded to the aforementioned models, such as P-value, area under the curve (AUC), variable importance in the projection (VIP), Gini index, and SVM importance. These models are performed in R (Version 4.1.2) using pROC, mix0mics, randomForest, caret, and kernlab packages. In validation, within each model and for each permutation, we construct PLS-DA prediction models. These models are created individually using the testing set and begin with the top n features. The known categories and predicted categories obtained from the PLS-DA models are used to generate a ROC curve. The AUC is employed to evaluate the prediction performance. Therefore, a table of AUC values [1000 × (N−1)] is created and converted into a vector. This vector is visualized as a curve by computing the mean of the 1000 permutations of AUC for top n, where n ranges from 2 to N. The model with the highest peak AUC value was selected and the corresponding top n value was noted as M feature candidates for the subsequent step. Next, to comprehensively represent the entire dataset, the training set and the testing set are merged into a single data matrix. This matrix has M features that are selected from the top M features based on the average ranking in the 1000 permutations of the last step. The matrix is used as input for the selected model in the last step, which is internally tested via a PLS-DA model. The testing process begins with m, denoting any variable within the range of 2 ~ M. The known categories and predicted categories obtained from the PLS-DA models are used to generate a 3D-AUC scatter plot for m features (here in Figure 1, each figure’s scatterplot is condensed into a boxplot for dimensionality reduction), where m ranges from 2 to M. The model with the most significant improvement in AUC value compared with the top-ranked method and the smallest corresponding m value was chosen. The significance is defined by the performance top-ranked feature beyond the interquartile range (IQR) of the corresponding sampling features. The corresponding m value was recorded as L, denoting the optimal size for the final feature selection.

Figure 1 Schematic Workflow of MOSS—illustrating the step-by-step process of feature subset selection.

Finally, we explored all possible combinations of L features out of the total M features and identified the subset that yielded the highest AUC value. The features in this subset were ultimately chosen as the most effective ones for accurately categorizing the two groups.

Internal and external validation

Through the implementation of the five commonly used models in foodomics/metabolomics, robust and reliable models were generated, and the corresponding variable importance scores were effectively utilized to rank the features. To assess the performance of the models, we conducted both internal and external validation procedures. According to previous research [4], in terms of performance metrics for external validation, sensitivity, specificity, and accuracy are commonly utilized, followed by classification accuracy, which signifies the correct proportion that is predicted by the model, based on the gold-standard label. These findings are consistent with another study that reported the use of ROC and accuracy as standard methods for assessing the predictive capabilities of models [10]. Thus, in this study, internal validation was evaluated using the AUC metric, while external validation was evaluated using the accuracy metric. Specifically, for PLS-DA model validation, the Euclidean norm of the first two components was used to achieve both internal validation and external prediction.

RESULTS

To demonstrate the performance of MOSS in feature selection and its advantages over traditional top-ranked methods, we conducted a thorough re-analysis using a dataset generated by LC-MS on bourbon samples from our previous study [11]. Our investigation encompassed three key perspectives of American bourbon whiskey: age, taste, and proof. We customized the molecular feature subset size to select significantly differential features for each perspective, ensured a comprehensive understanding of their distinct characteristics, and showcased the intuitive performance of MOSS and its ability to uncover relevant and differentiating features within the bourbon dataset.

In previous studies, the identified compounds have demonstrated an accurate prediction of aged and unaged bourbon. However, differentiating between younger and older bourbon, specifically in terms of the age difference, remains challenging. Hence, this study first focuses on precisely this aspect—distinguishing between younger and older bourbon based on all unidentified features. A series of permutations (10, 100, 1000) was conducted as shown in Figure 2A–C. Increasing the number of permutations results in smoother curves, which helped prevent local optima. Figure 2C shows that the curves under 1000 times of permutation were smooth enough to identify the location of best performance while providing reliable and stable results thus were also used in the following analysis. The trend observed in all the curves has a gradual increase followed by a plateau as the number of features (denoted as n) increases until it reaches the total N = 1406 features. In our comparative analysis of five ranking models (Student’s t-test, ROC analysis, SVM, RF, and PLS-DA), the PLS-DA model, on average, reached the best AUC performance at 0.67 with M = 18 features. Those top 18 features ranked based on their average weight of 1000 permutations were kept for further selection (Table S2, see Supplementary Data available online at http://bib.oxfordjournals.org/). Next, the training and testing datasets were combined for further analysis. Meanwhile, PLD-DA models were built on these datasets 1000 times for each feature subset, composed of randomly selected m (range from 2 to 18) features from M = 18 feature candidates. To visualize the change in the model’s fitting ability based on the selected features, a 3D plot was created with the feature ranking number m, the sum of feature ranking number Σ(m), and AUC as the x-, y-, and z-axes, respectively. Figure 2D demonstrates that an increase in both the feature ranking number m and the sum of feature ranking number Σ(m) led to a corresponding increase in AUC. However, it was observed that the sum of feature ranking number Σ(m) did not have a significant impact on the change in AUC for a specific feature ranking number m. Thus, Figure 2E aggregates data points across each feature ranking number m, providing information on the changes in feature ranking number m and AUC. In addition, the internal validation performance of traditional methods that used the top m features was applied as a benchmark to compare with those randomly selected features from top M = 18 feature candidates. Figure 2E indicated that the use of the top 4–14 features provided significantly better performance compared with randomly selecting the same number of features from the previous set of 18. Interestingly, the study observed that when building 15 feature models, randomly selected features from the previously enlisted 18 outperformed the top 15 features. As a result, the study selected L = 15 features as the optimal subset and optimized them further for the next step. In the final step of MOSS, the study selected the best L = 15 features from the original M = 18. This involved examining all possible combinations of L = 15 features and identifying the one that yielded the best performance with AUC = 0.908 through internal validation (Figure 2F). The selected combination was then considered the final optimized subset for the study. The results of the random internal permutation test revealed that the final optimized subset had a significant ability (sensitivity = ~0.8) to differentiate between younger and older bourbon (Figure 2G). Furthermore, in the external validation, when used for predicting the ages of eight unknown external samples, this subset achieved an accuracy of 88% (Figure 2H).

Figure 2 Binary category of feature selection model Input. Performance comparison of five ranking models using AUC as a metric versus number of features to determine the number top N with best performance of classifying younger and older bourbon under permutation at (A) 10 times, (B) 100 times and (C) 1000 times. (D) 3D representation of random top m features’ AUC performance in relation to feature number and sum of feature numbers to determine the number M with best improved performance of classifying younger and older bourbon. (E) The AUC performance comparison between top m features (red dots) and the average of random selected m features (box plot). (F) The AUC performance of traversal combination of L features. (G) Internal and (H) external validation of selected feature on classifying younger and older bourbon.

Next, we performed the MOSS-based molecular feature selection that can assist in predicting the taste profiles of our bourbon samples. We selected 19 representative tastes to investigate their related molecular features. Figure 3A displays a colormap visualization of the taste profiles for all the bourbons. Notably, the study identified several common taste combinations, such as oaky/woody + spicy/peppery, caramel + vanilla, and caramel + fruity in our samples. Five or more samples containing a specific taste were selected for further analysis (Figure S1, see Supplementary Data available online at http://bib.oxfordjournals.org/, and Figure 3B), and the MOSS method was utilized for molecular feature selection that represents the taste of bourbon. Eventually, the spicy taste was used as an example for the selection of its distinctive molecular features as shown in Figure 3. In Figure 3B, the PLS-DA model outperformed the other four models and achieved the maximum AUC of 0.57 when using 11 features on average (Table S3, see Supplementary Data available online at http://bib.oxfordjournals.org/). We kept the top M = 11 features on average from 1000 times permutations and built PLS-DA models using all samples with randomly picked feature subsets with m features, ranging from 2 to 11. The three-dimensional scatter plot presented in Figure 3C illustrates a positive correlation between both the feature ranking number m and the sum of feature ranking numbers with the observed trend. However, no distinct pattern was observed within a certain feature ranking number, showing that the apex dominance is no longer significant in these top-ranking features. In Figure 3D, the performance of the top m features demonstrates an unsmooth curve by feature number. Conversely, the average performance of randomly selected m features from Top M exhibits a smooth change and outperforms the top m features most when m equals four. Figure 3E displays the performance of the PLS-DA models constructed using a traversed selection of L = 4 features from the M = 11 previously selected features. The best subset of four features resulted in an AUC of 0.898.

Figure 3 Conversion of concrete category into binary category for binary feature selection model input. (A) A colormap visualization for taste profiles of all the bourbons. (B) The comparison of five ranking models using AUC as a metric versus number of features to determine the number top N with best performance of classifying spicy/unspicy under permutation at 1000 times. (C) 3D representation of random top m features’ AUC performance in relation to feature number and sum of feature numbers to determine the number M with best improved performance of classifying spicy/unspicy. (D) The AUC performance comparison between top m features (red dots) and the average of random selected m features (box plot). (E) The AUC performance of traversal combination of L features.

Last but not least, we also performed MOSS-based molecular feature selection to find combinations that may represent the proof of the studied bourbon samples. It is important to recognize that proof is a continuous variable. Therefore, the initial step in this analysis involved identifying the optimal cut-off point for dividing the samples into higher and lower proof samples. To accomplish this, an internal validation (by permutation) was conducted on PLS-DA models using all available features and samples. The results of this validation within a wide range of proofs from 80 to 135.4 are displayed in Figure 4A, which demonstrates that a cut-off proof value of 98 provides the optimal threshold for determining higher and lower proof samples. Like the analysis of age and taste, PLS-DA still outperformed the other four models when distinguishing higher proof and lower proof bourbon. Notably, Figure 4B highlights that optimal performance occurred at a large feature ranking number of M = 54. Additionally, Figure 4C provides evidence that top m features consistently outperformed randomly selected feature subsets, thereby eliminating the need for further feature selection.

Figure 4 Conversion of continuous category into binary category for binary feature selection model Input. (A) The AUC performance of classifying higher/lower concentrations of bourbon via PLS-DA model changed by proof cut-offs. (B) Random top m features’ AUC performance in relation to feature number and sum of feature numbers to determine the number M with best improved performance of classifying higher proof and lower proof. (C) Performance comparison of five ranking models using AUC as a metric versus number of features to determine the number top N with best performance of classifying higher proof and lower proof under permutation at 1000 times.

DISCUSSION

Significance of unarbitrary feature selection

It is a common practice that multivariate statistical models of foodomics/metabolomics data generally help to reduce the dimension of variables by eliminating molecular features that had little impact on the classification models. Although statistical models can provide values indicating how important they are for the classification, a natural filtering cutoff is undecided. Thus, three distinct strategies have been typically employed for determining the cut-off. The first strategy involves setting a predetermined threshold value, such as VIP > 1.0 or 1.5 or a P-value < 0.05. The second strategy entails selecting a fixed number of top features, for instance, the top 15. Finally, the third strategy involves selecting a specific percentage of high-value features, such as the top 10%. Feature selection unsupported by empirical evidence may result in redundancy or exclusion of critical features, ultimately impacting the accuracy of the model. Hence, it is imperative to employ an unbiased and objective method for determining the cut-off point of the importance indicator and, consequently, the size of the variable subset after dimension reduction and filtering with a certain threshold.

Previous efforts on reducing feature dimensions

To partially address this issue, sparse extensions on commonly used foodomics/metabolomics classification models have been proposed, which add a penalty to the loading scores forcing some of the variables to have zero weight in the final model. In PLS-based models, this approach is called sparse partial least squares (SPLS), where a penalty term is added to the PLS objective function to encourage sparsity in the loading scores. Alternatively, in an RF model, a common approach to incorporating sparsity constraints is to use sparse random projections, which randomly select a subset of the input features and apply a sparse projection matrix to reduce the dimensionality of the data [12]. In SVM, sparse extensions can be applied using sparse regularization techniques, such as L1 or Lasso regularization. These penalties encourage sparsity in the weight vector of the SVM, resulting in a model that only uses a subset of the input features for prediction [13]. Among these models, the non-zero features are regarded as the important features. Despite reducing trophy features, the sparse extension does not give a highly productive and simplified list of significantly differential features.

Rational for multiple screening and permutation tests

To achieve both the productive computational selection and a simplified list of significantly differential features, within the MOSS algorithm, three-step feature selection was implemented. First, cross-validation was employed to identify the most promising subset candidates. Second, the subset candidates were further narrowed down in terms of their dimensions using the entire dataset. Finally, all possible combinations of features were traversed to determine the optimal feature subset. We chose five commonly employed feature-ranking models to accommodate both univariate and multivariate analysis approaches [10, 14–17] in the metabolomics field. As for the predictive model in the validation step, PLS-DA was chosen for its effectiveness in handling complex feature relationships and enhancing prediction accuracy in multivariate analysis, particularly for tasks with logical or linear relationships, offering simplicity compared to RF and SVM, which excel in handling non-linear relationships. The purpose of cross-validation is to assess the efficacy of a predictive model by training a portion of the data in the model while the remainder evaluates the model. In conjunction with cross-validation, permutation testing can reduce model bias by generating an empirical null distribution of the test statistic based on repeated shuffling of the labels for the data points. This process, repeated 1000 times, can help to identify sample outliers and ensure an ambiguous but inclusive candidate list for subsequent analysis. Additionally, employing features that are ranked on average can enhance the accuracy of the model’s performance assessment. By increasing the number of permutations from 10 to 1000, the performance curve becomes smoother, indicating observed changes in performance are more reliable and less affected by chance due to the permutation test. The performance curve resulting from permutation testing exhibits a characteristic pattern whereby model performance improves rapidly as features are added, reaches an optimal level at a certain number of permutations, and then gradually declines thereafter. During the initial phase of rapid improvement, the top-ranking features are utilized to their fullest extent, leading to significant gains in model performance. However, as additional features are added, the performance gains begin to diminish, suggesting that some of the added features contribute more redundancy than information to the model. This phenomenon underscores the importance of selecting features based on their contribution to the model rather than an arbitrary cut-off. Neglecting this principle can lead to arbitrary selection of subsets of features, resulting in suboptimal model performance. Sometimes, due to the limitations in sample size, individual differences may be just as apparent as group differences. As a result, cross-validation may yield some false positive features. Moreover, the primary objective is to identify the significantly differential features for the entire dataset. Therefore, in the second step, the entire dataset, comprising the training and testing sets, was utilized to construct the model using diverse combinations of feature candidates generated in the preceding stage. As the major advantage has already been attained in the last step, our objective now shifts from pursuing the highest performance to finding a harmonious equilibrium between performance and the number of selected biomarkers (molecular features), which corresponds to computational cost. In practical scenarios, fewer features make the tool more user-friendly. Therefore, we conducted a performance comparison between our random feature combination approach and the conventional feature selection method of choosing the top m features. The goal was to determine the minimum number of features at which our approach could outperform the top m feature selection method. The final step is to find the best subset with one out of m features. The best feature selection method usually traverses all the possible subsets and thus has a high computational cost; however, the first step and the second step greatly narrow the variable selection range and consequentially reduce the computational cost. Meanwhile, it is guaranteed that the optimal subset selected by traversing an ambiguous candidate list can achieve the balance between a minimal set of predictors and high predictive accuracy.

Different models of classification

Given the constraints due to the limited availability of samples, e.g. some groups only had several samples in flavor analysis, it may not be feasible for all groups to allocate samples for additional validation. Therefore, utilizing the AUC metric on the original samples would be considered the optimal approach to demonstrate the model’s performance. However, it is important to note that ROC is a binary model that necessitates binary group classification. To address this limitation, our study outlines several strategies to convert other variable types into binary variables, which can subsequently be inputted into the MOSS algorithm. As an example, we utilized three characteristics of bourbon—age, proof, and taste—to demonstrate how numerical (concrete or continuous), categorical variable types can be converted into binary variables. These strategies enable us to use a binary classification approach to evaluate the performance of the model and showcase its effectiveness even with limited sample availability. In this process, arbitrary selection of cut-off was also prevented by selecting the best cut-off of distinguishing between high proof and low proof among all the possible cut-offs.

Age is typically considered a continuous numerical variable; however, here, the age information provided is all integer years. Thus, we considered it as a concrete variable here. According to previous research, 0 and 4 years are commonly used cutoffs for categorizing unaged, younger, and older bourbon [9, 11, 18]. In terms of tastes, a categorical variable, each bourbon has multiple tastes, which makes tastes one of the features of tastes. To make better use of limited samples, each taste was used as an independent categorical variable to classify bourbon with and without this taste. There is no universally accepted threshold for determining the proof of bourbon. Typically, high-proof bourbon is defined as having an alcohol by volume (ABV) of 50% (proof = 100) or greater, while low-proof bourbon has an ABV of lower than 50%‘ or ’50% (proof = 100) or greater, while low proof bourbon has an ABV of 40% (proof = 80). Nonetheless, this range can fluctuate among various brands and even within a single brand’s product line. Due to the lack of an established cut-off, we chose to use each proof measurement as a temporary threshold. We employed a PLS-DA model, in conjunction with cross-validation, to classify the samples above and below each respective proof threshold. This approach helped us to avoid using an arbitrary cut-off and to treat proof as a continuous variable. The optimal proof of 98 lies between the typical ranges for high and low-proof bourbon. Not surprisingly, the PLS-DA model outperformed the other four models that were evaluated. In total, 52 features in the PLS-DA model were retained, contributing significantly to its superior performance. This number is substantially larger than the candidate features identified by age and taste analysis. This suggests that proof may have a more prominent impact on the selected features compared to age and taste. Subsequent analyses confirmed this hypothesis, as the top m features consistently outperformed the average performance of randomly selected m features, across the range of m features tested. However, we did not perform the next step of the analysis, as our objective was to identify the best method of improving the top m approach while minimizing the number of selected features.

Breakthrough/limitation/prospect

MOSS significantly enhances bourbon distinguishability based on age, flavor, and proof and holds the potential for authentic bourbon identification in the future. Moreover, this method’s applicability extends to other metabolomics and foodomics classification tasks, adding value to various analytical endeavors. Our MOSS algorithm offers several advantages over traditional approaches used to determine the cut-off for significant differential feature numbers. Typically, researchers rely on experience or default software settings to select a feature subset, but this can be arbitrary and vary significantly depending on the study aims, such as age, taste, and proof. MOSS provides a robust tool to avoid arbitrary feature subset selection and allows for customization of the subset size that will contribute best to the model classification. While acknowledging the significance of the top m method, we tactfully preserved its topmost advantage by integrating permutation analysis and, concurrently, preventing local optimal results caused by contingency on top features method and may also occur with forward feature selection and backward feature elimination methods. Additionally, it overcomes the computational challenges of best subset selection methods by narrowing down the features to an ambiguous candidate list, which saves a large amount of computational power. Overall, through the three-step feature selection conjugated with the permutation test, MOSS achieves a balance between robustness, computational cost, and minimal selected features. Despite the limitation of MOSS, which restricts its application to binary classification problems, our approach of transforming multiple categorical and numerical variables into binary variables represents an innovative approach to sample classification. Among all the methods that we used to support MOSS in this study, PLS-DA outperformed other statistical models. It effectively handles multicollinearity and high-dimensional data, unlike ROC analysis and t-tests that consider individual dimensions without accounting for correlations between predictor variables commonly seen in foodomics/metabolomics data. PLS-DA identifies important variables for classification, while RF and SVMs provide feature rankings. However, RF and SVM rankings may be challenging to interpret. RF measures metabolite importance by assessing the decrease in prediction accuracy upon removal, while SVM considers the weight in the decision boundary. PLS-DA’s transparency and interpretability make it widely used in the foodomics/metabolomics field.

CONCLUSION

In conclusion, the MOSS algorithm offers several advantages over traditional methods for determining significant differential feature numbers. By providing a robust tool for feature subset selection, MOSS avoids arbitrary selection and enables the customization of the subset size to contribute optimally to the model classification. The integration of permutation analysis helps prevent local optimal results caused by contingency on top features method and overcomes computational challenges of best subset selection methods. While MOSS is limited to binary classification problems, its approach to transforming multiple categorical and numerical variables into binary variables expands its range of applications. Future work may involve altering the evaluation criterion of the model to R2 to eliminate the issue of sample classification.

Key Points

Conventional statistical approaches often yield arbitrary feature selections, rankings, and cut-off thresholds, prompting the development of the Multiple and Optimal Screening Subset (MOSS) method.

MOSS performs a three-step approach, including the promising candidate subset selection, refining the subset, and determining the best subset.

Under the MOSS workflow, five widely recognized feature-selection models were rigorously compared in a case study focusing on the characteristics of bourbons, including age, taste, and proof.

By using MOSS, bourbon age prediction improved to 88%, and the AUC performance for assessing bourbon sweetness reached 0.898 with just four predictors.

MOSS is capable of striking a balance between minimal predictors and high predictive accuracy resonates, selecting an optimal subset of features in food omics/metabolomics studies.

Supplementary Material

Bourbon_II_supplementary_information-11152023_bbae046

FUNDING

OSU EHE dissertation fellowship (to R.X.); National Institute of General Medical Sciences of the National Institutes of Health (Award Number R35GM133510).

AUTHOR CONTRIBUTIONS

Rui Xu (Methodology, Software, Data curation, Writing—original draft, Visualization, Investigation, Writing—review & editing); Huan Zhang (Sample Collection, Methodology, Writing—review & editing); Michael W. Crowder (Conceptualization, Supervision, Writing—review & editing); and Jiangjiang Zhu (Conceptualization, Writing—original draft, Visualization, Investigation, Supervision, Writing—review & editing)

DATA AVAILABILITY

The mass spectra data have been deposited to MassIVE database (https://massive.ucsd.edu/ProteoSAFe/static/massive.jsp) with access # MSV000092689. The code of our analysis algorithm has been made available in GitHub with the following link: https://github.com/Xrrrr98784/MOSS.

Author Biographies

Rui Xu is an PhD student at The Ohio State University. Her research interest is bioinformatic methods and tools development.

Huang Zhang is a postdoctoral researcher at The Ohio State University. She has been working on chemical analyses of food/wine components with an interest of identifying molecular biomarkers.

Michael W. Crowder is a Professor in Chemistry and Biochemistry and the Dean of the Graduate School at Miami University. Dr Crowder has research interests in chemical analysis of wine/food ingredients for authentication and quality control purposes.

Jiangjiang Zhu is an Associate Professor at The Ohio State University. His primary interest is the development and applications of metabolomics strategies for human health and nutrition/food research.
==== Refs
References

1. Schiffman C , PetrickL, PerttulaK, et al.  Filtering procedures for untargeted LC-MS metabolomics data. BMC Bioinformatics  2019;20 :1–10.30606105
2. Karnovsky A , LiS. Pathway analysis for targeted and untargeted metabolomics. Methods Mol Biol  2020;2104 :387–400. 10.1007/978-1-0716-0239-3_19.
3. Li S , ParkY, DuraisinghamS, et al.  Predicting network activity from high throughput metabolomics. PLoS Comput Biol  2013;9 :e1003123.23861661
4. Considine E , ThomasG, BoulesteixA-L, et al.  Critical review of reporting of the data analysis step in metabolomics. Metabolomics  2018;14 :1–16.29249916
5. Gu H , PanZ, XiB, et al.  Principal component directed partial least squares analysis for combining nuclear magnetic resonance and mass spectrometry data in metabolomics: application to the detection of breast cancer. Anal Chim Acta  2011;686 :57–63.21237308
6. Deng L , GuH, ZhuJ, et al.  Combining NMR and LC/MS using backward variable elimination: metabolomics analysis of colorectal cancer, polyps, and healthy controls. Anal Chem  2016;88 :7975–83.27437783
7. Zhu W . P < 0.05, < 0.01, < 0.001, < 0.0001, < 0.00001, < 0.000001, or < 0.0000001 …. J Sport Health Sci  2016;5 :77–9.30356881
8. Kennedy-Shaffer L . Beforep < 0.05 to Beyondp < 0.05: using history to contextualizep-values and significance testing. Am Stat  2019;73 :82–90.31413381
9. Yang K , SomogyiA, ThomasC, et al.  Analysis of barrel-aged Kentucky bourbon whiskey by ultrahigh resolution mass spectrometry. Food Anal Methods  2020;13 :2301–11.
10. Ghosh T , ZhangW, GhoshD, et al.  Predictive modeling for metabolomics data. Methods Mol Biol  2020;2104 :313–6. 10.1007/978-1-0716-0239-3_16.
11. Xu R , ChenL, ZhangH, et al.  Characterizing bourbon whiskey via the combination of LC-MS and GC-MS based molecular fingerprinting. Food Chem  2023;423 :136311.37167670
12. Li P , HastieTJ, ChurchKW. Very sparse random projections. In: Proceedings of the 12th ACM SIGKDD international conference on Knowledge discovery and data mining. Association for Computing Machinery (ACM), New York, NY, USA, 2006, 287–96.
13. Sun Z , QiaoY, LelieveldtBP, et al.  Integrating spatial-anatomical regularization and structure sparsity into SVM: improving interpretation of Alzheimer's disease classification. Neuroimage  2018;178 :445–60.29802968
14. Gromski PS , MuhamadaliH, EllisDI, et al.  A tutorial review: metabolomics and partial least squares-discriminant analysis–a marriage of convenience or a shotgun wedding. Anal Chim Acta  2015;879 :10–23.26002472
15. Liebal UW , PhanAN, SudhakarM, et al.  Machine learning applications for mass spectrometry-based metabolomics. Metabolites  2020;10 :243.32545768
16. Xia J , BroadhurstDI, WilsonM, WishartDS. Translational biomarker discovery in clinical metabolomics: an introductory tutorial. Metabolomics  2013;9 :280–99.23543913
17. Saccenti E , HoefslootHC, SmildeAK, et al.  Reflections on univariate and multivariate analysis of metabolomics data. Metabolomics  2014;10 :361–74.
18. Collins TS , ZweigenbaumJ, EbelerSE. Profiling of nonvolatiles in whiskeys using ultra high pressure liquid chromatography quadrupole time-of-flight mass spectrometry (UHPLC–QTOF MS). Food Chem  2014;163 :186–96.24912715
