Metabolomics-based approaches for food authentication and traceability

Rebeca Tatiana de Souto Santos Metabolomics-Based Approaches For Food Authentication and Traceability Universidade do Minho Escola de Engenharia April 2023 Rebeca Tatiana de Souto Santos Metabolomics-Based Approaches For Food Authentication and Traceability Minho | 2023 U Universidade do Minho Escola de Engenharia Rebeca Tatiana de Souto Santos Metabolomics-Based Approaches For Food Authentication and Traceability Doctoral Thesis Doctoral Program in Chemical and Biological Engineering Work developed under the supervision of: Professor Miguel Rocha Professor Marcelo Maraschin April 2023 DIREITOS DE AUTOR E CONDIÇÕES DE UTILIZAÇÃO DO TRABALHO POR TERCEIROS Este é um trabalho académico que pode ser utilizado por terceiros desde que respeitadas as regras e boas práticas internacionalmente aceites, no que concerne aos direitos de autor e direitos conexos. Assim, o presente trabalho pode ser utilizado nos termos previstos na licença abaixo indicada. Caso o utilizador necessite de permissão para poder fazer um uso do trabalho em condições não previstas no licenciamento indicado, deverá contactar o autor, através do RepositoriUM da Universidade do Minho. Licença concedida aos utilizadores deste trabalho Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International CC BY-NC-SA 4.0 https://creativecommons.org/licenses/by-nc-sa/4.0/deed.en ii Acknowledgements I would first like to thank European Social Fund under the scope of Norte2020 - Programa Operacional Regional do Norte for doctoral advanced training I was awarded (call NORTE-69-2015-15), without which this work would not be possible. I also thank my main host institution, the Centre of Biological Engineering at the University of Minho. I also wish to express my gratitude to my supervisors, Professor Miguel Rocha and Professor Marcelo Maraschin. I thank them for offering me this opportunity and for all their support. I am immensely grateful for the guidance and feedback they provided throughout this project. I would also like to thank Professor Eugénio Campos Ferreira for his availability, kindness, and opportunity to develop work in the BIOSYSTEMS group. I also thank to the collaborations for the kindly supplied of the metabolomics data used in this dissertation, mainly the Enology Laboratory of INIAV - Polo de Inovação de Dois Portos, the University of Aveiro, Portugal, the federation “Federação das Associações de Apicultores e Meliponicultores de Santa Catarina” (FAASC) Brazil, the company “Empresa de Pesquisa Agropecuária e Extensão Rural de Santa Catarina”(EPAGRI) Brazil and the Federal University of Santa Catarina (UFSC) Brazil. I would like to thank all of my colleagues from the BISBII research group for the friendly and collaborative working environment, and for their useful comments and suggestions, especially Sara Cardoso for her technical help with the R package specmine and Delora Baptista, Telma Afonso, Patrica Dias for their friendship till today. A special thanks to Saulo Rodrigues for the words of support and friendship that started with this work and that will continue for a lifetime. Finally, a special thanks to my close friends Cristina Costa, Manuela Tiago, Carlos Alberto, my beloved Parents and Brother and my dear boyfriend, Pedro Tiago Evangelista, without a doubt your support, encouragement, patience, kindness, friendship and love have been precious in last years. To you, I dedicate the commitment of these years. iii STATEMENT OF INTEGRITY I hereby declare having conducted this academic work with integrity. I confirm that I have not used plagiarism or any form of undue use of information or falsification of results along the process leading to its elaboration. I further declare that I have fully acknowledged the Code of Ethical Conduct of the Universidade do Minho. iv Resumo Abordagens Metabolómicas para Autenticação Alimentar e Rastreabilidade Nos últimos anos, a procura do consumidor por produtos alimentares naturais autênticos aumentou consideravelmente, assim como a ocorrência de eventos de adulteração. Paralelamente, a conscientização sobre a responsabilidade de produtores e vendedores no combate à fraude alimentar tem vido a crescer. O uso de ferramentas estatísticas (ST) e algoritmos de aprendizado de máquina (ML), em conjunto com abordagens metabolómicas, podem ajudar numa ampla gama de aplicações na área da autenticação alimentar. Os modelos de ML permitem identificar e entender quais os fatores com maior influência na autenticação alimentar, tais como as condições biológicas e ambientais. A Ressonância Magnética Nuclear (NMR) e metodologias metabolómicas espectrais são técnicas analíticas avançadas e abrangentes amplamente utilizadas na autenticação de alimentos. No entanto, existe a necessidade de construir modelos mais precisos de forma a prever a autenticidade dos produtos alimentares, garantindo a sua segurança e qualidade para o bem-estar dos consumidores e da economia. Esta dissertação foca-se no uso da análise estatística multivariada (MSA) em dados metabolómicos de produtos alimentares naturais de alto valor, como o vinho e o mel, com recurso à aplicação de modelos de ML a vários problemas de classificação, incluindo a determinação de origens botânicas e geográficas. Modelos de ML foram construídos para discriminar e prever variáveis relacionadas com os problemas de autenticação dos alimentos estudados, incluindo a previsão de diferentes influências climáticas e a previsão da idade de armazenamento com base em análises de dados metabolómicos. Diferentes estratégias de pré-processamento, seleção de variáveis, técnicas de redução de dimensionalidade e métricas de avaliação foram testadas. Foi avaliado o impacto de diferentes métodos utilizados no desempenho desses modelos. Com base nessas análises, foram identificados os modelos com melhor desempenho preditivo. A interpretação dos modelos foi realizada comprovando que com o uso de modelos de ML é possível associar características químicas importantes com a classificação de alimentos naturais no processo de autenticação. ST e ML com tratamentos de pré-processamento adequados e métricas de avaliação são um pipeline adequado para diferentes problemas em autenticação de alimentos, como demonstrado neste trabalho. Apesar de algumas limitações, tais como, o tamanho dos conjuntos de dados dos nossos estudos, as nossas descobertas são úteis para orientar o desenvolvimento de novas abordagens usando a metabolómica para diversos problemas existentes em autenticação alimentar conjugando com ST e ML em amostras de vinho e mel. Para realizar análises de dados metabolómicos, diferentes abordagens MSA e de ML foram usadas, testadas e validadas usando o package da linguagem de programação R specmine , um software desenvolvido pelo grupo de investigação, e um conjunto de bibliotecas da linguagem de programação Python . Palavras-chave: Autenticação de Vinho e Mel, Metabolómica, Ferramentas Estatísticas e Aprendizado Máquina. v Abstract Metabolomics-Based Approaches For Food Authentication and Traceability In recent years, the consumer demand for authentic natural food products has increased concurrently with an expand of adulteration events. At the same time, the awareness regarding the responsibility of retailers and producers to combat food fraud has also grown. The use of statistical tools (ST) and machine learning (ML) algorithms, in conjunction with highthroughput metabolomic-based approaches, can help in a wide range of applications in food authentication. ML models allow to identify and understand the most influential factors in food product authentication, such as biological and environmental conditions. Nuclear Magnetic Resonance (NMR) and other spectral metabolomic methodologies are advanced and comprehensive analytical techniques widely used in food authentication. Nonetheless, there is a need to construct more accurate models to predict the authenticity of food products to ensure their safety and quality for the well-being of the consumers and the economy. This dissertation is focused on the use of multivariate statistical analysis (MSA) of metabolomics data from high-value natural food products, as wine and honey, addressing the application of ML models to several classification problems, including the botanical and geographical origins. ML models were also built to discriminate and predict variables related to food authentication problems studied, including the prediction of different climatic influences on production samples, and the prediction of the age of storage based on metabolomics data analyses. Different preprocessing strategies, feature selection, dimensionality reduction techniques, and evaluation metrics were tested. The impact of different methods used on the performance of these models was evaluated. Based on these analyses, the models with better predictive performance were identified. Models interpretability was performed showing that by using ML models it is possible to associate important chemical features with food authentication classification. ST and ML with adequate preprocessing treatments and evaluation metrics are a suitable pipeline for different food authentication problems, as shown in this work. Despite some limitations of the datasets sizes of our studies, our findings are useful to guide the development of new metabolomics based approaches for different food authentication issues using ST and ML algorithms in wine and honey samples. To perform metabolomic-data analyses, different MSA and ML approaches were used, tested, and validated using the open source R package specmine , a software developed by the host group, and a set of open source Python libraries. Keywords: Wine and Honey Authentication, Metabolomics, Statistical and Machine learning tools. vi Contents List of Figures xi List of Tables xvii Acronyms xxi 1 Introduction 1 1.1 ContextandMotivation .............................. 1 1.2 ResearchObjectives................................ 3 1.3 Outline...................................... 4 2 Metabolomics 6 2.1 Definitions and general workflow . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.2 SamplePreparation................................ 9 2.3 Data Acquisition: analytical techniques . . . . . . . . . . . . . . . . . . . . . . 9 2.3.1 MassSpectrometry............................ 10 2.3.2 Nuclear Magnetic Resonance . . . . . . . . . . . . . . . . . . . . . . . 10 2.3.3 Infrared Spectroscopy . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.3.4 Ultraviolet-visible Spectroscopy . . . . . . . . . . . . . . . . . . . . . . 11 2.3.5 RamanSpectroscopy........................... 11 2.4 DataPreprocessing................................ 11 2.4.1 Missing values and Outliers . . . . . . . . . . . . . . . . . . . . . . . . 12 2.4.2 SpectraTreatment ............................ 12 2.4.3 Mean-centering.............................. 14 2.4.4 Scalingmethods............................. 14 2.4.5 First and second-order derivatives . . . . . . . . . . . . . . . . . . . . . 14 2.4.6 Data transformation . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.5 DataAnalysis................................... 15 vii List of Figures 29 The 20 most important features are shown and ordered according to their importance to the best model classification for each production region. The Shap values and their impact on model output is shown for each feature. Each point represents a sample in the sub-dataset, and the color of the points represents the value of a particular feature for that sample. A: SHAP values and impact of feature’s on model of LR to classify 1A. B: SHAP values and impact of feature’s on model of LR to classify 1B. . . . . . . . 110 30 Top 20 most important features ranked by mean absolute SHAP values. The features are ordered according to their importance of the best model classification for each each production region. The Mean Shap values and their average impact on model output is shown for each feature. A: SHAP values and the average impact of feature’s on model of PLS-DA to classify 2C. B: SHAP values and the average impact of feature’s on model of RF to classify 3B. C: SHAP values and the average impact of feature’s on model of LR to classify 3C........................................... 110 31 The 20 most important features are shown and ordered according to their importance to the best model classification for each production region. The Shap values and their impact on model output is shown for each feature. Each point represents a sample in the sub-dataset, and the color of the points represents the value of a particular feature for that sample. A: SHAP values and impact of feature’s on model of PLS-DA to classify 2C. B: SHAP values and impact of feature’s on model of RF to classify 3B. C: SHAP values and the impact of feature’s on model of LR to classify 3C. . . . . . . . . . . . . 111 32 Top 20 most important features ranked by mean absolute SHAP values. The features are ordered according to their importance of the best model classification for each each production region. The Mean Shap values and their average impact on model output is shown for each feature. A: SHAP values and the average impact of feature’s on model of SVC to classify 4A. B: SHAP values and the average impact of feature’s on model of PLS-DA to classify 4B. C: SHAP values and the average impact of feature’s on model of SVC to classify 5. .......................................... 111 33 The 20 most important features are shown and ordered according to their importance to the best model classification for each production region. The Shap values and their impact on model output is shown for each feature. Each point represents a sample in the sub-dataset, and the color of the points represents the value of a particular feature for that sample. A: SHAP values and impact of feature’s on model of SVC to classify 4A. B: SHAP values and impact of feature’s on model of PLS-DA to classify 4B. C: SHAP values and the impact of feature’s on model of SVC to classify 5. . . . . . . . . . . . . 112 xiv List of Figures 34 Top 20 most important features ranked by mean absolute SHAP values. The features are ordered according to their importance of the best model classification for flora origin. The Mean Shap values and their average impact on model output is shown for each feature. A: SHAP values and the average impact of feature’s on model of PLS-DA to classify “Uva do Japão”. B: SHAP values and the average impact of feature’s on model of SVC to classify Other. C: SHAP values and the average impact of feature’s on model of SVC toclassify“Silvestre”.................................. 113 35 The 20 most important features are shown and ordered according to their importance to the best model classification for flora origin. The Shap values and their impact on model output is shown for each feature. Each point represents a sample in the sub-dataset, and the color of the points represents the value of a particular feature for that sample. A: SHAP values and impact of feature’s on model of PLS-DA to classify “Uva do Japão”. B: SHAP values and impact of feature’s on model of SVC to classify Other. C: SHAP values and the impact of feature’s on model of RF to classify “Silvestre”. . . . . . . . . 113 36 Synopsis of the results whit higher performance and the data analysis on 1HNMR spectra for the regions honeys productions and flora. The squares present the classification model and the metric score of the best model performance to flora origins. And the colored lines represent the links between flora and production regions. . . . . . 115 37 Principal Components Analysis (PCA) 2D scores plot (PC1, PC2) from FT-IR spectra of red wine musts samples of varieties Aragonez, Syrah and Touriga Nacional (Touriga Nacional). Explain variance of Principal Componets 1 and 2:PC1 = 69.8% PC2 =21.6%. ....................................... 144 38 t-distributed Stochastic Neighbor Embedding (t-SNE) from FT-IR spectra of red wine musts samples of varieties Aragonez, Syrah and Touriga Nacional (Touriga Nacional) ...................................... 145 39 Dendrogram plot of the hierarchical clustering analysis (HCA), with Euclidean distance between samples, from FT-IR spectra of red wine musts samples of varieties Aragonez, Syrah and Touriga Nacional (Touriga Nacional). . . . . . 145 40 Principal Components Analysis (PCA) 2D scores plot (PC1, PC2) from FT-IR spectra of white wine musts samples of varieties Arinto and Viosinho. Explain variance of Principal Components 1 and 2: PC1 = 43.4% PC2 = 30.4%. . . . . . . . . . . . . . 146 41 t-distributed Stochastic Neighbor Embedding (t-SNE) from FT-IR spectra of white wine musts samples of varieties Arinto and Viosinho. .............. 146 42 Dendrogram plot of the hierarchical clustering analysis (HCA), with Euclidean distance between samples, from FT-IR spectra of white wine musts samples of varieties Arinto and Viosinho. ........................... 147 xv List of Figures 43 Principal Components Analysis (PCA) 2D scores plot (PC1, PC2) from 1H-NMR spectrum of honey samples of all SC Regions. Explain variance of Principal Components 1 and 2: PC1 = 52%; PC2 =29%;PC3=0.5%. . . . . . . . . . . . . . . . . . . . 148 44 t-distributed Stochastic Neighbor Embedding (t-SNE) from 1H-NMR spectrum of honey samples of all SC Regions. ......................... 149 45 Dendrogram plot of the hierarchical clustering analysis (HCA), with Euclidean distance between samples, from 1H-NMR spectrum of honey samples of all SC Regions. ....................................... 149 46 Principal Components Analysis (PCA) 2D scores plot (PC1, PC2) from 1H-NMR spectrum of honey samples of Flora origins. Explain variance of Principal Components 1 and 2: PC1 = 52%; PC2 = 29%; PC3 = 0.5%. . . . . . . . . . . . . . . . . . . . . 150 47 t-distributed Stochastic Neighbor Embedding (t-SNE) from 1H-NMR spectrum of honey samples of Flora origins. .......................... 151 48 Dendrogram plot of the hierarchical clustering analysis (HCA), with Euclidean distance between samples, from 1H-NMR spectrum of honey samples of Flora origins. ....................................... 152 xvi List of Tables 1Metabolite Annotation and Metabolite Identifications tools. .......... 16 2Available free tools for metabolomics and spectral data. ............ 23 3Metabolomics Approaches on Food Authentication. ............... 32 4Metabolomics Approaches for Wine Authentication and Traceability. ..... 36 5Metabolomics Approaches for Honey Authentication and Traceability. . . . . 38 6Metabolites and their 1H chemical shifts annotated by 500 MHz 1H-NMR with higher scores for frequency, organism and solvent. ............... 44 7ANOVA results - peaks with the lowest corrected p-values (< 0.01). The first column shows the peaks, the second one the respective p-value, the third column includes the corrected p-value with False Discovery Rate (FDR) method, and the final column the result of the Tukey’s HSD test, which consists of the pairs of harvest years that were significantly different in terms of means for each peak. . . . . . . . . . . . . . . . . . . . . . . . 51 8ANOVA results - amino acids with the lowest corrected p-values (< 0.05). The first column contains the considered amino acids, the second one the respective p-value, the third column includes the corrected p-value with False Discovery Rate (FDR) method, and the final column the result of Tukey’s HSD test, which consists of the pair of harvest years that were significantly different in terms of means for each amino acid concentration. 52 9ML models results for the 1H-NMR data and the amino acids concentration with cross-validation accuracy. The first column contains the considered input data as features of the ML model, the second one the applied ML model, the third one the result with the accuracy of the correspondent ML model computed using cross-validation, and the last column the features (1H-NMR resonances - ppm and amino acids) most relevant, with those marked with bold type indicating those shared with ANOVA results. . . . . . . . . . . . 55 10 Climatic Data by each harvest year ........................ 58 xvii List of Tables 11 Results for linear regression analysis of the climatic influence over 1H-NMR data. The first column contains the considered peaks, the second column presents the influence type, followed by the respective corrected p-value(<= 0,0001), adjusted R2(> 50%), and coefficient of the LR model. Bold type indicates 1H-NMR resonances with highest adjusted R2. ......................................... 58 12 Regions of 1H-NMR resonances (ppm) with predictive capability to distinguish samples selected by ANOVA and ML models. .................. 63 13 ANOVA results. Wavenumber (cm-1) with the lowest corrected p, values (< 0.01). The first column shows the t-test, which consists of the pairs of grapevine musts varieties that were significantly different in terms of means for each feature, the second one the respective sprectroscopic bands, the third column includes the chemical groups related to thewavenumbersidentified............................... 76 14 Performance scores of the best classification models for the grapevine musts varieties studied. All models had pre-process methods(SNV plus Scaling) with the application of Feature Selection (Drop of constant features with 90% of correlation). The MCC evaluation metric followed a threshold based on the calculated PR Threshold mean. Note: PR-AUC: Area Under the Precision-Recall Curve; ROC-AUC: Area Under the Receiver Operating Characteristic Curve; MCC:Matthews Correlation Coefficient; BA: Balanced Accuracy; PR Threshold: Precision-Recall mean Threshold; SD: Standard Deviation. . . . . . . . . 80 15 Comparative analysis between Statistical, Machine Learning approach and Shapley values results. Note: Arg: Aragonez; TN: Touriga Nacional; SY: Syrah; LR: Logistic Regression; SVC: Support Vector Classification; PLS-DA: Partial Least Squares-Discriminant Analysis; (-): Negative impact on the model output; (+): Positive impact on the model output. 87 16 Hyper-parameters used by the distinct models. ................. 93 17 Geographical information of the honey samples. ................ 99 18 Flora information of the honey samples. Uva mix samples are multiflora honeys with presence of “Uva do Japão” flora. Other honey samples are multiflora honeys with different types of flora as cinnamon, fruit trees, native trees, ”Canudo de Pito”tree, ”Bracantinga”tree, and other plants from each region. . . . . . . . . . . . . . . . . . . . . . . . . . . 99 19 Feature with highest variance and the correspondent correlated features group. First column presents the feature with the higher variance correlated to a set of peaks given a predefined threshold. These features are the selected to the new sub-dataset used at all models classification. And these features were used by Shapley method to interpret their impact on each model classification. Second column are the correspondent correlated features to each feature. Note: * are features with high variance and with lower threshold correlation (0.90)withotherfeatures................................ 103 xviii List of Tables 20 Performance scores for best models classification for all Regions. All models had pre-process methods with the application of Feature Selection (Drop of constant features with 90% of correlation). The MCC evaluation metric followed a threshold based on the calculated ROC Threshold mean. Rows with background color gray are the regions with higher classification models. Note: PR-AUC: Area Under the Precision-Recall Curve; ROCAUC: Area Under the Receiver Operating Characteristic Curve; MCC: Matthews Correlation Coefficient; BA: Balanced Accuracy; ROC Threshold: Receiver Operating Characteristic. . 104 21 Performance scores for best models classification for all Flora honeys. All models had pre-process methods with the application of Feature Selection (Drop of constant features with 90% of correlation). The MCC evaluation metric followed a threshold based on the calculated ROC Threshold mean. Rows with background color gray are the regions with higher classification models. Note: PR-AUC: Area Under the Precision-Recall Curve; ROCAUC: Area Under the Receiver Operating Characteristic Curve; MCC: Matthews Correlation Coefficient; BA: Balanced Accuracy; ROC Threshold: Receiver Operating Characteristic. . 105 22 Best seven Mean absolute SHAP values and their impact on model output of the best models classification for Geographic Regions. Note: (-): Negative impact on the model output; (+): Positive impact on the model output. . . . . . . . . . . . . . . . . 109 23 Best seven Mean absolute SHAP values and their impact on model output of the best models classification for Flora groups. Note: (-): Negative impact on the model output; (+): Positive impact on the model output. . . . . . . . . . . . . . . . . . . . 112 24 Hyper-parameters used by the distinct models. ................. 120 25 Feature with highest variance and the correspondent correlated features group for Red varieties. First column presents the feature with the higher variance correlated to a set of FT-IR bands given a predefined threshold. This features are the selected to the new sub-dataset used at all red models classification. And these features were used by Shapley method to interpret their impact on each wine grape must variety model classification. Second column are the correspondent correlated features to each feature. Note: 1015*, 1123*, 1127*, 1130* are the features computed by Shapley value with zero correlations to other features, in our work designated as unique features. The wasenumbers at 930**, 980**,1080** and 2639** were features computed by Shapley values for red and white varieties. . . . 153 xix List of Tables 26 Feature with highest variance and the correspondent correlated features group for White varieties. First column presents the feature with the higher variance correlated to a set of FT-IR bands given a predefined threshold. This features are the selected to the new sub-dataset used at all white models classification. And these features were used by Shapley method to interpret their impact on each wine grape must variety model classification. Second column are the correspondent correlated features to each feature. Note: 2924* is the only feature computed by Shapley value with zero correlations to other features, in our work designated as unique features. The wasenumbers at 930**, 980**, 1080**, 2639** were features computed by Shapley values for red and white varieties. . . . . . . . . . 154 27 Supplemental Table: Performance scores for all models classification for all wine musts varieties. All models had pre-process methods (SNV plus Scaling) with the application of Feature Selection (Drop of constant features with 90% of correlation). Note: PRAUC: Area Under the Precision-Recall Curve; ROC-AUC: Area Under the Receiver Operating Characteristic Curve; MCC: Matthews Correlation Coefficient; BA: Balanced Accuracy; PR Threshold: Precision-Recall Threshold; SD: Standard Deviation. . . . . . . . . . . . . . 155 28 Supplemental Table: Performance scores for all models classification for all Regions. All models had pre-process methods with the application of Feature Selection (Drop of constant features with 90% of correlation). Note: PR-AUC: Area Under the Precision-Recall Curve; ROC-AUC: Area Under the Receiver Operating Characteristic Curve; MCC: Matthews Correlation Coefficient; BA: Balanced Accuracy; PR Threshold: Precision-Recall Threshold; SD:StandardDeviation. ............................... 156 29 Supplemental Table: Performance scores for all models classification for all Flora honeys. All models had pre-process methods with the application of Feature Selection (Drop of constant features with 90% of correlation). Note: PR-AUC: Area Under the Precision-Recall Curve; ROC-AUC: Area Under the Receiver Operating Characteristic Curve; MCC: Matthews Correlation Coefficient; BA: Balanced Accuracy; PR Threshold: PrecisionRecall Threshold; SD: Standard Deviation. . . . . . . . . . . . . . . . . . . . . . . . 157 xx Acronyms 1H-NMR Hydrogen Proton Nuclear Magnetic Resonance 1D-NMR one-dimensional NMR spectrum 2D-NMR two-dimensional NMR spectrum Ala Alanine ANN Artificial Neural Network ANOVA Analysis of Variance BA Balanced Accuracy CE Capillary Electrophoresis CS Cabernet Sauvignon wines Cys Cysteine D2O deuterium oxide D-UPLS least-squares and discriminant unfolded partial least-squares dGMP 2-deoxyguanosine monophosphate EU European Union FIR Far-Infrared FSA Food Standards Agency FT-ICR Fourier transform ion cyclotron resonance FT-IR Fourier-transform infrared spectroscopy xxi GBC Gradient Boosting Classifier GC Gas Chromatography Glu Glutamic acid Gly Glycine HCA Hierarchical Clustering Analysis Hpro Hydroxyproline IM Ion Mobility IR Infrared Spectroscopy Iso Isoleucine KNN K-Nearest Neighbors LC Liquid Chromatography LDA Linear Discriminant Analysis Leu Leucine LnR Linear Regression model LR Logistic Regression model Lys Lysine MCC Matthews Correlation Coefficient MIR Mid-Infrared ML machine learning MRL maximum residue limits MS Mass Spectrometry MSA Multivariate Statistical Analysis MSC Multiplicative Scatter Correction NIR Near-Infrared NMR Nuclear Magnetic Resonance xxii OIV International Organisation of Vine and Wine OPLS-DA Orthogonal Projections to Latent Structures Discriminant Analysis PCA-DA Principal Component Analysis-Discriminant Analysis PCA Principal Component Analysis PDO Protected Designations of Origin PGI Protected Geographical Indication Phe Phenylalanine PLS-DA Partial Squares-Discriminant Analysis PLS-r Partial Least Squares Regression PR Threshold Precision-Recall threshold PR-AUC Area Under the Precision-Recall Curve Pro Proline RF Random Forest Classifier RMSE Root Mean Square Error ROC-AUC Area Under the Receiver Operating Characteristic curve SD Standard Deviation Ser Serine SHAP Shapley Additive exPlanations SIMCA Independent Modelling by Class Analogy SNV Standard Normal Variate SPME-GC-MS Solid-Phase Microextraction followed by Gas Chromatography-Mass Spectrometry ST statistical tools SVC Support Vector Classifier SVM Support Vector Machines t-SNE t-distributed Stochastic Neighbor Embedding xxiii 2 Metabolomics 2.1 Definitions and general workflow Different analytical techniques for the suitability of food authentication processes have been assessed throughout the years. Chromatographic and molecular techniques are the most common techniques used and presented in studies of food authentication. Nonetheless, there are non-invasive techniques able to achieve good performances which are being used for the food authentication process, like: isotopic, vibrational, UV-vis, fluorescence spectroscopy, elemental techniques, and nuclear magnetic resonance presented in Figure 1 [10]. This chapter describes briefly main analytical techniques employed for food authentication with emphasis in the methods used in this work, regarding wine and honey authentication. Figure 1: Percentage of publications distributed over different analytical techniques applied to food authentication until 2016. Adapted from [10]. The field of Metabolomics is one of the main omics areas, and consists of the exhaustive study of the whole metabolome composition of a particular system or organism and their environmental interactions, focusing on small chemical molecules involved in cell processes called metabolites [9], [11]. 6 2.1. DEFINITIONS AND GENERAL WORKFLOW There are a considerable variety of metabolites, such as carbohydrates, lipids, proteins, amino acids, amines, steroids, phenolic compounds, carotenoids, alkaloids, or volatile compounds, among others. Determining the authenticity of food products could involve different methodologies, depending on the purpose and the extension of the food analysis [11], [9], [10]. Therefore by choosing the most suitable metabolomics approach(s) and method(s), it is possible to clear up different issues and unsolved problems, such as the ones in food authentication [9]. The constant evolution and the new developments in analytical techniques, instrumentation, analytical software, statistical methods, and computational techniques allow to accelerate or improve data collection, data analysis, and data interpretation [12], [13]. Nowadays, the analytical techniques used in metabolomics approaches are robust, sensitive, reproducible, and more economical, meaning that it is possible to validate the results of more studies and also to standardize methodologies [12], [13]. Metabolomics analyses have been generally classified based on the selected methods used to study the metabolites. The targeted analysis focus on a specific group of metabolites with identification and quantification. Untargeted metabolomics studies follow a detection approach to a non-specific group of metabolites, to obtain patterns or fingerprints without the need for identification or quantification, as shown in Figure 2 [14], [9]. Based on the data manipulation, metabolomics studies can be categorized into four to four conceptual approaches: • Targeted analysis: follow the identification and quantification of a small set of known metabolites related to a specific metabolism interaction, also called targets. This approach uses analytical technique(s) detection of interest compounds with the best performance. • Metabolite profiling: the analysis focus on a large set of compounds related to a group of specific metabolites. The metabolites of this analysis do not require quantification and do not need to be identified. Semi-quantification is the procedure commonly performed. • Metabolic footprinting: an approach with a study focus on the external and/or secreted metabolites by the cells. This analysis can follow a targeted procedure at a specific group of metabolites or through the spectra not directed to any selective metabolites group. • Metabolic fingerprinting: This approach delivers a metabolites ”signature” or fingerprint spectra of a sample of interest. The generated data is then compared to a large set of other samples to screen for differences. This approach does not focus on providing information about specific metabolites and is common used in classification problems based on its spectrum, as in food authentication studies. Considering the specific objective of the analysis some metabolomics studies can be classified as 7 CHAPTER 2. METABOLOMICS discriminative, informative, and/or predictive as shown in Figure 2[15], which can be achieved by two most common approaches: the Target and the Untargeted as presented in Figure 2[15]. In Food authentication studies when spectral signals are able to discriminate between samples, then the metabolites can be identified and the biological relevance of that compound elucidated, saving valuable analysis time. It is commonly used in forensics, among other fields, having been applied for instance in the discrimination of food quality and in food authentication. This approach will be emphasized throughout this dissertation. Figure 2: More common Metabolomics Purposes and based Approaches. Adapted from [9]. An example of a discriminative purpose using metabolomics approaches is assigning specific labels and certifications to food products, confirming their authenticity [10]. So depending on the study purpose (Predictive, Discriminative or Informative) a different metabolic-based approach is followed [16]. In summary, metabolomics approaches have the potential for solving many problems in food authentication, namely in the issues regarding food origin determination (botanical and geographical), traceability, and testing of possible adulterations. Although the existence of different metabolomics approaches, the methodology steps are similar between studies. Often the metabolomics workflow follows a sequence of steps. The process starts with the Experimental design with the main question of the study, followed by Sample Preparation, then the selection of Data Aquisition, Metabolite extraction and Detection (when aplicable) followed by Data treatment (data preprocessing and data processing) preceded the ending of Data Analysis, as presented in Figure 3. However, several studies choose more than one analytical technique for the metabolomics study, different data analysis methods can be also combined according to the investigation issue to interpret the results and achieve the most accurate analysis [3],[9]. This phenomenon contributes to the lack of standardized guidelines for every step of the workflow in metabolomics studies. Therefore, to ease the data integration and avoid duplication efforts, the Metabolomics Society organized an initiative (the Standards Metabolome Initiative) to formulate a minimum set of reporting standards describing the experiments associated with the statistical and chemometric analysis[16]. 8 2.2. SAMPLE PREPARATION Figure 3: Standard workflow on metabolomics food authencication. Adapted from [9]. 2.2 Sample Preparation Preserving the quality and the diversity of the metabolites at the moment of the acquisition is a requirement in any metabolomics approach to achieve the most accurate analysis. For this reason in any metabolomics pipeline, Sample Preparation is a step with the highest relevance to the whole metabolomics approach. For that, it is necessary to avoid the cell’s natural degradation process, to ensure the physics-chemical features of the metabolites according to their natural properties and location. This method is called the quenching method. Also, in this step, it is imperative to ensure the sample pureness, by removing totally or partially the solvents of the samples. Removing water from aqueous samples through this method can also prevent heat degradation. It is important to note that this method is related to the analytical technique selected[17]. 2.3 Data Acquisition: analytical techniques This section presents the most common analytical techniques used in metabolomics, GC/LC-MS, NMR, IR, UV-vis, and Raman spectroscopies. 9 CHAPTER 2. METABOLOMICS 2.3.1 Mass Spectrometry Mass Spectrometry (MS) is one of the analytical techniques which allows the identification and quantification of metabolites by the determination of the mass charge ratio (m/z) of charges on compounds. This technique can be applied on single forms, as molecules, or on large complex mixtures of molecules or fragments, when the metabolites must be separated. This separation is performed before the MS and is achieved by combining other techniques such as Gas Chromatography (GC), Liquid Chromatography (LC), Capillary Electrophoresis (CE), or Ion Mobility (IM), according to the chemical properties and state of the samples (gas or liquid). After the separation, the different metabolites are then ionized in MS by adding a charge to determine their mass/charge ratio. The results are presented as spectral data, where the x-axis is the mass/charge ratio and the y-axis the intensity/quantity of the different ions in the metabolite [18]. 2.3.2 Nuclear Magnetic Resonance Another analytical technique is Nuclear Magnetic Resonance (NMR) spectroscopy, a non-invasive technique that identifies and quantifies even the metabolites with the same mass. NMR uses the magnetic properties of a selected atomic nucleus to determine the chemical and physical properties and location of those atoms or molecules. The variation of the external magnetic field causes the absorption and reemission of energy by the atomic nuclei to vary. This shift is calculated as the difference between the resonance and the reference substance frequencies, divided by the operating frequency of the spectrometer. Regarding the atom selection, different metabolomics data can be generated based on the atom selected. The most common atom chosen is hydrogen since it is an abundant atom in biological samples (1H-NMR). The results are shown as spectral peaks pattern signal where the x-axis represents the chemical shifts whit the y-axis as the signal intensity, revealing a pattern [19]. In metabolomics approaches usually is selected the one-dimensional NMR spectrum (1D-NMR), and the two-dimensional NMR spectrum (2D-NMR) is chosen only when the metabolites compounds are impossible to identify through the 1D-NMR spectrum because the 2D-NMR allows separating overlapping peaks [18]. 2.3.3 Infrared Spectroscopy Infrared (IR) spectroscopy is another non-invasive analytical technique, which reveals the chemical composition of the metabolites present in the samples. This spectroscopy technique uses the resulting vibration from the absorption of specific frequencies according to the features of the structure of the metabolite. A great advance in IR technique was the introduction of Fourier-Transform spectrometers, because this evolution improved the quality of infrared spectra by reducing the acquisition time, and also allowed the analysis of samples of small size without a sample preparation procedure. According with the analysis length, the results presented as spectrum can be divided in three main regions: the Far-Infrared (FIR), 10 2.4. DATA PREPROCESSING 400cm−1100cm−1; the Mid-Infrared (MIR), 4000cm−1400cm−1; and Near-Infrared (NIR), 13000cm−1 - 4000cm−1[18]. 2.3.4 Ultraviolet-visible Spectroscopy Ultraviolet-visible Spectroscopy (UV-vis) is another non-invasive technique that reveals qualitative information about the metabolites samples by the emission of high energy radiation in the UV(200-400nm) and the visible (400-700nm) range, causing an electromagnetic transition in the molecules at the electron’s level. The UV-vis spectrum results show absorbency values on the y-axis as the wavelength values (nm) on the x-axis [20]. 2.3.5 Raman Spectroscopy Raman spectroscopy is a quantitative or semi-quantitative analytical technique that uses vibrational, rotational, and other low-frequency modes to provide information about the chemical and physical forms of the metabolites. Essentially, this technique uses a tiny portion of light resulting from the interaction between the photon’s light with the matter. The technique uses a single frequency of radiation to irradiate the sample, and it is the radiation scattered from the molecule, one vibrational unit of energy different from the incident beam, which is detected. Most of the scattered light does not change its wavelength in the process, but part of it does, and such scattering is known as Raman scattering. The intensity of the scattered light is plotted against its frequency (cm−1) and the result is a Raman spectrum of the sample. This technique is less used due to problems with sample degradation and fluorescence. Therefore, recent advances in instrument technology have simplified the equipment and reduced the problems substantially. These advances, together with the ability of Raman spectroscopy to examine samples in a wide range of states and with minimal spectrum manipulation needed, have led to rapid growth in the application of the technique [21]. In this dissertation the focus is in 1NMR and FT-IR spectroscopy techniques used to performed metabolomics based approaches to authentication of wine and honey study cases. 2.4 Data Preprocessing Data pre processing is a important step in any metabolomics approach because it addresses the necessary changes to conduct over the data to make the samples analyzable and comparable. Raw data is usually complex and needs to be clean and treated for subsequent analysis. Incomplete records, atypical and inconsistent data can affect the classifier model’s performance and lead to misleading data interpretations [22],[23]. Also it is important to consider that data collection can be an expensive, and time-consuming. Thus to avoid missing data there are exist distinct preprocessing methods. 11 CHAPTER 2. METABOLOMICS In metabolic approaches, the most common pre processing methods include missing values treatment, outliers removal, normalization or scaling, and spectral processing strategies. With these considerations, data pre processing is a crucial step in any metabolomics approach because raw data has the information to answer the food authentication issues or questions. This section covers the principal methods employed in the metabolomics-based approaches to food authentication cases supported by recent studies regarding FT-IR and 1NMR spectroscopies, the techniques used in this dissertation. 2.4.1 Missing values and Outliers A missing value consists of a variable with no assigned value caused by three common mechanisms. The first mechanism is randomly missing values distributed in the data matrix with no dependency on variables. The second mechanism is randomly missing data associated with possible variable dependency to other known variables (X) but not to response variables (Y) [24]. The third mechanism originates nonrandomly missing data in which the missing values have a pattern and dependency within missing values in a variable [25]. There are two methodologies for handling missing data: the removal or replacement. Removing missing values means the elimination of one or multiple variable(s) or a sample(s) that contains the missing data. It is also possible to remove missing values by eliminating variable(s) based on the high percentage of missing values (e.g. more than 5% of missing data per variable(s) or sample)[25]. Replacement is a method in which the missing values are replaced by a numerical value from mean, median, or using other more sophisticated methods (e.g. K-Nearest Neighbor, linear approximation, or other substitute methods) [25]. Outliers are data values that are distant from the other observations and can occur from variations in data measurement or caused by experimental errors. These values are usually excluded from the dataset to avoid interference with the analysis results [26], [23]. 2.4.2 Spectra Treatment IR and 1NMR data processing methods consist of data treatment methods as baseline correction, and noise filtering processes as smoothing. In 1NMR data, peak detection and peak alignment, and data correction (which includes baseline, offset and background corrections) are also applied. Notwithstanding, the spectrum needs to be transformed from FID into frequency spectra using the preprocessing steps of zero-filling, apodization, Fourier transformation, and phase correction before the previously performing the aforesaid preprocessing methods [27]. Relative differences observed between samples are corrected using the methods of Normalization and Scaling. In IR data, mean centering and Standard Normal Variate (SNV) are often used, as well as first 12 2.4. DATA PREPROCESSING and second-order derivatives calculations. In some metabolomics studies with 1NMR data, logarithmic transformations are applied. However, this method does not deal with missing values, requiring a previous treatment of null data [28],[29]. Spectral peaks can suffer dislocations through the x-axis, which can occur due to changes in the chemical environment of the sample like solvents or pH [18]. Spectral alignment procedures are categorized into warping and segmentation. The warping method applied a non-linear transformation to the respective axis. The segmentation method applied a constant shift to all spectral points [18]. Another approach using fast Fourier transformation has been used to perform data alignment [30], [18], [26]. The process of spectral alignment allows to correct the position of matched peaks, that belong to the same metabolic feature(s), in a multiple spectra study [18]. Data correction, which can either be baseline, offset, or background, aims to eliminate the effect of certain signal variations through the spectra and background noise during the sample analysis. Distinct factors can produce these signal variations, as absorption associated with the sample holder and/or solvent used or by either experimental or instrumental variations [18], [11]. These correction methods are very important to avoid false results in data analysis and to not reveal metabolites that are not present in the samples [30]. The smoothing methods filter the spectra noise with significant impact in cases with low noise ratio signals or when the followed data analyses approaches are more sensitive to noise signals. These methods help with the visual interpretation and robustness of the analysis, therefore they require a balance between the noise reduction and the peak retention, especially in the peaks with reduced absorbancy [30]. The Savitzky-Golay filter is one of the most popular smoothing methods. This procedure fits successive sub-sets of adjacent data points with a low-degree polynomial by the method of linear least squares, using a process known as convolution [31]. Another approach is to treat relative differences between samples instead of the absolute values allowing different data to be comparable and consistent. In this approach Normalization, Data mean centering, Scaling, Standard Normal Variate (SNV), Multiplicative Scatter Correction (MSC), and Derivative spectroscopy are the most common metabolomics preprocessing methods used. Normalization methods contribute to a correct measurement of the features in the metabolomics analysis [11], [26], [18], [32]. The subtraction of the row mean to each data value and dividing by the standard deviation is the most common normalization method utilized in 1NMR data [11], [26], [18], [32]. MSC method is a relatively simple processing step that attempts to account for additive and/or multiplicative effects in spectral data. It does so by estimating light scattering or change in path length for each sample relative to that of an ideal sample. MSC is probably the most widely used pre-processing technique for NIR [33]. SNV method, a weighted normalization is performed (not all points contribute to the normalization equally). Therefore, the values are subtracted from their mean, and then the result is divided by the standard deviation. Spectra treated in this manner have always zero mean value and a variance equal to 13 CHAPTER 2. METABOLOMICS one and are thus independent of original absorbance values. SNV is the second most applied method for scatter correction of NIR data. [33]. 2.4.3 Mean-centering Data mean-centering is a column-wise transformation and consists of subtracting the mean spectrum from each sample, resulting in all columns with a mean of zero. This method accentuates the differences and not the similarities of the data samples, so it requires the combination with scaling methods [29], [26],[32]. 2.4.4 Scaling methods Scaling methods divide each variable by a factor, the scaling factor, which is different for each variable. Scaling methods allows samples with variables measured on different scales comparable [29]. Scaling has two subclasses considering the type of measure. Scaling data with dispersion such as the standard deviation as a scaling factor, this measure includes Auto-scaling, Range scaling, Pareto scaling, and Vast scaling. Scaling data using size as a measure, like the mean, includes the Level scaling. • Auto-scaling: Compare metabolites based on correlations. All metabolites become equally important. May inflate the measurement errors. • Range scaling: Compare metabolites relative to the biological response range. Both Auto and Range scaling methods treat metabolites with equal importance but also give a gain effect in measurement errors. Range scaling also inflames outliers values. To deal with this inflation standardization methods are applied as SNV. • Pareto scaling: Reduce the relative importance of large values, but keep data structure partially intact. This method keeps closer to the original measurement however is sensitive to large fold changes. • Vast scaling: Focus on the metabolites with minor fluctuations. Aim for robustness when using prior group knowledge. Therefore is sensitive to higher variations without group structure. • Level scaling: Focus on relative response. Indicate for identification of biomarkers. As auto and range, Level scaling can increase the measurement errors [29], [30]. 2.4.5 First and second-order derivatives First and second-order derivatives are the derivative spectroscopy most commonly calculated for metabolomics analyses. These methods allow to identify significant differences in the derivative mode present in spectra 14 2.5. DATA ANALYSIS samples that are very similar. The separation of overlapping peaks is achieved by differentiation of a zeroorder spectrum and obtaining consecutive derivative spectra. These methods increase the discrimination without the separation of the analytes. The first-order derivative consists of the rate of change of absorbance concerning wavelength. Using the first-order derivative the baseline shifts are eliminated because the resulted constant absorbance offset is zero. The second-order derivative has a very characteristic feature consisting of a negative band with a minimum at the same wavelength as the maximum on the original spectrum [34]. 2.4.6 Data transformation Data transformations, such as logarithmic and power, are non-linear conversions methods usually applied to correct data, as in the case of heteroscedasticity, which is the existence of absolute noise that increases with the rising of signal intensity. Despite both methods correcting the heteroscedasticity, logarithmic transformation does not deal with null values in the opposite of Power transformation that does not require this previous data treatment [28],[29]. 2.5 Data Analysis After the preprocessing methods, the metabolomics data is ready to be analyzed, interpreted, and compared. Metabolomics data reunite vast information and allow to perform different analyses according to the food authentication issue using chemometrics tools, such as Metabolites identification, Exploratory analyses, Univariate and Multivariate Statistical analyses, Unsupervised and Supervised analyses (machine learning) including Feature importance and interpretation. 2.5.1 Metabolite Annotation and Metabolite Identification The Annotation and the Identification of metabolites are two distinct methods used in metabolomics approaches achieved by using NMR and GC/LC-MS techniques. The Metabolite annotation is the assignment of metabolite candidates to the signal based on matching spectra of a studied sample with reference database or library entries spectra. Metabolite identification entails the comparison of the metabolite’s properties of an authentic standard obtained under identical analytical conditions to those of the identified metabolites. [35]. Here it discusses the metabolites annotation by 1NMR, a method utilized in this dissertation. The annotation of metabolites from 1NMR data is possible by matching the spectra of measured 1NMR peaks against the ones from each reference metabolite, acquired under similar conditions to those of the samples analyzed. Each metabolite has a unique and specific spectrum with spectral shifts or differences resulting from 15 CHAPTER 2. METABOLOMICS dependent upon a subset of the training data [51], [32]. To achieve a efficient and robust model, SVMs requires the optimization of parameters which also allow to avoid overfitting problem [32]. ANNs extract linear combinations of the inputs as derived features, and then model the target as a nonlinear function of these features. These models can predict quantitative or classes according to the type of network developed. Generally, the model structure used three layers of interconnected neurons. The first layer consists of an input layer with the correspondent independent variables with one neuron per column, the second (or more) layer(s) is hidden with k neurons that correspond to the weights which will have to be trained to build the model and the final layer correspondent to the output with the dependent variables with one neuron per column. The connections between different neural units can be enforcing or inhibitory in their effect on the activation state of connected neural units, through a limiting threshold function [52]. To use ANN models it is recommended to conduct several iterations because these are nonlinear stochastic methods meaning that in each model performance different results can be produced. These models should be used with caution however when applying modeling parameters robust models can be generated [52]. Decision trees are recursive data structures including decision rules, which may be inferred from the input data through appropriate algorithms [53]. A set of weak classifiers can be combined by using two major classes ensemble algorithms, to create a strong learner with improved performance. In this context ensemble corresponds to multiple models of the same type trained together. The ensembles help to minimize the effects of noise, bias and variance. In Bagging the training set is sampled with replacement, and each learner can be trained in parallel, while in Boosting a serial process is followed where the training data is sampled with repetition, weighting more heavily the samples the were previously misclassified by the previous models in the sequence Random Forest classifiers utilize the Bagging algorithm, and build a large collection of trees, and then averages their output to improve the classification performance [52]. The Gradient boosting (GB), utilizes a gradient descent based method is used to minimize an arbitrary differentiable loss function, while employing a decision tree as base learner. As the name implies this trees are trained using a boosting algorithm, guided by the negative gradient of the loss function being minimized [54]. Some of these and other machine learning models used in metabolomics data analysis applied to food authentication studies available in the literature are listed on the fifth column of Table 1 in chapter of food authentication. 2.5.6 Model Interpretability Machine learning models are generally considered ‘black boxes’, but the ability to easily interpret the results of a predictive model would improve our understanding of the factors underlying the differences between variables, rather than merely predict them. This information can be useful in metabolomics to improve the authentication process and avoid events of fraud. To assess the feature importance, the Shapley Additive 22 2.5. DATA ANALYSIS exPlanations (SHAP) value estimation method was adopted from the field of game theory. This method can be applied for classification problems, as in this dissertation, by computing the average marginal contribution of a feature value across all possible coalitions. [55], [56], [57],[58]. By computing the Shapley values of the best classification models, the selected features with the higher impact on the model performance and their contributions to the model performance are analyzed. An important note to consider is that to compute Shapley values, KernelShap needs to use a reduced number of features due to the time required to work, and also correlated features should be previously eliminated to avoid weighting the features wrongly, leading to an incorrect interpretation of the model performance. 2.5.7 Tools and databases There are several tools and web tools which allow the processing and analysis of metabolomics data. Specifically, webtools have the advantage to help on the comprehension and extraction of knowledge from metabolomics data, without the need of programming skills. To performmetabolite annotation or identification , pathway analysis, enrichment analysis and biomarker identification, the researcher use metabolomics databases such as the presented bellow on 2. Nowadays, with the increasing number of open and web accessible tools, the performance of the analysis and interpretation have also become more precise, robust and also reproducible. Of the web-based tools for metabolomics and spectral data analysis, MetaboAnalyst is the most comprehensive available tool. The specmine R package, developed by the host group is another freely available tool that provides a set of methods for metabolomics data analysis, including data loading in different formats, preprocessing, metabolite annotation, univariate and multivariate data analysis, machine learning, and feature selection. Table 2: Available free tools for metabolomics and spectral data. Tool Type Sepctral Data URL chemometrics R package Chemical data https://CRAN.R-project.org/package=chemometrics ChemoSpec R package NMR, IR and Raman https://CRAN.R-project.org/package=ChemoSpec COLMAR Web tool NMR https://spin.ccic.ohio-state.edu/index.php/colmar hyperSpec R package UV-vis, IR, NMR, MS, Raman, ... https://CRAN.R-project.org/package=hyperSpec MeltDB Web tool GC-MS and LC-MS https://meltdb.cebitec.unibielefeld.de/cgi-bin/login.cgi MetaboAnalyst Web tool NMR, LC-MS and GC-MS http://www.metaboanalyst.ca 23 CHAPTER 2. METABOLOMICS metabolomics R package NMR, GC-MS, LC-MS and MS https://CRAN.R-project.org/package=metabolomics MetaboMiner Java software NMR https://wishart.biology.ualberta.ca/metabominer metaP-Server Web tool NMR, GC-MS, LC-MS and MS http://metap.helmholtzmuenchen.de/metap2/ muma R package NMR, GC-MS, LC-MS and MS https://CRAN.R-project.org/package=muma MVAPACK Toolkit NMR and MS https://bionmr.unl.edu/mvapack.php OpenMS C++ library LC-MS https://www.openms.de specmine R package NMR, GC-MS, LCMS, UV-vis, IR, Raman, ... https://CRAN.R-project.org/package=specmine 2.5.8 specmine and Webspecmine The specmine is an R package created by the host research group for the integrated analysis of metabolomics and spectral data which fills an unmet need of metabolomics. The need for an easy way to process disparate data sources in the open source statistical computing program R. The package handles datasets in different formats (GC/LC-MS, NMR, IR and UV-Vis), helps with the pre-processing data, metabolite annotation, and provides machine learning support [26]. The WebSpecmine (Figure 4) is a website for metabolomics data analysis and mining, which provide simple, interactive front-end and easy-to use tools to cover the main steps of the metabolomics data analysis data, such as NMR, MS, IR, Raman, UV-Vis and concentrations data. It contains modules for data reading and dataset creation, data pre-processing and analysis, as well as metabolite identification, by using implemented functions from the R package specmine. The webspecimine is easy to use, the user does not need to know any kind of programming language, and the results are presented using graphical visualization, so that they can be easily interpretable. Finally, the platform allows the download of any format text, table or graphical. 24 2.5. DATA ANALYSIS Figure 4: Graphical interface of webspecime front-end. 2.5.9 Worflow and Data Analysis pipeline The data analysis workflow in this dissertation encompassed the data preprocessing, statistical and machine learning modeling, and model interpretability. In Figure 5 is presented the general workflow used in the three works presented in this dissertation. Figure 5: General workflow used for the data analysis pipeline. Two distinct pipelines were employed: The first pipeline utilized the specmine R package, developed by the host research group, for spectral data treatment of untargeted metabolomics data of wine samples 25 CHAPTER 2. METABOLOMICS from 1H-NMR and amino acid concentration present in Chapter 4. This approach also included exploratory, statistical, unsupervised, supervised, feature importance and metabolite annotation analysis. The second pipeline was used for the data analysis of untargeted metabolomics data from wine musts samples using FT-IR (Chapter 5) and honey samples using 1H-NMR (Chapter 6), employing freely available Python libraries. Notably, the spectral signal preprocessing of honey samples using 1H-NMR was performed using the specmine R package. The data analysis in these pipelines included exploratory, statistical, unsupervised, supervised, and model interpretability analyses, performed using Python scripts and libraries such as Pandas, Numpy, ResearchPy, Scipy stats, Scikit-learn and SHAP. The comprehensive utilization of these pipelines in the present study filled a critical need in the field of metabolomics. More details are presented in methods section in each Chapter of the presented works. In Figure 6 a schematic link of the worflow steps and the respective Python libraries is shown. Figure 6: Used Worflow on this work and the correspondent Python libraries. 26 3 Metabolomics for authentication in wine and honey food products 3.1 Food Authentication and Traceability The concerns with food authenticity started after the end of World War I, with the term “food security” associated with the traumatic experience lived by European countries who suffered from control and domination of food supplies [2]. In the present days, the Agricultural industry in the European Union (EU) is the second largest industrial sector, and with the globalization of the food market, the demand for the acquisition of certain food products despite the season, location, and availability leads to an increase of the production and distribution guided according with restricted regulation. Economic concerns, and more importantly, safety issues are the main reasons for the regulation of the food systems applied in the EU, and integrated within those guidelines is the directive of the consumer’s rights to receive truthful information about the food they buy (EU regulation No. 178/2002). The main aim of these guidelines is the prevention of fraudulent or deceptive practices, the adulteration of food, and the exercise of other practices which may deceive the consumers [59]. Currently, in Europe, there is also specific legislation for food products according to geographic origin - one of the main issues of authentication. This regulation protects high-quality food products with specific features derived from special geographical conditions or with traditional methods of production [60]. This special certification tries to ensure the safety of the food production system, and the commercial valuation for consumers and producers; aiming to reduce the cases of counterfeit food production and commercialization. To better understand the present work, it is important to present a brief explanation of some concepts and definitions. Food authenticity is an ancient field that tries to validate the veracity of the description of food products and their ingredients, according to the label description [61]. The food product needs to follow the regulation and the standard procedures applied to the food system (production, distribution, and commercialization), adopted by the country to ensure quality and safety. These concerns are integrated with the 27 CHAPTER 3. METABOLOMICS FOR AUTHENTICATION IN WINE AND HONEY FOOD PRODUCTS Food authentication process [62]. This process may differ according to the regulation adopted by the country, which may lead to a panoply of analysis tests to apply. The aforementioned process may also include distinct research topics, such as food origin or traceability, food control, or adulteration. Thus, having these concerns in mind, there is a need to select the proper analytical test and the best chemometric data analysis approach. White the increase of tests also increases the number of studies, in Figure 7 is showing the countries with the most publications related to the authentication process. Figure 7: Top 10 countries on food authentication publications. Adapted from [10]. 3.1.1 Food Control Food Control is one of the most important topics within the food analysis field, and it has to be applied in most stages of production to ensure the quality and the safety of the food products. These concerns encompass the verification of the quality and quantity of certain ingredients, and other elements connected with the legislation applied to the last stage of the production before reaching the final consumers. Notwithstanding, to ensure the safety of the product the analytical exams are performed with the purpose to assess: (1) the presence of selected compounds in foods that may be present below certain thresholds (maximum residue limits - MRL); (2) the detection of microbial-related spoilage; (3) the determination of allergens; (4) the detection of environmental contaminants, as well as banned external compounds; and (5) the assessment of the occurrence of natural toxins [3]. To certify the quality of the product – which is the verification of quality and quantity of certain ingredients and other elements, - analytical tests are conducted to check the presence of food adulterants. Adulteration is described as the deliberate alteration of the product by adding or removing ingredients from the food product, aiming to get an increase in profit. According to Food Standards Agency (FSA) (United Kingdom Agency), there are two main types of adulteration: (1) the transaction of food which is unfit and potentially harmful; and (2) the deliberate misdescription of food [63] [64]. Fortunately, most of the food 28 3.1. FOOD AUTHENTICATION AND TRACEABILITY adulteration cases do not have the intention to harm the consumers wellbeing [65] [66], but have as main consequence to reduction of consumer confidence leading to loss of product value [66] [10]. Adulteration is not a new practice, in 1820, German chemist Fredrick Accum published his book “A Treatise on Adulterations of Food and Culinary Poisons” in which he described methods of adulteration and analytical techniques to detect them, applied to tea, coffee, bread, beer, and pepper. He also highlighted the quality of the water used for food and drug production. His work still serves as a key element in the quality and safety schemes, at present days [67] [68][10]. 3.1.2 Food Origin and Traceability In Europe, the main concerns related to food authenticity are food origin and traceability. Food origin is concerned with the geographical place or region where the food or its ingredients were produced. This also includes the floral species, the traditional methods of production, and the animal diet and/or activity [60]. To protects the food product according to their geographical origin, EU is applied specific legislation for these authentication issues: (1) the protected designation of origin (PDO); (2) the protected geographical indication (PGI) labels for “mountain products” and “island farming products” (EU regulation No.1151/2012); and (3) “traditional specialty guaranteed” (TSG), a previously introduced term (EU regulation No. 509/2006) [10]. This special certification has the purpose to protect the reputation of the regional foods and promoting good practices in rural and agricultural activity. Those practices help the producers to obtain the best prices for their authentic products, reducing the unfair and misleading competition from non-genuine products, usually with inferior quality. With this regulation, there is a reinforcement on the consumer perception about the special quality of the product [69]. Another important definition is the term traceability, meaning the path covered from its first production, the geographical origin, to the final consumer. In developed countries, traceability has the legal obligation to follow the regulation (EU regulation No.178/2002) directed to ensure the quality management schemes and standards (i.e., HACCP, ISO 22000:2005, FSSC 22000, EurepGAP, BRC, etc.). The large applications of the Food Authentication process with their economical profit managed to gain exponential attention in the field, especially for producers, manufacturers, sellers, consumers, agencies, and governments. The food authentication processes managed to lead the attention of producers, manufacturers, sellers, consumers, agencies, and governments due to the large sums of money involved at the industrial level. In the last decades, there has been an increasing interest in scientific research supported by the coordination actions of the European Union, including “Food Integrity”, “MoniQa” and “TRACE” within HORIZON 2020. The process of food authenticity reveals the unique features or the several attributes of food products in evaluation. These specific features add high value to the place and/or the production method, recognized by proper certifications granted by the aforementioned labels PDO, PGI and TSG as presented in Figure 29 CHAPTER 3. METABOLOMICS FOR AUTHENTICATION IN WINE AND HONEY FOOD PRODUCTS 8 and nowadays there are a growing number of products with the authenticity certification as shown in Figure 9. Figure 8: Official GIs logos of the European Union (PDO, PGI and TSG) according to Commission Implementing Regulation (EU) No. 668/2014. Adapted from [10]. Figure 9: Collective number of PDO, PGI and TSG products per EU country. Adapted from [10]. The creation and implementation of specific laws, regulations, and guidelines regarding food authentication in conjunction with the exponential growth of the food industry to profits, did not prevent the occurrence of events concerning food contamination, adulteration, and fraud. In the last years, the horse meat scandal, the melamine adulterated baby milk formula, and more recently the eggs contaminated with pesticides, lead to the alertness of the consumers concerning the quality and authenticity of the food [1]. These events are forcing governments, food agencies, and producers to inquire about new authentication approaches to apply to the food system to deal with food crime [1]. Due to these types of criminal events, the fields of authentication and traceability are under constant pressure regarding the continuous improvement and development of robust, efficient, sensitive, and costeffective methodologies [3]. 30 3.1. FOOD AUTHENTICATION AND TRACEABILITY In the last decades, there has been continuous work on developing, improving, and standardizing methods and practices to assure food authenticity. Nowadays, new studies have proved the potential of metabolomics based approaches in the field of food science and nutrition [3]. Some of these examples are presented in Table 3. The combination of analytical techniques in conjunction with bioinformatic tools for data analysis based on Metabolomics gives the capacity to face some of the important challenges in the field [3]. As presented in the previous chapter, metabolomics, is one of the main areas in the field of -omics techniques, and consists of the exhaustive study of the whole small metabolite composition of a particular system or organism. In practice, food metabolomics aims to analyze the chemical composition of the food metabolome, meaning some quite diverse compounds. For this, several studies choose more than one analytical approach and combine the interpretation of several results. This is the main reason for the lack of standard guidelines described for every step of the workflow in metabolomics studies as mentioned and explained in section 2.1 of chapter metabolomics. Despite, there is no standard workflow due to the non-applicability of the aforementioned steps in all studies [9], in the food authentication process metabolomics follow a common number of steps where the analytical techniques and chemometric approach are selected according to the purpose of the study, as presented in Figure 3 (Chapter 2). [9]. Most food authenticity studies are based on Metabolomics approaches and are applied to food products such as alcohol drinks (beer, wine, spirits), honey, olive oil, milk and cheese products, meat, spices, and juices [70]. The studies reviewed in this work have focused on particular techniques, confirming the proof-of-concept, and also presented robust statistical results based on a large number of samples, contributing to creating predictive models. These studies are presented in this subsection along with the analytical techniques, data analysis tasks, the authentication purpose, the food product tested, and their references in Table3. 31 CHAPTER 3. METABOLOMICS FOR AUTHENTICATION IN WINE AND HONEY FOOD PRODUCTS 3.3.3 Classification for Floral and Geographical Origin Another recent study used H-NMR metabolomics approaches and chemometric tools to determine the floral origin of Finland honey and also their geographical origin, presented by Kortesniemiand their colleagues in 2016. In this investigation, key markers of authentication for different honey types were identified (buckwheat, clover, dandelion, heather, Himalayan balsam, honeydew, linden, lingonberry, multiflora) and the fingerprints of some samples with full discrimination were accomplished with unsupervised (PCA) and supervised (OPLS-DA) models. In the PCA, melezitose, glucose, and fructose explained 83% of the variation, while the origin-specific non-saccharide metabolites were the discriminatory compounds in OPLSDA loadings [127]. To classify floral and geographical origin and to identify fraudulent honey productions, there are other studies with relevant results using metabolomics data approaches combining chemometric tools for data analysis where the data came from other analytical techniques. Gok in 2015, differentiated and classified 120 honey Turkish types by using FTIR spectroscopy techniques combine with unsupervised models (PCA) and hierarchical cluster analysis [128]. Woodcock in 2007 used NIR spectroscopy to evaluate the potential of the technique to confirm the geographical origin of several samples of honey (Irish, Mexican, Irish, Argentinean, Czech, Hungarian, and Spanish honey types) [129]. In Table 3.4 is presented a selection of the most relevant studies in wine authentication metabolomicsbased approaches using NMR, MS, and IR techniques followed by data analysis. Table 5: Metabolomics Approaches for Honey Authentication and Traceability. Purpose of analysis Techniques Preprocessing Data Analysis Reference floral origin 1H-NMR Shift correction normalization, 1st derivative, mean-centering, scaling PCA, HCA, KNN, SIMCA, PLSDA [130] HPLC-UV Integration (s/n <3), scaling to zero mean and unit variance PCA, LDA [131] NIR 1st and 2nd derivative, Savitzky-Golay smoothing bucketing, normalization PCA, MD-DA, BP-ANN [132] 1H-NMR bucketing, normalization PCA, O2PLS-DA, SVM [126] 1H-NMR bucketing normalization, DmodX PCA, PLS-DA [133] FT-IR - PCA [134] 2D H-NMR Normalization PCA, GDA [135] NIR - PLS, BONN [136] SIFT-MS - SIMCA [137] GC-MS/sensory analysis - PCA, CA [138] 1H-NMR - PCA [139] 1H-NMR - GP [140] 1H-NMR - PCA [141] 1H-NMR and - PLS [142] and voltammetric-electronic-tongue (VET) - Bot. and Geo. origin GC-MS Normalization OPLS-DA, OPLS-HCA, SIMCA [143] Floral origin Raman MSC PCA, PLS-DA, SVM [144] NIR 1st and 2nd derivative, SNV SIMCA, UNEQ, POTFUN [145] FT-IR SNV, 1st and 2nd derivative FADA, PLS, PCA [146] NIR None, diferent derivatives, SNV PLS [147] 1H-NMR Mean centering, unit variance PCA, PLS-DA [148] H-NMR None, diferent derivatives, SNV PLS-DA, GP, PLS-GP [149] FT-IR SNV, 1st and 2nd derivative PLS, FDA, SIMCA [150] 38 3.4. CHALLENGES, LIMITATIONS AND POSSIBLE PATHS FOR FURTHER IMPROVEMENT. GC-GC-TOFM - LDA, DPLS, SVM [151] Adulteration Raman AirPLS PLS-DA [152] 1H-NMR None FA, GDA [153] NIR Smoothing, SNV PCA, LS-SVM, SVM, BP-ANN, LDA, KNN [154] FT-IR - PCR, PLS [155] FT-IR Baseline correction, normalization, peak correction, outlier removal PCA, Linear Discriminant Analysis (LDA) [156] 3.4 Challenges, limitations and possible paths for further improvement. The rapid growth of metabolomics-based approaches in food authentication in conjunction with advances of the analytical techniques and data analysis, permits to inform, discriminate and predict several characteristics of food products. These progresses can help to solve problems in safety and quality arenas and at the same time provide important information to the food industry. Despite a distinct range of metabolomics analyses performed a diversified set of food items, most of the studies can be only considered as discriminative (Table 4 and Table ) with very few compounds identified. Nowadays, there are more sensitive, precise and robust techniques which can analyse several features, however like mentioned before the large amount of data can make difficult the analysis. It is necessary not only to implement new data analysis tools but also to validate and to reproduce and standardize the procedures. The major hurdles in metabolomics approaches concerning profiling and fingerprint in food authentication are: the reduced number of public databases with the information concerning the distinct food metabolites composition, and the respective steps of data analysis. Another important challenge is the lack of data regarding metabolites produced based on external environmental factors, such as the weather changes and the different seasons conditions, which may affect the selection of metabolites used as biomarkers for the food authentication. 3.5 Food Products studied in Data Analyses and their Authentications issues In the context of this dissertation, chemometrics tools such as preprocessing methods, multivariate statistical analyses (MSA), and several machine learning (ML) approaches were utilized on metabolomics data from two high-value food products (Wine, Wine musts and Honey samples), and with PDO importance. These tools aided in chemical identification, standardization, classification, and prediction of these food product variables for authenticity issues, mainly related to origin classification problems. 39 CHAPTER 3. METABOLOMICS FOR AUTHENTICATION IN WINE AND HONEY FOOD PRODUCTS In the first study the purpose of the data analyses was to classify Cabernet Sauvignon wine samples according to their harvest year, to determine which and how the climatic influence factors influence the 1H-NMR data, and to predict the age of bottle storage of the samples. In the second study the focus was on the variety classification of wine musts samples of four native Portuguese grape varieties and one international red grape variety, all produced in two Portuguese PDO regions. This study also aimed to interpret how the features with higher importance modulate the classification model. The last study classified geographical and flora Brazilian honey samples from eleven regions with different flora origins. Model interpretability for the model’s whit higher performance and the relation between floral and geographical origins were explored. The main studies in literature related to these analyses are presented and discussed in each of the presented study in their correspondent chatper. 40 4 Harvest, Age of botle and Climatic factors prediction using statistical and machine learning approaches on metabolomics data of Brazilian Cabernet Sauvignon wines Abstract Wine possesses a complex matrix resulting from the combination of several factors such as grape variety(ies) and terroir, among others. This study contributes to understand the influences of the age of bottle storage and of climatic factors on the wines’ chemical composition. A multidimensional metabolomics data set of 1H-NMR spectra and metabolites concentrations from 30 Brazilian Cabernet Sauvignon wine samples, with distinct harvest years, were collected and analyzed. Metabolomics based-approaches were combined with quantitative statistical and machine learning tools. Qualitative analysis revealed that wines with different harvest years presented similar composition of metabolites. Univariate ANOVA shows that 30 relevant 1H-NMR peaks are distributed in three regions, e.g., phenolic acids, sugars, and amino acids and 5 amino acids allowed to discriminate wine samples according to their harvest years. Supervised models were used to predict the harvest years from 1H-NMR data, with an accuracy of 92%. Using linear regression models, the data allowed to discriminate 1H-NMR resonances that were most influenced by the age of bottle storage and the climatic factors, being the temperature the factor with most influence on this case study. 4.1 Introduction Winemaking has a high cost and is a complex process based on the combination of grape varieties, terroir effects, fermentation, harvest timing/ methods, storage conditions and age, as well as many other conditions capable of influencing the wine quality and authenticity [1], [157], [158]. Complementary to the challenges of the winemaking process and wine metabolome analysis, the map of wine production has expanded to regions, such as Argentine, Brazil, Chile, China, United States, South Africa, New Zealand 41 CHAPTER 4. HARVEST, AGE OF BOTLE AND CLIMATIC FACTORS PREDICTION USING STATISTICAL AND MACHINE LEARNING APPROACHES ON METABOLOMICS DATA OF BRAZILIAN CABERNET SAUVIGNON WINES and others, considered as the New Worlds for wine production [159]. Cabernet Sauvignon (CS) is one of the best-known red wine varieties and it is produced in many regions of the world. This cultivar is known for having a long life-cycle (around 214 days), requiring higher temperatures in all phenological stages. Following climatic conditions similar to the ones in Europe, the produced wines have a good quality [93]. Based on different climatic conditions, such as high insulation, humidity, rain, and temperature, maturation parameters can be changed, influencing the times of harvest and storage and ultimately the wine quality [93]. Complementary to the control of the wines’ quality and their authentication, studies focusing on the factors influencing wine’s chemical composition have been carried out, following the global climatic changes in the wine production regions [157]. Due to the complexity of winemaking, several chemical and omics approaches have been successfully applied to the wine science to better understand the mutual wine composition interferences [85], [157], [160]–[163]. Regardless of these advances on wine sciences, the influence of each individual (a)biotic factor in the wine chemical profiles remains, in most cases, unclear and not well established [85], [157], [160], [161]. The combination of different metabolomics tools, such as Nuclear Magnetic Resonance (NMR) and other spectral techniques, leads to robust, efficient, sensitive, and cost-effective analytical methodologies, resulting in large amounts of data on wine metabolomics profiles. 1H-NMR has the sensitivity to identify chemical shifts in each metabolite per sample, as a result of structural and environmental changes. Even small environmental changes can lead to a significant effect on the spectrum. This means that this technique has the capability to distinguish two similar molecules with identical spectra [164]. The resulting data provide an attractive resource for studying the coordinated behavior of wine chemical profiles. However, to cover the needs of data processing and data mining it is necessary to combine metabolomics based approaches with powerful computational tools, such as multivariate data analysis or machine learning tools. These techniques have shown the potential to face challenges in food science, and, more specifically, in wine science. Statistical multivariate analysis and machine learning employ several methods that allow simultaneously the study of many metabolomics features to identify relations between them, which may increase the quality of the data evaluation [18]. In our study, an analysis was performed considering the influence of combined factors such as the harvest years, the age of bottle storage, and a selected set of climatic factors on metabolomics datasets of CS wines. For that, a comprehensive analysis was conducted using bioinformatics tools, including multivariate statistical analysis and machine learning, over untargeted metabolomics data from 1H-NMR spectra and targeted metabolomics data from Gas Chromatography (GC) (amino acids) for metabolite quantification. The wine samples were produced in the same location in southern Brazil, in distinct years, with similar climatic conditions, and different ages regarding the bottle storage. Based on the combination of the metadata available and metabolomics datasets, the study intends to discriminate and predict wine samples with distinct harvest years and also to identify the influence of the age of bottle storage and of a set of selected climatic factors. 42 4.2. RESULTS AND DISCUSSION 4.2 Results and Discussion This study contemplates the analysis of a multidimensional metabolomics data set of 30 Brazilian CS wine samples with distinct harvest years (i.e., 1992, 1994, 1996, 1999, 2002, and 2005) and also distinct ages of storage (i.e., 1, 4, 7, 10, 12, and 14 years). Firstly, an untargeted and qualitative metabolomics analysis was applied to the 1H-NMR spectra data set, followed by a target and quantitative analysis of 14 specific metabolites by GC-FID (amino acids). As mentioned in the methods section below, each group of wine samples was produced in distinct years, in the same geographic origin and kept under similar conditions of light, humidity, and temperature over the time of storage. Metabolomics based-approaches were used in all samples at the same time and the resulting data were combined with univariate and multivariate statistical analyzes, and machine learning tools to answer to the aims of the study. 4.2.1 Metabolites annotation from 1H-NMR resonance peaks – qualitative and untargeted metabolomics analysis To understand the effect of the harvest year on the wine’s chemical composition two qualitative analyses were carried out: one using the 1H-NMR data from all wine samples (considering all the harvest years) and the other focused on each wine sample group (wine samples with the same harvest year). The qualitative analysis of the 1H-NMR profiles of the CS wine study model revealed that samples produced over the years investigated showed the same set of annotated metabolites. A total of 37 metabolites with a high score for annotation were detected in the wine samples, being all of them important for wine composition. The annotation of these metabolites was based on a new algorithm developed by the host group [165], by comparing the wine samples spectra with the spectral references available in 1H-NMR libraries from free databases. The algorithm takes into account a set of parameters for 1H-NMR spectra acquisition, as the type of solvent, the frequency of acquisition, and the organism. Additionally, glutathione, carnitine, 2’-deoxyguanosine monophosphate (dGMP), isoleucine, thiamine monophosphate, deoxycarnitine, and some biogenic amines (e.g., putrescine, cadaverine, spermine spermidine, histamine, tyramine, and tryptamine) were tentatively annotated according the algorithm used, but require further investigations by 2D-NMR experiments and other analytical techniques to prove its occurrence in the wine samples studied. The main metabolites found in this work have been previously annotated in CS wines, such as alcohols, sugars, acids, phenolic compounds, esters, aldehydes, and some amino acids [85], [166], [167]. From the annotated metabolites, several of them are relevant for wine composition and they have been previously reported in wine samples using other analytical techniques. Nonetheless, they are not so commonly found by using 1H-NMR spectroscopy. For instance, sorbitol and myo-inositol, naturally present in grapes, are polyhydric or sugar alcohols which may play a sensory and potential balance role in winemaking [93] particular in dessert wines [168] and are frequently found in wine by other analytical techniques, but only in minor amounts [93]. 43 CHAPTER 4. HARVEST, AGE OF BOTLE AND CLIMATIC FACTORS PREDICTION USING STATISTICAL AND MACHINE LEARNING APPROACHES ON METABOLOMICS DATA OF BRAZILIAN CABERNET SAUVIGNON WINES Several free amino acids were annotated in our study, such as lysine, alanine, leucine, serine, glutamic acid, tyrosine, cysteine, phenylalamine, aspartic acid, threonine, histidine, methionine, glycine, aminobutyric acid, proline, and others. These metabolites are important for the fermentation process allowing the normal activity and grow of yeast, also to originate vitamins and to integrate the flavor and aroma of wine [169]. With these results, our study brings a new highlight to the use of the 1H-NMR technique for the annotation of relevant metabolites with lower concentrations to improve the wine chemical profile characterization, classification and the winemaking process as well. In Table 6, the annotated metabolites with higher score values concerning importance and relevance in wine composition are presented. Bold type indicates amino acids quantified in this study. Table 6: Metabolites and their 1H chemical shifts annotated by 500 MHz 1H-NMR with higher scores for frequency, organism and solvent. Compound Molecular formula δH, ppm (multiplicity, J, Hz, assignment)) Lysine C6H14N2O21.37; 1.40; 1.43; 1.46 (m, γ-CH2); 1.49; 1.52; 1.55; 1.65; 1.68; 1.71 (m, δ-CH2); 1.74; 1.77; 1.80; 1.83; 1.86 (m, β-CH2); 1.90; 2.98; 3.03(t, ϵ-CH2); 3.64; 3.67; 3.71 (t, 6.07, αCH) Serine C3H7NO33.80; 3.83 (dd, 5.51/3.79, α-CH); 3.86; 3.89; 3.92; 3.95; 3.98; 4.01 (m, β-CH2) Methionine C5H11NO2S 2.05; 2.08; 2.11 (m, β-CH2); 2.14 (s, SCH3); 2.17; 2.60; 2.63 (t, 6.2, γ-CH2); 2.66; 3.77 (t, 7.57, α-CH); 3.80; 3.83 Tryptophan C11H12N2O23.25; 3.28 (dd, 15.31/8.08, β-CH); 3.34; 3.43; 3.46 (dd, 15.33/4.77, β0-CH); 3.49; 3.52; 4.01; 4.04 (dd, 8.12/4.02, α-CH); 7.17 (m, C5H-ring); 7.28 (m, C6H-ring); 7.30 (s, C2H-ring) Phenylalanine C9H11NO23.06 (d, 7.88, β0-CH2); 3.22 (m, β0-CH); 3.25; 3.28; 3.95; 3.98 (dd, 7.86/5.31, α-CH); 4.01; 7.28; 7.30 (d, 6.96, C2H, 6H-ring); 7.35 (m. C4H-ring); 7.37 (m, C3H, 5H-ring) Tyrosine C9H11NO33.00 (dd, 14.6, β0-CH); 3.03; 3.06; 3.19 (dd, 7.88, β0-CH); 3.22; 3.89 (dd, 5.16, α-CH); 3.92; 3.95; 7.17 (d, 7,97, C2H-ring) 44 4.2. RESULTS AND DISCUSSION Sucrose C12H22O11 3.43; 3.46 (t, 9,29, G4H); 3.49; 3.52; 3.55 (dd, 9,89/3.87, G2H); 3.57; 3.64 (s, F1H); 3.71; 3.74 (t, 9.55, G3H); 3.77; 3.80; 3.83; 3.86 (dd, 6.46/3.25, F6H); 3.89 (dd, 6.51/3.65, F5H); 3.92; 4.01; 4.04 (t, 8.57, F4H); 4.08; 4.21 (d, 3,88, F3H) Fructose C6H12O63.52; 3.55; 3.57 (m, C1H); 3.60; 3.64; 3.67; 3.71; 3.74; 3.77 (dd, 12.72/1.05, C6H); 3.80; 3.83; 3.86; 3.89 (dd, 9.98/3.45, C4H); 3.92; 3.98; 4.01; 4.04 (m, C5H); 4.08 (m, C5H); 4.12 (m, C3H) Glycerol C3H8O33.52; 3.55 (m, C1H2); 3.57; 3.60; 3.64 (m, C3H2); 3.67; 3.74; 3.77 (t, 6.51, C2H); 3.80 Leucine** C6H13NO20.91; 0.94 (d, 6.1, δ-CH); 0.98 (d, 6.1, δ-CH); 1.62; 1.65; 1.68 (m, γ-CH); 1.71 (m, β-CH2); 1.74; 3.67; 3.71 (t, α-CH); 3.74 Myo-inositol C6H12O63.22 (t, 9,49, C5H); 3.25; 3.28; 3.49; 3.52; 3.55; 3.57 (dd, 9.97/2.88, C1H, C3H); 3.60 (t, 9,98, C4H, C6H); 3.64; 4.01; 4.04 (t, 2,85, C2H); 4.08 Proline C5H9NO21.94; 1.96; 1.99; 2.02 (m, γ-CH); 2.05; 2.08 (m, β-CH); 2.11; 2.29; 2.33; 2.36; 2.35 (m, β0CH); 2.39; 3.28; 3.31; 3.34 (m, δ-CH); 3.37; 3.40 (m, δ0-CH); 3.43; 4.08; 4.10 (m, α-CH) Mannose C6H12O6.34; 3.37 (dd, 9.61/6.36, CH5) 3.40; 3.52; 3.55; 3.57 (t, 9.71, CH4); 3.60; 3.64 (m, CH3); 3.67; 3.71; 3.74; 3.77; 3.80; 3.83; 3.86; 3.89 (dd, 4.15/2.13, CH6); 3.92; 3.95 Maltose C12H22O11 3.22; 3.25; 3.28 (dd, 9.45/8.06, CH3); 3.37; 3.40 (t, 9.57, CH12); 3.43; 3.52; 3.55; 3.57 (m, CH6/CH10); 3.60; 3.64 (m, CH1); 3.67 (m, CH11); 3.71 (m, CH13); 3.74; 3.77 (m, CH2); 3.80; 3.83 (m, CH16); 3.86; 3.89 (dd, 12.30/1.95, CH7); 3.92; 3.95 (m, CH2); 3.98 (m, CH2); 4.01; 4.62; 5.21 (d, 3.81, CH4); 5.25 45 CHAPTER 4. HARVEST, AGE OF BOTLE AND CLIMATIC FACTORS PREDICTION USING STATISTICAL AND MACHINE LEARNING APPROACHES ON METABOLOMICS DATA OF BRAZILIAN CABERNET SAUVIGNON WINES Adenosine C10H13N5O43.80; 3.83; 3.86 (dd, 12.89/3.53, CH15); 3.89; 3.92 (dd, 12.80/2.77, CH15); 3.95; 4.28; 4.31 (q, 3.32, CH5); 4.41 (dd, 5.01/3.46, CH4); 4.46 Uridine C9H12N2O63.77; 3.80 (dd, 12.73/4.41, CH5); 3.83; 3.89; 3.92; 3.95; 4.12; 4.21 (dd, CH3); 4.25; 4.31; 4.34 (dd, 12.02/3.99, CH4); 4.38 Guanosine C10H13N5O53.80; 3.83 (d, 4.07, CH1); 3.86; 3.89; 3.92; 3.95; 4.25 (q, 3.75, CH2); 4.28; 4.41 (dd, 5.11/3.88, CH3); 4.46 Cytidine C9H13N3O5 3.77; 3.80 (dd, 12.90/4.53, CH12); 3.83; 3.89; 3.92 (m, CH12); 3.95; 4.08; 4.12 (m, CH9); 4.17; 4.21; 4.28; 4.31 (t, CH8) 2’-deoxyuridine C9H12N2O52.36; 2.39 (m, CH10); 2.42; 2.44; 3.71; 3.74; 3.77 (dd, CH12); 3.80; 3.83 (dd, CH12); 3.86; 4.01; 4.04 (m, CH8); 4.08; 4.41; 4.46 (dt, 6.45, 4.03, CH9); 4.49 Arabitol C5H12O53.55; 3.57 (dd, 8.29/1.55, CH1); 3.60; 3.64; 3.67 (m, CH6); 3.71; 3.74; 3.77; 3.80; 3.83 (dd, 11.74/2.63, CH3); 3.86; 3.89; 3.92 (m, CH2); 3.95 Sorbitol C6H14O63.57; 3.60; 3.64 (m, CH4); 3.67; 3.71 (d, 286, CH6); 3.74 (m, CH3); 3.77; 3.80 (d, 2.91, CH1); 3.83 (m, CH2/CH5); 3.86 Trehalose C12H26O13 3.40; 3.43 (t, 9.45, CH4/CH10); 3.46; 3.60; 3.64 (dd, 9.91/3.83, CH2/CH8); 3.67; 3.71; 3.74; 3.77 (m, CH2/CH16); 3.80; 3.83; 3.86; 3.89; 5.21 (d, 3.79, CH1/CH7) Salicin C13H18O73.46; 3.49 (m, CH2); 3.52; 3.57; 3.60; 3.64; 3.67; 3.71 (dd, 12.00/5.40, CH11/CH12); 3.74; 3.77; 3.80; 3.89 (dd, 12.01/2.10, CH11/CH12); 3.92; 3.95; 5.09; 5.12; 7.13; 7.17 (dd, CH14); 7.35 (dd, 7.43/1.33, CH20); 7.37 Methionine C5H11NO2S 2.05; 2.08; 2.11 (m, β-CH2); 2.14 (s, SCH3); 2.17; 2.60; 2.63 (t, 6.2, γ-CH2); 2.66; 3.77 (t, 7.57, α-CH); 3.80; 3.836 46 4.2. RESULTS AND DISCUSSION Argininosuccinic acid C10H18N4O61.62; 1.65; 1.68 (m, γ-CH2); 1.71; 1.74; 1.77; 1.96; 2.36; 2.39; 2.48; 2.51(dd, 16.13/9.61, CH14); 2.54; 2.57; 2.63; 2.66; 2.69; 2.78; 2.81 (dd, 16.10/3.43, CH14); 2.84; 2.87; 3.25; 3.28 (t, 6.87, δ-CH2); 3.34; 3.60; 3.64; 3.67; 3.74; 3.77 (t, 6.11,α-CH); 3.80; 4.21; 4.25 (dd, 9.58/3.30, CH4); 4.28; 4.31 2’-deoxyinosine C10H12N4O42.78; 2.81; 2.84 (m, CH3); 2.87; 3.74; 3.77; 3.80; 3.83; 3.86 (m, CH16); 4.17 (m, CH5); 4.21; 4.62 (m, CH4) Valine C5H11NO20.98 (d, 7.0, γ-CH3); 1.04 (d, 7.0, γ0-CH3); 1.07; 2.21; 2.27 (m, β-CH); 2.29; 3.57 (d, αCH); 3.60 Cysteine C3H7NO2S 2.98; 3.00; 3.03; 3.06 (m, β-CH2); 3.92; 3.95; 3.98 (dd, 7.88/3.11, α-CH) Alanine C3H7NO21.46 (d, 7.26, β-CH3); 1.49; 3.74; 3.77 (q, 7.22, α-CH) Choline C5H14NO 3.19 (s, N(CH3)3); 3.49 (dd, 5.79/4.14, CH3); 3.52; 4.01; 4.04 (ddd, CH2); 4.08 Threonine C4H9NO31.30 (d, γ-CH3); 3.55; 3.57 (d, α-CH); 4.21; 4.25 (m, β-CH2); 4.28 Glutamate C5H9NO41.99; 2.02 (m, β-CH2); 2.05; 2.08; 2.11; 2.14 (m, β0-CH2); 2.17; 2.27; 2.29 (m, γ0-CH2); 2.36 (m, γ-CH2); 2.39; 2.42; 3.71 (t, 7.16, αCH); 3.74; 3.77 N-Acetyl-L-aspartic C4H7NO42.02 (s, CH3); 2.48; 2.51 (dd, 15.64/9.98, βCH2); 2.54; 2.66; 2.69 (dd, 15.64/9.98, β0CH2); 2.72; 2.75; 4.38; 4.41 (d, 6.75, N(CH2)) Histidine** C6H9N3O23.19 (dd, 15.53/7.71, β-CH2); 3.22 (dd, 16.08/4.89, β-CH2); 3.25; 3.28; 3.95; 3.98 (dd, 7.71/4.97, α-CH) ; 4.01 Ethanolamine C2H7NO 3.80 (m, OCH2); 3.83 Alpha-aminobutyric C4H9NO20.94; 0.98 (t, 7.56, γ-CH 3); 1.90 (m, β-CH2); 3.71 (dd, α-CH); 3.74 Asparagine C4H8N2O32.81 (m, β0-CH2); 2.84 (m, β-CH2); 2.87; 2.98; 3.95; 3.98; 4.01 (dd, 7.67/4.25, α-CH) Glycine C2H5NO23.52 (s, α-CH) 47 CHAPTER 4. HARVEST, AGE OF BOTLE AND CLIMATIC FACTORS PREDICTION USING STATISTICAL AND MACHINE LEARNING APPROACHES ON METABOLOMICS DATA OF BRAZILIAN CABERNET SAUVIGNON WINES [A] [B] Figure 14: A: PCA 2D scores plot (PC1, PC2) from amino acids concentrations. B: PCA plot of the first five PCs from amino acids concentrations. The wine samples are presented in different colors corresponding to the harvests groups: 1992 (dark-green), 1994 (orange), 1996 (purple), 1999 (pink), 2002 (light-green), and 2005 (yellow). Clustering gathers samples into a defined number of groups based on the used variables by utilizing a distance measure. Unlike discriminant analysis, the number and composition of such groups is often unknown, while techniques such as using elbow graphs can be utilized to estimate a possible number of clusters. Clustering approaches, such as Hierarchical Clustering Analyses (HCA) involve grouping biological samples using all peaks from the spectra or metabolite concentrations. From the results of HCA, it is possible to observe that the samples from 1H-NMR spectra were grouped in the considered harvest years quite well, with a small overlap between the years of 1994 and 1996. The samples from the 1992, 1999, 2002, and 2005 harvests were well discriminated (Figure15 A). When performing the same method over the amino acids contents, the results revealed that only the 2002 harvest was well discriminated, while the remaining ones had significant overlap, as presented in Figure 15 B. [A] [B] Figure 15: A: Dendrogram plot of the hierarchical clustering, with Euclidean distance between wine samples from 1H-NMR data. B: Dendrogram plot of the hierarchical clustering, with Euclidean distance between wine samples from amino acids concentrations. The wine samples are presented in different colors corresponding to the harvests groups: 1992 (black), 1994 (red), 1996 (green), 1999 (blue), 2002 (light-blue), and 2005 (pink). 54 4.2. RESULTS AND DISCUSSION The overlap and the proximity of the samples from the harvests 1994 and 1996 can be related to the influences regarding the proximity of the age of storage (the wines have 12 and 10 years of age of bottle storage, respectively), as their harvest years had similar influences from three climatic factors (temperature, humidity and insulation, present in supplemental material). To better understand those influences, other analyses were made using an univariate linear regression model and are presented in section: Wine age of storage and climatic factors influence. Afterwards, three distinct machine learning (ML) classification models were built, to identify which features were more relevant to discriminate the wine samples according to the harvest year including: Random Forests (RF), Partial Least Squares regression (PLS), and K-Nearest Neighbors algorithm (KNN). The results, using 1H-NMR data as input features, show that these models have a very good performance, reaching 86-92% of accuracy. Utilizing only the data for amino acid concentrations as input for classifying the wines according to their harvest, the best model (RF) had a performance of 66%. Table 9 presents the results from all models. For the 1H-NMR data, the variables, i.e. resonances, with more importance in the classification by the RF model were similar to the ones returned by ANOVA (1.37, 1.40, 1.77, 2.05, 2.14, 2.21, 2.29, 2.39, 2.57, 3.37, 3.4, 3.60, 4.12, 4.28, and 7.13 ppm), while PLS and KNN models considered four other important resonances (2.11, 3.43, 3.95, and 3.98 ppm), still in the same spectral regions. Indeed, as in the previous results from ANOVA, the main features accounting for the classification models were resonances occurring at the amino acids, sugars, and phenolic regions of the 1H-NMR spectra. For the amino acid concentrations as input data, the variables with more importance in the classification by RF were again similar to the ones returning from ANOVA and Tukey test analyses,i.e., Ala, followed by Pro, Lys and Tyr. Only Ile was considered by the model as important features not previous selected by ANOVA and Tukey test analyses, what may point to non-linear relationships or relevant interactions among different variables regarding these amino acids. Table 9: ML models results for the 1H-NMR data and the amino acids concentration with cross-validation accuracy. The first column contains the considered input data as features of the ML model, the second one the applied ML model, the third one the result with the accuracy of the correspondent ML model computed using cross-validation, and the last column the features (1H-NMR resonances - ppm and amino acids) most relevant, with those marked with bold type indicating those shared with ANOVA results. Input data (features) Model Accuracy Features Amino acids concentration RF 66% Alanine,Proline, Isoleucine, Lysine,Tyrosine 1H-NMR spectra PLS 92% 0.85; 1.37; 1.40; 1.77; 2.21; 3.60; 3.86; 3.95; 3.98; 4.01; 4.28; 4.41; 5.21; 5.32; 7.13 1H-NMR spectra RF 90% 1.37; 1.40; 1.77; 2.05; 2.14; 2.21; 2.29; 2.39; 2.57; 3.37; 3.40; 3.60; 4.12; 4.28; 7.13 1H-NMR spectra KNN 86% 0.85; 1.40; 1.77; 2.05; 2.11; 2.14; 2.29; 2.39; 3.40; 3.43; 3.95; 4.01; 4.12; 5.25; 7.13 55 CHAPTER 4. HARVEST, AGE OF BOTLE AND CLIMATIC FACTORS PREDICTION USING STATISTICAL AND MACHINE LEARNING APPROACHES ON METABOLOMICS DATA OF BRAZILIAN CABERNET SAUVIGNON WINES 4.2.3 Wine Age of Bottle Storage and Climatic factors influence In our study, regression analysis is used to identify which peaks and amino acid amounts are correlated with the age of bottle storage and climatic factors individually, considered as continuous (numerical) variables. Linear Regression (LR) models were also used to infer correlation between the independent (age of bottle storage and climatic factors) and dependent variables (1H-NMR resonances and amino acids concentrations). Note that regression analyses by themselves only reveal relationships between a dependent variable and a collection of independent variables in a fixed dataset [183]. With these considerations, other influences not considered in study may affect the metabolomics wine profiles. 4.2.3.1 Assessing the influence of Age of Bottle Storage Bottle storage of wine is a long term process with an important influence in wine preservation, being also used to improve the quality and value of certain wines for a long term aging process, under adequate conditions [93], [169]. Some studies provide important analyses to elucidate the impact of the storage conditions influences on wine composition and quality, such as humidity, temperature, light, and packaging [184]–[188]. Other ones also reported the type of storage influence such the use of bottles or oaks on wine composition and quality [93], [169], [174], [189], [190]. However, the majority of the reported studies only comprise some months or a few years, such as two to five years, and are also focused in a certain type of compounds, such as phenolics and mainly anthocyanins [191], [192]. Due to the high storage cost, only few wine varieties with certain features benefit from the aging process [93], [193]. CS is one variety with good features for the aging process, such as the high tannin content, a certain range of pH values and other ideal features for storage conditions [174], [177], [191], [192], [194]. The pH values of the studied wine samples, before the experiment analyses, were between 3.64 and 3.90, as others also reported similar values of pH for the same wine variety produced in similar climatic conditions [180]. In our study, as mentioned before, the wine samples have distinct years of bottle storage, between 1 to 14 years. To understand the effect of the aging on the wine’s chemical profiles, a supervised approach was applied to the 1H-NMR data set and the metabolite concentrations. In this case, a LR analysis was performed to better understand which peaks and metabolites were influenced directly by the age of bottle storage. The results of the LR model over all 1H-NMR spectra for the wine’s age of bottle storage influence revealed that two 1H-NMR resonance peaks at 3.74 ppm (coefficient: -0.17) and 4.28 ppm (coefficient: -0.16) presented significant p-values (<0.01) and high R2(>55%). Both resonances decrease their intensities levels with the increase of the age of bottle storage, as represented in figures (Figure 16A, 16B). Looking to amino acids concentration influenced by age of bottle storage, the LR models reveal one affected metabolite, Pro (with significant p-value <0.01, R2= 0.32 and coefficient:-65.35 ), whose concentration decreases with the age of bottle storage, as represented in Figure 16(c).These results should be taken with a grain of salt due to the adjusted R2explaining only about one third of the variability in the data. 56 4.2. RESULTS AND DISCUSSION In wines, the concentration of different amino acids depends on the methods used for their determination. Usually, Pro, Ala, Arg, Glu, Ser, and Thr are the amino acids with highest concentration in wine samples [195]. Indeed, because Pro is not consumed in fermentation and maturation processes [196] its concentration is the highest among the remain amino acids and is also higher in some varieties like CS [197], being this confirmed in our data. [A] [B] [C] Figure 16: Plots for the linear regression models showing the influence of age of bottle storage on 1H-NMR resonances and one amino acid with highest adjusted R2.A: 3.74 ppm resonance. B: 4.28 ppm resonance. C: Pro amino acid concentration. 4.2.3.2 Assessing the influence of Climatic conditions Vineyard climatic conditions are some of the environmental factors which modulate the wine chemical composition and, consequently, give a specific fingerprint in the observed metabolomics profile. Furthermore, grape cultivars may respond differently to yearly climatic changes, so it is important to assess which are the climatic factors with highest impact on modulation of the wine metabolomics profile. Usually, CS 57 CHAPTER 4. HARVEST, AGE OF BOTLE AND CLIMATIC FACTORS PREDICTION USING STATISTICAL AND MACHINE LEARNING APPROACHES ON METABOLOMICS DATA OF BRAZILIAN CABERNET SAUVIGNON WINES wines are produced in European regions with maturation average temperature of 17-20ºC [198], [199]. In our study, the CS wines were produced in southern Brazil, under the same winemaking conditions over the harvests, with maturation average temperature of 15-26.4ºC, having microclimatic variations among different years of production. In our study, we considered temperature variations (ºC) (minimum temperature (TempMin), median temperature(TempMean), and maximum temperature(Tempax)), including the minimum mean night temperatures (ºC) (Cool Night Index (Cool Index)), and variations for Precipitation (mm) (PP), Humidity, (%) and Insulation levels ((kWh/m2) per day). Nonetheless, the analysis was conducted using the aforementioned climatic factors to answer the main question of the individual climatic influence. The climatic data used in this study is shown in Table 10. Table 10: Climatic Data by each harvest year Year CoolIndex Ins Hum PP TempMean TempMin Tempax 1992 15.9 206.7 76.5 152.6 20.6 16.1 26 1994 15.7 210.7 75.5 134.0 20.6 16.1 25.5 1996 16.5 225.2 73.3 168.0 20.6 16.1 25.8 1999 18.7 223.9 72.3 83.9 20.5 16.0 26.0 2002 17.0 216.7 80.0 140.1 19.1 15.0 24.7 2005 17.0 240.3 70.3 99.1 20.4 15.6 26.4 The previous results showed differences in the 1H-NMR profiles and amino acids concentrations of wine samples according to the different harvests, so it is important to know which are the climatic factors that influence the wine metabolomics profile, considering the case study. The results of the LR model to study the influence over all 1H-NMR peaks and amino acid concentrations reveal some combinations with a significant p-value (<= 0,0001) and an adjusted R2above 50%. The LR models reveal 12 peaks (Table 11) showing a significant influence by climatic variables. Table 11: Results for linear regression analysis of the climatic influence over 1H-NMR data. The first column contains the considered peaks, the second column presents the influence type, followed by the respective corrected p-value(<= 0,0001), adjusted R2(> 50%), and coefficient of the LR model. Bold type indicates 1H-NMR resonances with highest adjusted R2. Peak (resonance, ppm) Climatic influence p-value Adjusted R2Coefficient 2.05 T Max 9.69E−11 77% 1.63 2.14 T Max 1.29E−06 56% 1.39 2.17 T Max 1.06E−06 56% 1.40 2.17 T mean 4.82E−06 51% 1.32 2.27 T mean 1.07E−07 63% -1.45 2.27 T min 3.04E−06 53% -1.79 58 4.2. RESULTS AND DISCUSSION 2.29 T mean 2.69E−09 71% -1.54 2.29 T min 5.44E−07 58% -1.87 2.39 T Max 1.37E−07 62% 1.47 2.57 T mean 5.42E−07 58% -1.40 3.00 T Max 1.73E−06 55% 1.38 3.60 Cool Index 8.44E−11 78% 0.89 3.60 Precipitation 4.07E−06 52% -0.02 3.89 Humidity 4.06E−07 59% 0.24 3.89 T Max 9.23E−08 63% -1.48 4.28 Insulation 5.65E−07 58% 0.07 5.32 Humidity 1.76E−06 54% 0.24 The results indicate that temperatures (maximum, mean, minimum, and Cool Night Index) are independent variables significant for 10 of the 12 selected peaks, while insulation, humidity, and precipitation influenced only four peaks. A peak at 3.60 ppm is the 1H-NMR resonance most influenced by the Cool Night Index climatic factor(R2of 78%). This peak increases its intensity level with the increase of the minimum mean night temperature as presented in Figure 17, but also decreases its intensity with the increase of precipitation levels (Table11,line:13). The resonance at 2.05 ppm is influenced by maximum temperature (explained variance of 77%), meaning that the peak increases its intensity with the higher maximum temperatures (Figure18A). The signal at 2.29 ppm is influenced by both mean temperature (Figure 18B) and minimum temperature (Table11,line:8). Figures 17 and 18 presents the plots for the linear regression models showing the influence of the climatic factors on the three mentioned 1H-NMR peaks with the highest adjusted R2. To note that, the 1999 harvest year is depicted as the isolated group of peaks with the highest value of Cool Night Index (Figure 17). Analogously, the isolated group of peaks regarding the mean temperature (T mean) concerns the 2002 harvest year (Figure 18B). These results show the influences of the studied climatic factors set on NMR resonances. 59 CHAPTER 4. HARVEST, AGE OF BOTLE AND CLIMATIC FACTORS PREDICTION USING STATISTICAL AND MACHINE LEARNING APPROACHES ON METABOLOMICS DATA OF BRAZILIAN CABERNET SAUVIGNON WINES Figure 17: Plot for the linear regression model showing climatic influences on 3.60 ppm resonance with a highest adjusted R2.The wine samples are presented in different colors corresponding to the harvests groups: 1992 (green), 1994 (yellow), 1996 (brown), 1999 (light-blue), 2002 (red), and 2005 (dark-blue). 60 4.2. RESULTS AND DISCUSSION [A] [B] Figure 18: Plots for the linear regressions models showing climatic influences on 1H-NMR resonances with highest adjusted R2.A: 2.05 ppm resonance. B: 2.29 ppm resonance. The wine samples are presented in different colors corresponding to the harvests groups: 1992 (green), 1994 (yellow), 1996 (brown), 1999 (light-blue), 2002 (red), and 2005 (dark-blue). The majority of the remaining peaks influenced by the climatic factors are related to temperature. The results show that with higher temperatures most peaks have their intensities increased, while the resonances at 2.27, 2.29, 2.57 and 3.89 ppm decrease their intensity with higher temperatures. Considering the other climatic factors, with the increase of insulation, the signal at 4.28 ppm increases its level (this 61 CHAPTER 4. HARVEST, AGE OF BOTLE AND CLIMATIC FACTORS PREDICTION USING STATISTICAL AND MACHINE LEARNING APPROACHES ON METABOLOMICS DATA OF BRAZILIAN CABERNET SAUVIGNON WINES peak is also known for decreasing its intensity with the age of storage), and with the increase of humidity the resonances at 3.89 ppm and 5.32 ppm increase their intensity levels. Based on LR models for CS wines produced in the region under study, temperature seems to be the climatic factor with higher influence on 1H-NMR profiles, specifically the minimum mean night temperature. The 1H-NMR resonances concentrations that did not show a direct linear relationship with climate metadata were not selected. However, this does not mean that other non-linear relations cannot happen between these data. An important note is that the identified resonances influenced by wine age of bottle storage are not the same as the peaks influenced by climatic factors, except for the 4.28 ppm one. This means that the two influential conditions have two distinct groups of peaks. With these results, these two groups of peaks can be used as metabolomics 1H-NMR wine markers to assess the age of storage and climatic influence on wine CS produced in southern Brazil. The main studies reported to understand the climatic influences on red wine quality concern to the influence of one factor or some related environmental factors on some wine metabolites, such as phenolics or other aromatic compounds, making difficult to understand the final metabolomics wine profile as a result from the climatic factors combinations. Another study revealed that light exposure of Merlot berries increased the skin and pulp flavonols, histidine, and valine contents, and reduced the organic acids, GABA, and alanine contents [200]. A larger study used three varieties, Merlot , Cabernet Franc and Cabernet Sauvignon to study the influence of climate factors for a period of six years, reporting the effects of climate as highly significant regarding to vine behavior and berry composition based on the anthocyanin concentration [201]. Our study used only one wine variety and the corresponding amino acids concentrations, thus it is not possible to compare our results with other investigations. To overcome this limitation, the pipeline of our study should be applied on more wine varieties in conjunction with extra metabolomics data. This allows getting a more detailed characterization of wine profiles under the influence of the factors, such as age of storage and climatic factors. Considering the actual climatic changes on the world, vine management practices should be adjusted appropriately to regional growing conditions and grape cultivars. The results allowed us to putatively identify the relevant 1H-NMR peaks in CS wines samples produced in different years under similar climatic influences. 4.2.4 1H-NMR resonances with capability to discriminate wine samples by harvest, age of bottle storage and climatic factors Qualitative analyses of the identified 1H-NMR resonances were performed after the statistical and chemometrics analysis, based on ANOVA, LR models and machine learning results. The main selected 1H-NMR resonances with best predictive capability to distinguish wine samples based on harvest corresponding to the regions of amino acids, phenolic acids, and sugars (Table 12). 62 4.2. RESULTS AND DISCUSSION Table 12: Regions of 1H-NMR resonances (ppm) with predictive capability to distinguish samples selected by ANOVA and ML models. Phenolic acids Sugars Amino acids 5.32, 7.13 3.60, 3.74, 3.89, 4.12, 4.28 1.37, 1.40, 1.77, 2.05, 2.14, 2.17, 2.21, 2.29, 2.39, 2.57, 3.00, 3.37, 3.40 Based on the resonances obtained from the previous results, it is possible to observe that there are 22 peaks utilized in the discrimination of the harvest year. However, there are eight resonance peaks with the capability to discriminate harvests which are not influenced linearly by age of bottle storage, or the studied climatic factors set (in green in Figure 19). The resonance peak of 4.28 ppm has capability to discriminate wine samples by three factors under study (harvest, age of bottle storage, and the insulation climatic factor). There are twelve resonance peaks that allow to predict only the studied climatic factors and not age of bottle storage (in yellow in Figure 19), using a linear model. There is one resonance peak which can predict the age of bottle storage (3.74 ppm) (in blue in Figure 19). 63 5 Multivariate statistical analysis and machine learning approaches to classify musts of Portuguese native grapevine varieties using FT-IR Abstract The identification of musts of grapevine varieties has high interest and social-economic value, particularly for Portuguese native wines. These wines possess unique organoleptic properties, resulting from specific grape cultivars. Globally, Portuguese wines made with native grapes had a growth consumption and economic value in the last years. Fourier Transform Infrared spectroscopy (FT-IR) allows fast and robust spectral recording without sample pre-treatment, but requires spectral signal treatment. Due to the complexity of the generated signals, specific spectral and data pre-processing treatments are needed. Often, these procedures conjugate statistical tools and machine learning approaches to maximize the accuracy of wine must variety classification. This study aims to assess the capability of common high-speed analytical devices present in the wine industry, such as FT-IR, used to control wine must maturation (grape juices before fermentation), to classify Portuguese native musts of grapevine into categories according to their varieties. Another aim is to describe which specific spectral bands allow discrimination of a given variety and how they influence on the model output. Our data analysis was applied on one red international must of grapevine variety (Syrah) and four native Portuguese ones (white varieties: Arinto and Viosinho, red varieties: Aragonez and Touriga Nacional) used for wine production from two PDO provenances in the region of Lisbon. The analytical approach was divided into several stages: exploratory and statistical analysis, variety classification, feature analysis, and a comparative analysis between statistical and Machine Learning approach results. Despite the statistical similitude of the FT-IR spectra among different musts of grapevine varieties (observed through mean spectral differences and ANOVA analysis), Partial Least Squares-Discriminant Analysis (PLS-DA), Logistic Regression (LR) and Support Vector Classifiers (SVC) were the supervised models with highest performances in classifying samples ranging from 0.69 to 0.98 of PR-AUC score for 70 5.1. INTRODUCTION the different varieties. Arinto and Viosinho were the varieties with the best model classification followed by the red variety of Syrah under this metric. Finally, features influencing the model output were identified using Shapley values. These features were linked to the correspondent functional chemical groups. Unique features, spectral bands correlated only with the target variety. Aragonez presents two unique features with positive impact on the model output, the bands at 1127 cm-1 and 1130 cm-1, related to aliphatic compounds. Syrah has one feature related to other three features, with highest impact for this variety whit a positive impact on the model output, the band at 1161 cm-1, associated to strong vibration in acyl groups of esters compounds. Touriga Nacional possesses three unique features, the bands at 999 cm-1, 1015 cm-1 with a positive impact and 1123 cm-1 with a negative impact on the model output. These features are related to monosubstituted alkenes. The white varieties have also a unique feature, the band at 2924 cm-1, related to carbon dioxide compounds with a O=C=O bonds pattern vibrations, with a positive impact to classify Viosinho, and a negative one to classify Arinto. In this context, the identified spectral bands are especially significant in wine authentication, providing useful information for musts of grapevine varieties classification of international and Portuguese native varieties. 5.1 Introduction Wine is one of the most internationally traded agricultural product with a continuous increase in market value. After 1990, wine consumption had an exponential growth worldwide due to the globalization effects [206]. That unprecedented boom led to improvements in the quality and diversity of available wines [206]. Winemaking has a high cost and is a complex process based on the combination of grape varieties, terroir effects, harvest timing, harvest methods, fermentation, storage conditions, and ageing, as well as many other conditions capable of influencing the wine quality and authenticity [1]. Since 2017, according to data published by INE (Portuguese National Statistics Institute), Portuguese wines had exponential growth in consumption. This rise also continued during the global pandemic challenges of 2020. The exportation growth led to an increase in the average price. This event also occurred on PDO (Protected Designations of Origin) Portuguese wines made from native grapevine varieties. Portugal has more than 250 native grapevine varieties and has 31 wine regions with PDO, defined as geographical areas where certain types of grapes are allowed, obeying constraints such as regulated winemaking processes with a maximum of vine yields. These constraints aim to ensure the quality and origin of that beverage, according to the classification of the European Union (OIV, CONGRESS 2011). The Portuguese native grapevine varieties give rise to variety-based wines with unique organoleptic characteristics such as color, aroma, flavor, and alcohol/ acidity balance [207], [208]. It is paramount to ensure the authenticity of the grapevine must of given varieties to avoid not allowed grapevine must blending or other events such as mislabeling. Thus, the capacity to confirm grapevine must variety in a winemaking process has social and economical importance (INE), by managing and preserving the 71 CHAPTER 5. MULTIVARIATE STATISTICAL ANALYSIS AND MACHINE LEARNING APPROACHES TO CLASSIFY MUSTS OF PORTUGUESE NATIVE GRAPEVINE VARIETIES USING FT-IR heritage, the terroirs and the PDO regions. Several methodologies based on isotopic, molecular and spectroscopic techniques, combined with various chemometric methods, have proved to be helpful in winemaking to discriminate distinct wine grape varieties [209], [210], [211], [212], [213]. In addition, the most “easy to use”, robust, economic, and globally applyed analytical technique is Fourier Transform Infrared (FT-IR) spectroscopy, coupled with chemometrics tools, by international organizations, companies, and also in academic research [210], [212]. FT-IR spectroscopy is used in wine research and winemaking, as it provides a simple, fast, costeffective, and non-destructive way to analyze samples without the need for chemicals or preparatory steps [210], [212]. Nonetheless, this analytical technique has also commonly been used in winemaking for grapevine must maturation control, being the most common tool in the wine industry. However, there is still a reduced number of publications and applications regarding the differentiation of samples based on wine must variety, using FT-IR[210], [212]. There is even less research work directed to the classification of Portuguese wines produced with native grape varieties, which achieved increased market value in the last years, especially those which are produced in PDO regions (INE). In this study, machine learning algorithms were applied to a FT-IR dataset coupled with chemometric algorithms aiming to discriminate musts of grapevines according to their varieties. For that, the dataset used came from routine analysis used by industry to control the quality of the wine must before the fermentation process. Thus, these methods permit to reuse the data set for authentication and characterization purposes. These study is the first one to include the red and white Portuguese grape varieties produced in PDO regions of Portugal, besides an international red variety in a metabolomics dataset, combined with statistical, machine learning tools and model interpretability. The application of data analysis tools on metabolomics wine must data provides a comprehensive discussions capable of bringing new highlights to the field of wine production and quality, helping to characterize the Portuguese native wines produced in Portugal, and consecutively contributes as another tool that can be used in their authentication process. The samples of grapevine musts investigated in this work are from Portuguese native grapes, Arinto and Viosinho (white grapes), Aragonez and Touriga Nacional (red grapes), as well as Syrah, an international red grape, from two vineyards of the geographical area of Lisbon district (PDO - Óbidos and, PDO - Alenquer). The grapes used in our study are considered high-quality ones. The major differences between the two white must grapes used in our study are the maturation time, and the bleeding process. Viosinho has a reduced maturation time compared to Arinto, as Viosinho is one of the selected mono-varietal grapes considered to the production and preparation of blending wines, and also considered to add structure and flavor to the fortified Port wines, as well as being used to blend with other grape varieties [208]. The red varieties Aragonez and Touriga Nacional present an average maturation time, but Touriga Nacional grapes produce a wine must with a slightly higher sugar content compared to Aragonez ones [214]. The international Syrah grape variety is well adapted and allowed to be cultivated in the two Portuguese PDO 72 5.2. RESULTS regions considered in this study. It can be used to produce mono-varietal wines and also to blend with other varieties, more commonly with Touriga Nacional [215]. The study approach is divided into the following stages: exploratory and statistical analysis, variety classification feature analysis, and a comparative analysis between statistical and Machine Learning results. The first objective deals with the complex discrimination of varieties based on the visual analysis of the FT-IR spectra and statistical results to identify and understand the chemical discrepancies between the investigated grape varieties. The second goal is to create accurate classification models of wine must variety using five machine learning approaches. These models are combined with evaluation metrics (Precision-Recall Area Under the Curve, Receiver Operating Area Under the Curve, Mathew Correlation Coefficient, and other metrics), feature selection methods, and preprocessing data treatments. The final goal of the work is to provide access to feature importance and enable structural interpretation of model classification. FT-IR coupled with chemometrics data analysis tools may provide an accurate and reliable method to discriminate grape must varieties based on some identified spectroscopic bands, as shown by the results. 5.2 Results This study contemplates the analysis of FT-IR metabolomics data set of 96 Portuguese wine must samples from five wine grape varieties produced in Lisbon region, during two consecutive harvests. Two datasets were organized from the samples based on the wine type of must (red and white). The data were preprocessed as mentioned in Materials and Methods section. Firstly, a exploratory and statistical analysis was conducted, followed by an untargeted and qualitative metabolomics analysis with ML models to classify each grape must variety, followed by analysis of the features with importance on the model classification using Shapley values to better understand the impact and how the features (FT-IR spectroscopic bands) influence on the model performance. As mentioned in the Materials and Methods section, each group of wine must samples was produced in two places of the same PDO region in Portugal and kept under similar conditions of light, humidity, and temperature over the storage time till analysis. Metabolomics based-approaches were used in all samples during the control maturation time and the resulting data were combined with univariate and multivariate statistical techniques and machine learning tools to answer to the aims of the study. 5.2.1 Exploratory analysis Visual observation of the FT-IR spectra of grapevine musts samples, does not allow to discriminate the samples according to their grape must variety (Figure 20A and B). This difficulty is expected due to the type of grapevine musts, the proximity of production region and the biological origin of the grapevine varieties, and winemaking process used on the samples of this study. Looking at the fingerprint region interval [926-1450 cm-1], it was possible to observe that some band 73 CHAPTER 5. MULTIVARIATE STATISTICAL ANALYSIS AND MACHINE LEARNING APPROACHES TO CLASSIFY MUSTS OF PORTUGUESE NATIVE GRAPEVINE VARIETIES USING FT-IR regions with different intensity levels. This fingerprint region is known due to the bending and structural vibrations, which are particularly sensitive to large wavenumber shifts, and permits the identification of specific functional groups. As an example, it was observed spectral differences at the absorption bands around 960 cm-1 and 990 cm-1. These bands are related to aromatic compounds and present an angular deformation at CH sp3 outside the aromatic ring. Usually, absorptions bands between 900-1500 cm-1, with CO bonding are related to organic acids such as the phenolic ones. In the spectral window at 10001100 cm -1 it was also noted some vibrations with different intensity levels associated with CO group present in glucose and fructose molecules. In the 1050-1150 cm-1 spectral window related to alcohols, some CO stretching bands at 1061 cm-1 and 1065 cm-1 bands were also found. Another interval region between 1200-1450 cm-1 showed intense bands differences related to ring vibrations and bending CH sp2 also reported by Murru et al., 2019[212]. Finally, spectroscopic discrepancies were also found at 1800-1880 cm-1, and 2900-3000 cm-1 regions assigned, respectively, to alkenes with streching of double bond CO sp2 and alkanes with streching of CH group at 2950-2960 cm-1 [216], [217]. Besides visual inspection of the FT-IR spectra, in order to explore and identify metabolomic differences concerning the grapevine musts varieties, a statistical analysis was performed. Figure 20: FT-IR spectra plot for wine must samples grape varieties without water region (1450-1800cm-1). A: Plot of FT-IR spectra of red musts of grapevine varieties samples. B: Plot of FT-IR spectra of white musts of grapevine varieties samples. 5.2.2 Statistical analysis To understand if there were statistical differences between the spectral information of each grape variety group an ANOVA followed by the t-Test and Tukey’s HSD test were applied to the FT-IR dataset to determine spectral differences between each grape variety pair with the information regarding the features (i.e., spectral bands) which allow performing linear discrimination of each grape variety pair. This analysis gives an indicative result to estimate the ML algorithms performance to classify each grape must variety. By knowing which variety pair have more features with the ability to discriminate according to the target variety, one will have an idea of which varieties will be easier to achieve a good classification score and the same as the inverse results. With these results one can test and select the 74 5.2. RESULTS best feature selection, evaluation and scoring methods to the models variety classification. By using ML algorithms combined with feature selection, a sub-data set of features is used to classify each variety, but here the selected features can have non-linear relationships between other features. Thus, the other import reason to perform ANOVA analysis is to check if the identified features with linear capability to discriminate varieties are the same as the calculated important features used on model variety classification. In Table 13 the ANOVA results with t-test shows the spectral bands (cm-1) identified with the lowest corrected p values (< 0.01) for each pair variety as well as the functional chemical groups related to the identified spectroscopic bands. 75 CHAPTER 5. MULTIVARIATE STATISTICAL ANALYSIS AND MACHINE LEARNING APPROACHES TO CLASSIFY MUSTS OF PORTUGUESE NATIVE GRAPEVINE VARIETIES USING FT-IR Table 13: ANOVA results. Wavenumber (cm-1) with the lowest corrected p, values (< 0.01). The first column shows the t-test, which consists of the pairs of grapevine musts varieties that were significantly different in terms of means for each feature, the second one the respective sprectroscopic bands, the third column includes the chemical groups related to the wavenumbers identified. Tukey test (pairs of grapevine musts varieties) Wavenumber (cm-1) Chemical Groups Aragonez vs Touriga Nacional 2203, 2207, 2211 Spectral window 2010-2342 cm-1:Wavenumbers related to to sp C-X triple bonds and stretching of carbon-carbon. (3 bands) 1802, 1806 Spectral window 1802-1806 cm-1:related to aromatics compounds with single CH bend. Wavenumbers at 1770-1800 cm-1 and 1735-1810 cm-1 are assigned to conjugated acids and small rings of esters, respectively. (2 bands) 988, 992, 995, 999, 1003, 1007, 1011 Spectral window 988-1011 cm-1:Wavenumbers at 6001000 cm-1 and 1000-1350 cm-1 are related to alkenes sp2 bend and amines with N-C bend, respectively. Wavenumbrs at 10001260 cm-1 are related to alcohols with C-OH bend. And the wavenumbers 988 cm-1, 992 cm-1, 995 cm-1 and 999 cm-1 are related to aromatics compounds with single C-H bend.(7 bands) Aragonez vs Syrah 1123*, 1127*, 1165, 1169, 1173, 1177, 1181, 1184, 1188, 1192, 1196, 1200, 1204, 1208, 1211 Spectral window 1123-1200 cm-1:Wavenumbers at 10001350 cm-1 are related to amines with N-C bend. Wavenumbers at 1000-1260 cm-1 are assignd to alcohols with C-OH bend. Wavenumbers at 1040-1250 cm-1 are related to ethers. Wavenumbers 1123 cm-1 and 1127 cm-1 are related to alphatic compounds. (15 bands, 2 unique features) Syrah vs Touriga Nacional 992, 995, 999, 1003, 1007, 1011, 1015* Spectral window 992-1015 cm-1:Wavenumbers at 6001000 cm-1 are related to alkenes sp2 bend, aromatics compounds with single C-H bend. (7 bands, 1 unique feature) Arinto vs Viosinho 995, 999, 1003 Spectral window 992-1003 cm-1:Wavenumbers at 6001000 cm-1 are related to alkenes sp2 bend, assigned to aromatics compounds with single C-H bend. (3 bands) The results based on the spectral absorbance profiles reveal that Aragonez vs Syrah was the pair of grape must variety with higher discrepant FT-IR bands i.e., 15 signals. Interestingly, all of these bands were detected in the region at 1123-1200 cm-1, which is related to amines with N-C bend, alcohols with C-OH bend, and ethers. The spectral bands at 1123 cm-1 and 1127 cm-1 are related to aliphatic compounds and these features will also be selected as important ones in the red varieties classification models. The variety pair Aragonez vs Touriga Nacional was the second one with more features, i.e., 12 FT-IR bands with linear capability to discriminate the pair variety. The majority of these bands was found to be at 9881011 cm-1 region that has been related to the alkenes sp2 bend and associated to aromatics compounds with single C-H bend. The remain features were detected at 2010-2342 cm-1 spectral window, which is related to stretching bonds of carbon, and also at 1802-1806 cm-1, where the bands have been assigned to aromatics compounds and to conjugated acids. FT-IR bands of the spectral windows at 1735-1810 cm-1 are related to sp2 C-X double bonds related to ester, ketones and anhydrides compounds. The variety pair Syrah vs Touriga Nacional showed the less number of discriminant features, only 7 76 5.2. RESULTS spectral bands at 992-1015 cm-1 region. These features have been related to alkenes sp2 bend and to aromatics compounds with single C-H bend. Considering the variety pair of white grapevine musts varieties, Arinto vs Viosinho, the ANOVA results revealed only 3 features, with spectroscopic bands at 995-1003 cm-1, which have been associated to alkenes sp2 bend and to aromatics compounds with single C-H bend. The interesting results from ANOVA were the identification of some resonances with capability to distinguish the variety pairs Touriga Nacional vs Aragonez (12 bands), Touriga Nacional vs Syrah (7 bands), and Arinto vs Viosinho (3 bands). These results allow us to perceive that some FT-IR spectral windows can linearly discriminate grapevine musts varieties They also reveal that the majority of the FT-IR spectral signals have a common expression in all varieties, as previously noted by the visual inspection of the FT-IR spectroscopic profiles. To overcome these constraint, machine learning models were constructed with the main goal to classify the grapevine musts varieties samples according to the variety. 5.2.3 Supervised Data Analysis: Musts of grapevine variety classification To overcome the challenge to positively predict wine grape must varieties, machine learning models were constructed with the main goal to classify the musts of grapevine samples according to the variety. Firstly, unsupervised models were also executed (PCA, t-SNE, HCA), but they did not show a good cluster separation regarding the grapevine musts varieties investigated. The overall results showed overlap between all varieties, as expected based on the previous statistical analyses and, for this reason, these results are not presented and not discussed in this chapter. The results concerning the non-supervised approach are shown in the supplemental material (Figures S.37, 38, 39, 40, 41, 42) for the interested reader. 5.2.4 Variety classification Based on the previous results of the univariate and the multivariate analyses, it was observed that by using statistical tools there were only a few number of features which allow to linear discriminate grapevine musts varieties. Besides, looking to the results from the unsupervised cluster analyses performed there were a high overlap between the grapevine musts samples from all varieties investigated. To overcome the difficulty of FT-IR complex matrix data, five supervised classification models, GBC, LR, PLS-DA, RF, and SVC were calibrated with hyper-parameters optimization and evaluated, leading to 15 and 5 models for red and white grape varieties, respectively. The performance of all produced models is presented in the supplementary material (in Supplementary Tables) for the interested reader. As previously mentioned, the main objective of this exploratory study is to build classification models with the highest performance and accuracy to classify each wine must variety, and then use interpretability tools to understand how the features with higher importance explain those models. 77 CHAPTER 5. MULTIVARIATE STATISTICAL ANALYSIS AND MACHINE LEARNING APPROACHES TO CLASSIFY MUSTS OF PORTUGUESE NATIVE GRAPEVINE VARIETIES USING FT-IR However, according to the literature, to achieve the most accurate classification and to analyze the importance of features without misleading interpretation, several treatments of data should be careful conducted: • Water region of the FT-IR spectra of grapevine musts samples should be eliminated to avoid noises which influence the models performance. In our data, the water region was eliminated according to the reference of the WineScan manufacturer. • The treatment of FT-IR spectral signals of grapevine musts samples need to be included to avoid noises and other misclassified points [210], [211], [213], [218]–[220]. In our study, after testing several methods combinations, SNV plus standardization achieved the best results on model performance. • To improve the model classification performance, feature selection methods should be used to reduce the problem dimensionality and to facilitate and improve the computation of feature importance [221]. In our study, two filter feature selection methods were used. Two methods allowed to detected and treat constant, quasi-constant, duplicated and correlated features. Only correlated features were detected and treated. The selected method was used to find groups of correlated features and to select from each group the feature with highest variance. As result, the feature selection method returned a sub-dataset of features that are uncorrelated with each other based on the utilized threshold. As mentioned above, the similarity of samples among different grape varieties generated a high number of correlated features. However, in our study if all correlated features were eliminated, only a reduce number of variables would be available and useful information could be lost. Thus, to overcome this constraint, several thresholds of variance were tested among with the default 0.80 (0.70, 0.75, 0.90, and 0.95). The use of Pearson correlation coefficient with the threshold of 0.90 was the method that achieved the best results on the models performance. This method reduce the 433 features to 20 features for the red varieties and 433 to 17 features for the white varieties. As consequence, the computation of features importance and the interpretability of the model output magnitude was also improved in all variety models, by eliminating features with a correlation equal or higher than 90%. The Tables 25 and 26 present in the supplemental material, shows each group of correlated features and the selected feature with the highest variance for red and white grape varieties. • To evaluate and compared the models performance, there are several available approaches and different scoring metrics. In authentication problems as our study, it is necessary to adequate the validation strategy to evaluate the model performance to ensure that the model is able to generalize well for every possible target. Binary classification model was the select approach to predict the target variety and due to the reduced number of samples, repeated nested k-fold cross-validation was the selected method to infer the model performance (details are shown in the Materials and Methods section). To compare each model with one another, we selected different scoring metrics. 78 5.2. RESULTS Since the main aim of the study is to classify a target must variety from any other source. This mean to classify the entire possible positive class target and to take attention to sensitivities. Precision-Recall curve (PR curve) and Receiver Operating Characteristic curve (ROC curve) are well suitable approaches to evaluate binary classification outcomes. However, for class imbalance, as the Aragonez variety, PR curve is more sensitive than ROC one, by changing the curve shape more drastically. Other metrics are suitable for our problem as Matthews correlation coefficient and Balanced Accuracy. These metrics also provide accurate evaluations for imbalance classes in binary classification models. However, these metrics require a defined decision threshold. To evaluate and compare our models results, several metrics (as shown in the supplementary material (in Supplementary Tables) were tested, but focus was putted on the PR-AUC, ROC curve and also at the MCC, with a threshold calculated based on the PR threshold mean. Because this is an exploratory study with contributions towine authentication, the score metric on the PR threshold mean was computed instead of the default threshold, due to the capability of PR curve metric to better test all possible positive class target with attention to the sensitivity of the model. Nevertheless, the selected threshold can be changed according to the sensitivity and cost established or required for the classification model, according to the defined wine authentication regulations. Based on those assumptions, the models from the cross-validation with higher score metric of PR-AUC were chosen to classify each given wine must variety. Table 14 presents the results of classification models with cross-validation for all varieties with the highest performance scores. The results using FT-IR spectral data set as input features show that PLS-DA, LR, and SVM were the models with the highest performance to classify the varieties of our studied grapevine musts samples. For Syrah variety the best model was PLS-DA, reaching 0.97 of PR-AUC. The red Portuguese native varieties Aragonez and Touriga Nacional achieved reduced performances on their classification, as the models with the higher performances were different from the Syrah. For Touriga Nacional, the best model was LR reaching 0.88 of PR-AUC, and for Aragonez the best model was SVC, reaching 0.69 of PR-AUC. The white varieties also achieved higher classification performance with PLS-DA model, reaching 0.98 and and 0.94 of PR-UC for Arinto and Viosinho varieties, respectively. 79 CHAPTER 5. MULTIVARIATE STATISTICAL ANALYSIS AND MACHINE LEARNING APPROACHES TO CLASSIFY MUSTS OF PORTUGUESE NATIVE GRAPEVINE VARIETIES USING FT-IR Figure 23: Top 17 most important features ranked by mean absolute SHAP values. The features are ordered according to their importance of the best model classification for each white wine must variety. The Mean Shap values and their average impact on model output are shown for each feature. A. SHAP values and the average impact of feature’s on model of PLS-DA to classify Arinto. B. SHAP values and the average impact of feature’s on model of PLS-DA to classify Viosinho. Figure 24: The 17 most important features are shown and ordered according to their importance to the best model classification for each white wine must variety. The Shap values and their impact on model output are shown for each feature. Each point represents a sample in the sub-dataset, and the color of the points represents the value of a particular feature for that sample. A. SHAP values and impact of feature’s on model of PLS-DA to classify Arinto. B. SHAP values and impact of feature’s on model of PLS-DA to classify Viosinho. 5.2.7 Comparative analysis between Statistical, Machine Learning approach and Shapley values results Table 15 contains the features computed by shapley values with highest impact on each wine grape must variety, and some of them are unique features, as previous mentioned, these features are spectral bands 86 5.2. RESULTS correlated only with the target variety. It is worth mentioning that some of the features with importance were previously identified by ANOVA. These results allow to understand that statistical tools are useful to discriminate wine grape must varieties, but the use of ML algorithms (with adequate evaluation methods, features treatments, and score metrics), combined with the Shapley values computation, gives more information to discriminate and classify grapevine musts varieties. Table 15: Comparative analysis between Statistical, Machine Learning approach and Shapley values results. Note: Arg: Aragonez; TN: Touriga Nacional; SY: Syrah; LR: Logistic Regression; SVC: Support Vector Classification; PLS-DA: Partial Least Squares-Discriminant Analysis; (-): Negative impact on the model output; (+): Positive impact on the model output. Feature ANOVA Group Varieties ML Variety Classification and Shapley values Correlated Features 999 TN vs SY TN variety: LR model [0.25 impact, (+) impact, 1st position (other varieties have reduced impact, lower ranking order, and (-) impact)] 3 correlated features and TN vs Arg 1015 - unique feature TN vs SY TN variety: LR model [0.14 impact, (+) impact, 2nd position (other varieties have reduced impact, lower ranking order, and (-) impact)] 0 correlated features 1123 - unique feature Arg vs SY TN variety: LR model [0.07 impact, (-) impact, 8th position (other varieties have reduced impact, lower ranking order, and (+) impact)] 0 correlated features 1127 - unique feature Arg vs SY Arg variety: SVC model [SVC model: 0.17 impact, (+) impact, 2nd position(other varieties have reduced impact, lower ranking order, and (-) impact)] 0 correlated features 1130 - unique feature —Arg variety: SVC model [ 0.12 impact, (+) impact, 5th position (other varieties have reduced impact, lower ranking order, and (-) impact)] 0 correlated features 1161 — SY variety: PLS-DA model [0.30 impact, (+) impact, 1st features position (other varieties have reduced impact, lower ranking order, and (-) impact)] 3 correlated features 2924 - unique feature —Arinto variety: PLS-DA model [reduce impact, (-) impact, 15th position (Viosinho variety has the same impact order but a (+) impact)] 0 correlated features In this context, after performing a comparative analysis based on the ANOVA, ML and Shapley values, it can be observed that all varieties, except Syrah, possess unique features with distinct impacts and are related to their chemical compounds. • Aragonez presents two unique features with positive impact on the model output, the bands at 1127 cm-1 and 1130 cm-1, related to alphatic compounds. • Syrah does not present any unique feature. However, it has one feature correlated to other 3 ones, with the highest impact ranking for this variety classification and with a positive impact on the model output, the FT-IR band at 1161 cm-1. This signal has been associated to strong vibration of acyl groups (C-O bond) in esters compounds. • Touriga Nacional possesses two unique features, 1015 cm-1, with a positive impact on the model output, and the signal at 1123 cm-1 with a negative impact on the model output. These features 87 CHAPTER 5. MULTIVARIATE STATISTICAL ANALYSIS AND MACHINE LEARNING APPROACHES TO CLASSIFY MUSTS OF PORTUGUESE NATIVE GRAPEVINE VARIETIES USING FT-IR are related to monosubstituted alkenes with C-O bend of alkoxy groups from esters or alcohols compounds. • Viosinho and Arinto have also a unique feature to distinguish them, the signal at 2924 cm-1, assigned to the O=C=O bonds pattern vibrations. 5.3 Discussion The ANOVA analysis identified Syrah and Aragonez, as the variety pair with more features with the capability to discriminate samples. However, despite the statistical similarities of FT-IR spectra of the varieties investigated, ANOVA analysis identified some features able to discriminate all the red grape variety pairs, and three features able to distinguish the white wines varieties. As expected, the similarity of the FT-IR spectra of the samples used in our study does not allow the discrimination of the grapevine musts by varieties using multivariate statistical analyses. Three unsupervised models, i.e., PCA, HCA, and t-SNE were tested, but did not achieve a good separation of the data related to the varieties. As mentioned above, the results are presented in the supplemental material. Five supervised machine learning models were tested for the variety classification task using PR-AUC as evaluation metric. These different models have been evaluated to understand and then select the model with the best behavior or the most adequate one to classify each variety. In the end of the analyses, the best approaches to classify each variety were the supervised model combined with feature selection and with adequate metric scores. With the tested pipeline, the main goal to classify each grapevine musts variety was achieved, and by using the Shapley values to help interpret the model outputs and feature impacts. There are no other works with model interpretation in the field of wine variety classification. Despite of being an exploratory analysis, the results found are interesting and can be extended to more wine varieties and a greater number of samples. In further works, our pipeline can be reproduce on other metabolomics data set from different techniques to compare the models performance and features with importance. Other previous studies (mentioned bellow) have developed models for grape variety classification based on must, which achieved a performance with accuracy above 80%. However, these did not take attention to the class size and in some works the signal spectra treatment was not mentioned. Although distinct results can be found in the literature, they are not directly comparable to ours as the carried spectral technique, grape varieties, samples type, and the classification pipeline are different. For example Arana et al., 2005 classified two white grape varieties (Viura and Chardonnay) using NIR data, achieving 97.2% accuracy using PLS regression model. In our work, we also applied binary classification to the white wine grape varieties as Arana and collaborators, but we also include three more red grapevine musts varieties. Another difference in our work is that all grape varieties were produced in PDO regions of Portugal. Cozzolino et al. 2012, have used NIR and MIR-ATR spectroscopies to perform a two-case classification of wine variety for Chardonnay and Riesling. They attain 86% accuracy using Linear Discriminant Analysis 88 5.3. DISCUSSION to classify the grape varieties. As in our work, one of the grape variety was produced in two places, which lead to a drop in the model performance. More recently, Murru et al., 2019 performed a multiclass variety classification using Artificial Neural Networks with AUC measure, having a success rate of 83-95%, using also grape varieties data, based on FTIR-ATR. The study closest to our work is a multi-classification model of white grape must varieties based on FT-IR (Roussel et al., 2003) where authors used pre-processing techniques based on Genetic Algorithms combined with a PLS-DA, a multivariate classification achieving 90.4% accuracy with unbalanced classes. In this context, with our data, after performing several models combined with different treatment and selection methods, Syrah, Arinto and Viosinho were the varieties that have the highest discrimination. This fact may be explained due to the biological origin of the varieties and the size of the classes. The computation of Shapley values allowed to understand the impact of the important features identified, but also in a comparative analysis it also allow identifying unique features. All varieties, except for Syrah, possess unique features with distinct impacts on the classification models and are related to their chemical composition. The overall synopsis of the process and the main results are shown in Figure 25. Figure 25: Synopsis of the results with higher performance and the data analysis on FT-IR spectrum. The squares present the classification model and the metric score of best model performance to each wine must variety. 89 CHAPTER 5. MULTIVARIATE STATISTICAL ANALYSIS AND MACHINE LEARNING APPROACHES TO CLASSIFY MUSTS OF PORTUGUESE NATIVE GRAPEVINE VARIETIES USING FT-IR 5.4 Conclusions Our study addressed two main objectives: the variety classification and the feature analysis with classification models interpretation, by reusing the available data avoiding other costs to explore Portuguese grapevine musts varieties. Despite some statistical similarities detected in the FT-IR spectra of the grape varieties, the PLS-DA, LR and SVM models resulted to be good to classify the grapevine musts varieties investigated. The calculation of Shapley values allowed identifying the most influential features and their effect on the models output. We also relate these identified features, i.e., the FT-IR bands, with their correspondent functional chemical groups. These outcomes highlight the importance and contribution of machine learning models, evaluation metrics, and optimization approaches on chemometrics to food analysis and authentication. To infer the influence of maturation, geographic origin, and other factors, further studies should be carried out with more data related to these issues, because our results regarding the grapevine musts variety classification could be affected/influenced by some other factors not herein considered. The identified spectral regions in this work raise further questions, namely if other experimental techniques would detect the same chemical groups or correlated compounds. Follow-up studies should focus on verifying if the same FT-IR bands or spectral regions are identified as relevant features in the wine musts varieties classification based on different techniques to validate our case study. To enrich the classification models of our pipeline more data from other musts of grapevine varieties should be added. Still, the proposed pipeline can serve as a basis for future works in the authentication of grapevine musts varieties, while reducing costs and time by reusing FT-IR data from routine analysis on grapevine musts varieties control. 5.5 Materials and Methods 5.5.1 Samples Ninety six samples of grapevine musts, resulting from the crushing of berries during the control of maturation, from five varieties (16 Arinto, 22 Viosinho, 26 Syrah, 18 Touriga Nacional, and 14 Aragonez), were kindly supplied by the Enology Laboratory of INIAV-Polo de Inovação de Dois Portos, to be used for grapevine musts variety classification using FT-IR spectroscopy data and machine learning models. The samples were produced in two consecutive harvests of 2017 and 2018, in two Portuguese PDO regions around Lisbon disctrict, (Alenquer: 39◦3’ 59.50”N, 9◦6’ 53.74”W, and Óbidos: 39◦11’ 45.69”N, 9◦ 9’ 54.74”W and, 39◦21’ 7.57”N, 9◦1’ 6.26”W). The control of maturation process, conducted by INIAV, used the same conditions over the grapevine musts according to the wine type (red or white), following reference methods. 90 5.5. MATERIALS AND METHODS 5.5.2 Methods 5.5.2.1 FT-IR spectral measurement There was no sample preparation to perform the FT-IR analysis. Before each experiment, the samples were kept at room temperature. A WineScan FT 120 instrument (Foss Electric, Denmark) was used to obtain the FT-IR spectra. The instrument was equipped with a model 5027 autosampler (64 tray, 40 mlcups). A sample volume of 7 ml (standard setting) was pumped through the cuvette (optical path length 37 mm), which is located at the heater unit of the instrument. The temperature of the samples was set to 40◦C. Analysis time took 30 s/sample. Cleaning was automatically programmed to occur every 5 min. The instrument was zeroed before any set of analyses with the zeroing solution (S-6060, Foss Electric). The instrument was standardized before the initial calibration with FT-IR equalizer solution (537811, Foss Electric) and repeated at least once a month. Samples were scanned from 926 to 5012 cm-1 at 4cm-1 intervals. The number of scans generated per sample, the selection of wavenumbers, and the processing of spectra have been fixed by the manufacturer and are not accessible to change by the user. For every wine must sample, two replicate spectra were acquired. A detailed discussion regarding the analytical procedure to carry out the FT-IR measurements can be found in the work reported by Patz and co-workers [220], [213]. 5.5.3 Statistical Analysis and Chemometrics All processing and data analyses were performed by using the Python language, and the libraries pandas, mumpy, ResearchPy, Scipy, statsmodels, and scikit-learn. 5.5.3.1 Spectral Data Processing and Data pre-processing Two datasets were organized from the samples based on the wine type of must (red and white). All samples that had replicates were transformed by computing the respective spectra mean values. The water region [1451-1799cm-1] was removed from all spectra due to noise as advised by the Winescan manufacturer. No missing values and no outliers were detected. Prior to data analyses, it was necessary to used appropriate pre-processing data and spectral data signal treatments to minimize physical effects as presented bellow [33]. SNV plus Standardization was the selected treatment utilized to perform the data analyses. So, multiplicative interferences in baseline shift and curvilinearity were first corrected by SNV [218]. This was followed by z-Score Standardization by removing the mean band values and dividing by the corresponding band standard deviation for each band individually. 91 CHAPTER 5. MULTIVARIATE STATISTICAL ANALYSIS AND MACHINE LEARNING APPROACHES TO CLASSIFY MUSTS OF PORTUGUESE NATIVE GRAPEVINE VARIETIES USING FT-IR 5.5.3.2 Data analysis Exploratory and statistical analysis A preliminary analysis included statistical univariate techniques and unsupervised multivariate statistical models to understand the samples relationships among the different varieties according to the wine type (red or white wine type). To understand if there are statistical differences between different combinations of wine must varieties a parametric test was conducted. ANOVA and t-test, between each varieties pair were carried to identify which specific spectroscopic bands (cm-1) could be useful to discriminate between the investigated classes. The paired differences with p-value < 0.01 were considered to had statistically significant resonances differences. The unsupervised multivariate statistical models applied to the FT-IR dataset in this study were Principal Component Analysis (PCA), Hierarchical Clustering Analysis (HCA) and t-Distributed Stochastic Neighbor Embedding (t-SNE). PCA is a widely used tool for interpreting and visualizing the main feaures of FT-IR data for classifying samples. In thi is study, PCA was used to investigate if there was a clear discrimination between the samples with different grapevine musts varieties in the latent space. HCA used an Euclidean distance between samples and the complete linkage method. This alternative was chosen because it tends to find compact clusters of approximately equal diameters and does not force clusters together due to single elements being close to each other, like in single linkage clustering. In it’s turn, t-SNE is a method which applies a non-linear dimensionality reduction technique to keep the similar data points close in a lower-dimensional space. This method uses a student t-distribution to compute the similarity between two points in lower-dimensional space. Variety classification and features analysis To classify each variety a one-vs-all approach was followed, employing binary classification to separate a target class from any other set of grapevine musts. The main objective was not to identify the grapevine musts present in the data set, but to classify a target variety from any other source. Since this is an exploratory work with potential to add more varieties in further studies, the binary classification is the best approach to apply into this scenario. The following models were used: Random Forest Classifier (RF), Gradient Boosting Classifier (GBC), Logistic Regression (LR), Partial Least Squares-Discriminant Analysis (PLS-DA), and C-Support Vector Classification (SVC). Several performance metrics were tested and compared (Area Under the Precision-Recall Curve (PR-AUC), Area Under the Receiver Operating Characteristic Curve (ROC-AUC), Matthews Correlation Coefficient (MCC), Balanced Accuracy (BA), Precision, Recall, Accuracy), being the PR-AUC the selected metric to evaluate the models performance. Due to the reduced number of samples, repeated nested k-fold cross-validation was utilized to infer the model performance. The outer loop was executed with k = 10, while the inner cross-validation loop was executed with k = 5. The inner cross-validation loop serves to fine-tune the model hyper-parameters, while the outer loop to estimate the model performance. It is important to be aware that all the models 92 5.5. MATERIALS AND METHODS are trained on the inner loop step. The whole process was repeated 10 times by reshuffling the data set. The mean and standard deviation of performance of the best models were computed (concerning the outer loop). Table 16 shows the hyperparameters tested for each model. Table 16: Hyper-parameters used by the distinct models. Model Hyper-parameters Random Forest Classifier N of Estimators: 10, 50, 100; Max N of Features: auto, 2, 5; Max. of depth: none, 2, 5, 10; Class weigth: balanced. SVC C-Support Vector Classification C: 0.1, 1, 10, 100; – Support Vector Machine Classification Gamma: 1, 0.1, 0.01, 0.001, 0.0001; Gamma: scale, auto; Kernel: rbf, linear; Class weigth: balanced Gradient Boosting Classifier N of Estimators: 10, 50, 100; Learning rate: 0.1, 0.01, 0.001, 0.0001; Max N of Features: auto, 2, 5; Max. of depth: none, 2,5,10; Class weigth: balanced Logistic Regression C: 0.00001, 0.001, 0.1, 1, 10, 100; Class weigth: balanced Partial least squares-discriminant analysis N of Components: 1-11. For each training set and test set, the following preprocessing was carried on: each model was calibrated using the treated raw spectra with the SVN method, using the aforementioned preprocessing pipeline. The calibration was performed with explicit feature selection. The explicit feature selection method consisted in dropping correlated features above a given threshold (in this work 90%). First, the correlation matrix of the spectral resonances was computed. Next, the correlation matrix was further processed by forming clusters of bands with correlation values above the threshold. For each cluster, the band with the highest variance is kept and the remaining ones were dropped. Feature analysis To assess the feature importance, the best models from each class were chosen based on the superior performance in the aforementioned cross-validation scheme for each target variety. Afterwards, each band influence on the target classes was inferred using the Shapley values. These values are the average marginal contribution of a feature value across all possible coalitions and can be applied for classification problems, as in our study [56] [58]. For binary classification, SHAP values for classes 0 and 1 are symmetrical. The selected feature contributes to a certain amount towards class 1, at the same time reduces the probability of being class 0 by the same amount. 93 6 Machine learning approaches to classify honeys based on flora and geographic origin from Santa Catarina, Brazil using 1H-NMR Abstract The identification of the flora and geographic origin of honey has high interest and economic value, particularly for the honey with Geographical Indication (GI), which possess unique organoleptic properties, resulting from unique fauna and flora available in native forests and protected areas. Globally, honeys produced in Brazil had an exponential growth in production and consumption since 2020. Nuclear Magnetic Resonance (NMR) allows fast and robust spectral recording without sample pretreatment, but requires spectral signal processing. Due to the complexity of the generated signals, specific data pre-processing methods need to be applied, which are typically followed by a conjugation of statistical tools and machine learning approaches to maximize the accuracy of classification models of honey as to their flora and geographic origin. This exploratory study aims to assess the capability of high-speed analytical devices present in the authentication and quality-control honey industry, such as 1H-NMR, used to classify Santa Catarina (SC) kinds of honey, according to their region of production (from the eleven agroecological regions) and flora group origin. Another aim is to describe which specific features allow for the discrimination of a given region and flora group and how they influence the model output. Our data analysis contemplated 56 honey samples from the eleven agroecological regions of Santa Catarina state, southern Brazil. These samples have been produced in geographic regions where seven flora types are predominant for bee foraging, as follows: “Uva do Japão” ( Hovenia dulcis ), “Silvestre”, Eucalyptus ( Eucalyptus spp ), “Uva mix”, and Other (representing several multifloras). The devised approach encompasses the following stages: (i) exploratory and statistical analysis of the 1H-NMR dataset; (ii) origin classification (region and flora) and feature analysis, and (iii) comparative cross-analysis between the results of geographic and flora classifications. Despite the statistical similitude of the 1H-NMR data from different types of honeys, machine learning 94 6.1. INTRODUCTION models reached good discrimination performances. Logistic Regression (LR), Support Vector Classifiers (SVC), Random Forest (RF), and Partial SquaresDiscriminant Analysis (PLS-DA) were the supervised models with the highest performance ranging from 0.60 to 0.91 for the Area Under the Curve of Receiver Operating Characteristic (ROC-AUC) for the different geographic regions origin. PLS-DA and SVC were the supervised models with the highest performances ranging from 0.56 to 0.72 of ROC-AUC for the different flora flora origins. Finally, features (i.e., 1H-NMR resonances)influencing the model output were identified using Shapley values. These features were linked to the correspondent molecules found in the honey investigated, and a comparative analysis based on the feature importance results for region of production and flora origin was conducted. It was observed that there are shared important features for the sample classification as to their geographic regions and flora groups. These results thus highlight some common features (resonances) with influence either on flora and the correspondent producer region of honey at the same time, which are especially significant in honey authentication, providing useful information for honeys flora and geographic classification of honeys produced in Santa Catarina, Brazil. 6.1 Introduction Honey is a natural food produced by bees from the nectar of flowers (floral honey), but it can also be obtained from living plant exsudates or from excretions of plantsucking insects producing honeydew, such as the “Bracatinga” one, a typical honeydew produced in southern Brazil [222], [223]. Honey is one of the most internationally consumed agricultural food, with a continuous increase in market value, and since 2009, its production and consumption have increased, generating strong interest in the honey trade globally [224], [225]. According to the statistical report by M. Shahbandeh, the average annual production volume of natural honey in the world reached a peaked in 2017, with 1.88 million metric tons and has reduced since then with the registered production of around 1.77 million metric tons in 2020. Despite the challenges of the Covid-19 pandemic in 2020, the interest and the consumption of natural honey has grown and its market value has also increased, reaching around eight billion dollars in that year [226]. Certified honey is a food product resulting from biological flora and fauna areas without any chemical treatments. Brazil has the largest extension of natural forest in the world, with rich native flora, good resources for the development of the bees, and tropical weather in most of its areas, without pesticides and/or other chemical residues [227]. The high biodiversity of flora and fauna of several protected geographical areas offers the conditions to produce honey with high quality and in quantity. With the regulation of standard requirements established by CODEX Alimentaris and the Brazilian government, the Brazilian honeys increased their commercialization around the world and in UE markets. The certifications issued to that food allow controlling its quality, 95 CHAPTER 6. MACHINE LEARNING APPROACHES TO CLASSIFY HONEYS BASED ON FLORA AND GEOGRAPHIC ORIGIN FROM SANTA CATARINA, BRAZIL USING 1H-NMR To overcome the difficulty of classifying honey samples from such a complex dataset matrix, five supervised classification models, GBC, LR, PLS-DA, RF, and SVC were calibrated with hyper-parameters optimization and evaluated, leading to 55 models for region of production and 25 models for floral groups. The performance of all produced models is presented in the supplementary material (Table 28 and Table 29) for the interested reader. As mentioned in the introduction section, the main objective of this exploratory study is to build classification models with the highest performance and accuracy to classify each sample regarding its production region and flora origin, and then use interpretability tools to understand how the features with higher importance explain those models. However, according to the literature, to achieve the most accurate classification and also to analyze the important features without misleading interpretation, several treatments should be carefully conducted: • Solvents region of the 1H-NMR spectra signal of honey samples should be eliminated to avoid noises that influence the models’ performance, in some cases can increase the models’ classification with misleading results. In our data, the D2O region was eliminated according to the literature [130]. • The treatment of the spectra signal of honeys needs to be included to avoid noises and other misclassified points [130]. In our study, after testing several method combinations, baseline correction, missing values treatment, followed by standardization achieved the best results on model performance. • To improve the model classification performance, feature selection methods should be used to reduce the problem of dimensionality and to facilitate and improve the computation of feature importance [221]. In our study, two filter feature selection methods were used. Two methods allowed the detection and treatment of constant, quasi-constant, duplicated, and correlated features. Correlated features were detected and treated by eliminating features with a correlation equal to or higher than 90%. The selected method was used to find groups of correlated features and to select from each group the feature with the highest variance. As result, the feature selection method returns a sub-dataset of features that are uncorrelated with each other based on the utilized threshold. The parameters of these methods are shown in the methods section. The feature selection allows reducing the 36 features (resonances, ppm) to 29 features. Table 19, shows each group of correlated features and the selected feature with the highest variance for the target variable. • An adequate validation strategy to evaluate the model performance is necessary to ensure that the model can generalize well for every possible target. In our study, we selectedthe binary classification model approach to predict the target variable and, due to the reduced number of samples, with leave-one-out cross-validation to infer the model performance (details are shown in the methods section) in the trainning set and a test set to validate the generalization capabilities with unforseen data. To compare each model with another, different scoring metrics were used. The main aim of the study is to correctly classify target honey by Region and Flora origin from any other source. 102 6.2. RESULTS Table 19: Feature with highest variance and the correspondent correlated features group. First column presents the feature with the higher variance correlated to a set of peaks given a predefined threshold. These features are the selected to the new sub-dataset used at all models classification. And these features were used by Shapley method to interpret their impact on each model classification. Second column are the correspondent correlated features to each feature. Note: * are features with high variance and with lower threshold correlation (0.90) with other features. Feature with highs variance(Peak (ppm)) Correlated features group (Peak (ppm)) 1.34* 1.34 2.07* 2.07 3.23* 3.23 3.26* 3.26 3.39* 3.39 3.41 3.41, 3.44 3.47* 3.47 3.50* 3.50 3.53* 3.53 3.56* 3.56 3.58* 3.58 3.69 3.66, 3.69 3.72 3.81, 4.02, 4.12 3.75* 3.75 3.78* 3.78 3.84* 3.84 3.90* 3.90 3.99* 3.99 4.04* 4.04 4.17* 4.17 4.21* 4.21 4.31* 4.31 4.64 4.65, 4.64 4.96 5.38, 4.96 5.12* 5.12 5.24* 5.24 5.31* 5.31 5.41* 5.41 In our case, the Area Under the Curve of Receiver Operating Characteristic (ROC-AUC) is suitable to evaluate binary classification outcomes for imbalanced classes as in our study. This metric is less sensitive than Area Under the Curve of Precision and Recall (PR-AUC). Other metrics suitable for our problem are Matthews correlation coefficient and Balanced Accuracy. These metrics also provide accurate evaluations for imbalanced classes in binary classification models, however, they require a defined decision threshold. To evaluate our models’ results we tested several metrics (as shown in the supplemental material (Table 28 and Table 29). However, because this is an exploratory study with contributions to honeys classification, we focus the scoring metric on the ROC-AUC and also at the MCC, with a threshold calculated based on the ROC threshold mean due to the main importance given to the global importance of the classifications. Nevertheless, the selected threshold can be changed according to the sensitivity and cost established or required for the 103 CHAPTER 6. MACHINE LEARNING APPROACHES TO CLASSIFY HONEYS BASED ON FLORA AND GEOGRAPHIC ORIGIN FROM SANTA CATARINA, BRAZIL USING 1H-NMR classification model, in accordance with the defined honey production regulations. Based on those assumptions, the models with higher score metric of ROC-AUC in the test set were chosen to classify each given Region and Flora group. 6.2.4.1 Geographic Regions classification Table 20 presents the results of classification models for all regions origins with the highest performance scores. Of the eleven agroecological regions of SC state, the Region 1B “Litoral de Florianópolis e Laguna” was the honey production zone with the best model classification performance with the LR reaching 0.91 of ROC-AUC. In Region 3C “Noroeste Catarinense” the highest performance score was achieved with the LR model 0.84 score of ROC-AUC, and Region 5 “Planalto Serrano de São Joaquim” with 0.69 score of ROC-AUC with the SVC model. These two regions were the ones with higher classification model scores followed by the other five regions performances. Region 4A “Campos de Lages” achieved the higher score of classification with a performance of 0.73 with SVC model. Region 3B “Planalto Norte Catarinense” with the RF model achieved 0.70 of ROC-AUC performance classification. Region 4B “Alto Vale do Rio do Peixe e Alto Irani” with PLS-DA achieved a performance classification score of 0.69 ROC-AUC. Region 1A “Litoral Norte, Vale dos Rios Itajai and Tijucas” reaching 0.65 score of ROC-AUC with LR model classification. And Region 2C “Vale do Rio Uruguai” with the PLS-DA model achieved a classification performance 0.60 of ROC-AUC. Table 20: Performance scores for best models classification for all Regions. All models had pre-process methods with the application of Feature Selection (Drop of constant features with 90% of correlation). The MCC evaluation metric followed a threshold based on the calculated ROC Threshold mean. Rows with background color gray are the regions with higher classification models. Note: PR-AUC: Area Under the Precision-Recall Curve; ROC-AUC: Area Under the Receiver Operating Characteristic Curve; MCC: Matthews Correlation Coefficient; BA: Balanced Accuracy; ROC Threshold: Receiver Operating Characteristic. Target Region ML Model ROC-AUC PR-AUC MCC BA ROC Threshold 1A LR 0.65 0.12 0.15 0.52 0.08 1B LR 0.91 0.23 0.30 0.71 0.10 2A RF 0.48 0.12 -0.07 0.45 0.09 2B GBC 0.43 0.06 -0.09 0.43 0.00 2C PLS-DA 0.60 0.26 0.30 0.61 0.33 3A PLS-DA 0.47 0.10 -0.02 0.49 0.15 3B RF 0.70 0.17 0.18 0.65 0.06 3C LR 0.84 0.47 0.20 0.67 0.03 4A SVC 0.73 0.13 0.16 0.64 0.10 4B PLS-DA 0.69 0.12 0.09 0.58 0.11 5SVC 0.82 0.23 0.41 0.78 0.12 104 6.2. RESULTS 6.2.4.2 Flora groups classification Although it is known that SC state has several native plants in the majority of its regions, spread in the Atlantic Rainforest and Araucaria biomes, for instance, creating a unique type of flora, there are some flora species in each region that allow producing honey with some peculiar traits. Here, instead of classifying each honey sample by the geographic region of production, we looked at the type of biological origin. Our data set resulted in honey samples with a vast combination of flora (different types of trees, bushes, and small flowers), but there are some samples where the major flora origin is known and registered by the producers. Because our samples came from controlled honey production system where the geographic localization and the major pollen donor plants are known for each apiary, we reduced the flora information to five groups of flora. Each group is designated with the main composition of existing flora origins. Table 21 presents the results of classification models with validation for all flora groups with the highest performance scores. The “Uva do Japão” was the flora group with the best model classification performance, with the PLS-DA reaching 0.72 of ROC-AUC. There were other flora group with similar performance score for the models’ classification the group Other (monofloral honey samples) with the SVC model with achieving a performance of classification of 0.67 of ROC-AUC. ”Silvetre”as native flora group with the SVC model reached the classification performance of 0.56 ROC-AUC. These three flora groups were the ones with higher classification model scores followed by the other two flora groups with reduced performances under 0.50 of ROC-AUC (in the Test set), not achiving a positive classification for Uva Mix and Eucalyptus flora groups. Table 21: Performance scores for best models classification for all Flora honeys. All models had pre-process methods with the application of Feature Selection (Drop of constant features with 90% of correlation). The MCC evaluation metric followed a threshold based on the calculated ROC Threshold mean. Rows with background color gray are the regions with higher classification models. Note: PR-AUC: Area Under the Precision-Recall Curve; ROC-AUC: Area Under the Receiver Operating Characteristic Curve; MCC: Matthews Correlation Coefficient; BA: Balanced Accuracy; ROC Threshold: Receiver Operating Characteristic. Target Flora ML Model ROC-AUC PR-AUC MCC BA ROC Threshold Uva do Japão PLS-DA 0.72 0.41 0.22 0.64 0.21 Eucalyptus RF 0.44 0.09 -0.05 0.46 0.04 Silvestre SVC 0.56 0.14 0.09 0.57 0.12 Uva Mix SVC 0.47 0.07 0.04 0.54 0.07 Other SVC 0.67 0.70 0.37 0.67 0.57 In our analyses there are some considerations to note that can influence the classification of the models’ performance: • Seasonality, altitude, and climatic changes are some of the influences on food product production like honey. 105 CHAPTER 6. MACHINE LEARNING APPROACHES TO CLASSIFY HONEYS BASED ON FLORA AND GEOGRAPHIC ORIGIN FROM SANTA CATARINA, BRAZIL USING 1H-NMR • Pollution, fertilizers, and human activities like fires also configure critical pressions on honey bees production [233], [234], [235], [236]. • The vast diversity of native plants, animals, and or insects interactions are other influences factors to difficult the honey classification. An example of a beneficial relationship is the “Bracatinga” honeydew. This honeydew is the result of plant-sucking insects which infest the tree species named “Bracatinga” (Mimosa scabrella) . This honey has an IG classification with production in Santa Catarina state since plant specie is more prevalent in highlands in SC. To better understand which of these influence factors have more impact on the model’s Santa Catarina region classification and how these influences can modulate the region classification further studies should be carried out. However, despite the complexity of the data and confounding influences, the results reveal that LR, SVC, RF, and PLS-DA, were the models with better performance to classify each Geographic Regions classification, with ROC-AUC scores reaching 60-91%. PLS-DA, and SVC were the models with better performance to classify each flora group of honey, with good ROC-AUC scores reaching 56-72%. 6.2.5 Features analysis of models classification ML models are generally considered ‘black boxes, but the ability to easily interpret the results of a variety of classification models would improve our understanding of the factors underlying the differences between each production honey region, rather than merely classify them. This information can be used in the geographic honey authentication process to avoid events of fraud. Using the 1H-NMR dataset of the honey samples produced in southern Brazil, itiis worth mentioning there are no studies addressed to the classification models interpretation, despite a large number of works to classify honeys according to the geographic and biological origin [125],[229], [237], [83]. Our work followed this direction, to find out which resonances are the most predictive to classify each production region, and how these features are related to the model performance. Thus, by determining which important features (resonances, ppm) contributed most to a certain region would help to uncover potential biomarkers resonances or to identify groups of chemical compounds to improve honey classification. To assess the feature importance from the 1H-NMR dataset, we used the Shapley Additive exPlanations (SHAP) algorithm, a value estimation method from the field of game theory. This method can be applied for classification problems, as in our study, by computing the average marginal contribution of a feature value across all possible coalitions. [57], [56], [58]. By computing the Shapley values of the best classification models for production region, it was analyzed the selected features with the higher impact on the model performance and how they contribute to each variety classification. An important note to consider is that to compute Shapley values, KernelShap needs to use a reduced number of features due to the time required to work, and also correlated features 106 6.2. RESULTS should be previously eliminated to avoid weighting the features wrongly, leading to an incorrect interpretation of the model performance. As mentioned above in geographic region classification topic and below in the methods section, correlated features were identified and treated. 107 CHAPTER 6. MACHINE LEARNING APPROACHES TO CLASSIFY HONEYS BASED ON FLORA AND GEOGRAPHIC ORIGIN FROM SANTA CATARINA, BRAZIL USING 1H-NMR 6.2.5.1 Features analysis of the best models performance for Geographic Regions classification The computed mean Shap values resulted in twenty features an average range of impact from 0-0.04. The results for the mean Shap values and their average impact on model output are shown in Figures 28, 29, 30, 31, 32, 33. The first seven computed mean Shap values with high importance and their impact on the models are shown in Table 22, for the calculated regions. By looking at these results the impact and order of features are different in each region model. 1A and 4A shared the feature 3.9ppm both with the first order of importance. However, in 1A their impact is higher and is positive compared to region 4A where the shap value for the feature is reduced and with a negative impact on the model. Regions 1B and 5 also shared a common feature with the highest importance in the model performance, in both models’ region classifications feature 5.41ppm has a positive impact but the impact is higher in 1B than in the model to classify region 5. The remaining region classification models have different features at the top of important features with different impacts on their classification models. 108 6.2. RESULTS Figure 28: Top 20 most important features ranked by mean absolute SHAP values. The features are ordered according to their importance of the best model classification for each production region. The Mean Shap values and their average impact on model output is shown for each feature. A: SHAP values and the average impact of feature’s on model of LR to classify 1A. B: SHAP values and the average impact of feature’s on model of LR to classify 1B. Table 22: Best seven Mean absolute SHAP values and their impact on model output of the best models classification for Geographic Regions. Note: (-): Negative impact on the model output; (+): Positive impact on the model output. LR LR PLS-DA RF LR SVC PLS-DA SVC Order 1A 1B 2C 3B 3C 4A 4B 5 1st 3.9 5.41 3.69 5.31 3.56 3.9 3.99 5.41 mean value 0.025 0.040 0.030 0.025 0.035 0.0005 0.05 0.009 impact (+) (+) (-) (+) (+) (-) (+) (+) 2nd 5.31 3.47 3.47 5.41 3.99 3.84 5.41 3.5 mean value 0.018 0.027 0.022 0.020 0.027 0.00045 0.035 0.0089 impact (-) (+) (-) (+) (-) (-) (-) (+) 3rd 3.75 1.34 4.17 3.75 1.34 5.41 2.07 3.56 mean value 0.017 0.018 0.022 0.018 0.023 0.00044 0.030 0.0069 impact (-) (-) (çç) (+) (-) (+) (+) (-) 4th 4.96 3.72 2.07 3.9 3.58 3.47 3.56 3.84 mean value 0.0165 0.0125 0.020 0.0125 0.022 0.00043 0.025 0.0069 impact (-) (-) (+) (+) (-) (+) (-) (+) 5th 4.04 3.78 3.56 3.84 5.41 4.96 3.9 3.26 mean value 0.016 0.0120 0.020 0.012 0.022 0.00025 0.025 0.0068 impact (-) (-) (+) (-) (-) (-) (-) (+) 6th 2.07 4.21 3.9 4.17 4.31 2.07 1.34 3.9 mean value 0.0125 0.011 0.017 0.011 0.021 0.00022 0.024 0.0067 impact (-) (-) (-) (-) (+) (-) (-) (-) 7th 3.84 3.58 3.72 4.31 3.9 3.23 3.23 3.58 mean value 0.0120 0.010 0.016 0.010 0.021 0.00021 0.020 0.005 impact (+) (+) (-) (+) (+) (-) (+) (+) 109 CHAPTER 6. MACHINE LEARNING APPROACHES TO CLASSIFY HONEYS BASED ON FLORA AND GEOGRAPHIC ORIGIN FROM SANTA CATARINA, BRAZIL USING 1H-NMR Figure 29: The 20 most important features are shown and ordered according to their importance to the best model classification for each production region. The Shap values and their impact on model output is shown for each feature. Each point represents a sample in the sub-dataset, and the color of the points represents the value of a particular feature for that sample. A: SHAP values and impact of feature’s on model of LR to classify 1A. B: SHAP values and impact of feature’s on model of LR to classify 1B. Figure 30: Top 20 most important features ranked by mean absolute SHAP values. The features are ordered according to their importance of the best model classification for each each production region. The Mean Shap values and their average impact on model output is shown for each feature. A: SHAP values and the average impact of feature’s on model of PLS-DA to classify 2C. B: SHAP values and the average impact of feature’s on model of RF to classify 3B. C: SHAP values and the average impact of feature’s on model of LR to classify 3C. 6.2.5.2 Features analysis of the best models performance for Flora groups classification Regarding flora groups with highest ROC-AUC classification models, Shapley values calculated an average range of impact from 0-0.025. The results for the mean Shap values and their average impact on model output are shown for each feature of Flora groups with higher ROC-AUC classification models in Figures 34, 35. The computed mean Shap values resulted in twenty features. The impact and order of importance are different in each flora 110 6.2. RESULTS Figure 31: The 20 most important features are shown and ordered according to their importance to the best model classification for each production region. The Shap values and their impact on model output is shown for each feature. Each point represents a sample in the sub-dataset, and the color of the points represents the value of a particular feature for that sample. A: SHAP values and impact of feature’s on model of PLS-DA to classify 2C. B: SHAP values and impact of feature’s on model of RF to classify 3B. C: SHAP values and the impact of feature’s on model of LR to classify 3C. Figure 32: Top 20 most important features ranked by mean absolute SHAP values. The features are ordered according to their importance of the best model classification for each each production region. The Mean Shap values and their average impact on model output is shown for each feature. A: SHAP values and the average impact of feature’s on model of SVC to classify 4A. B: SHAP values and the average impact of feature’s on model of PLS-DA to classify 4B. C: SHAP values and the average impact of feature’s on model of SVC to classify 5. model. The first seven computed mean Shap values with high importance and their impact on the models are shown in Table 23, for the calculated Flora groups. No common feature with importance at the first seven positions was observed in Flora group classification models. 111 CHAPTER 6. MACHINE LEARNING APPROACHES TO CLASSIFY HONEYS BASED ON FLORA AND GEOGRAPHIC ORIGIN FROM SANTA CATARINA, BRAZIL USING 1H-NMR 6.5.2 Methods 6.5.2.1 1D Nuclear Magnetic Resonance spectroscopy (1H-NMR) A 50 mg honey sample was added 700 uL D2O, mixed (1 min, vortex) and filtered on a cellulose membrane under reduced pressure. Before 1H-NMR analysis, 3-trimethylsilyl propionic-2, 2, 3, 3-d4 acid sodium salt (TSP, 98 atoms% D, 0.024g%) was added to samples as internal reference standard. The 1H-NMR analysis of the samples was performed at Brazilian Biosciences National Laboratory (Brazilian Center for Research in Energy and Materials, Campinas, SP), in a NMR spectrometer operating at a Larmor frequency of 499.726 MHz for 1H. The 1H-NMR spectra were recorded on a Varian Inova spectrometer, using the VNMRJ software. All experiments were performed at 300 K and non-spinning. For the high-resolution 1HNMR spectra 32K data points (time domain) were collected, with an acquisition time 4.8 s, delay time 1.0 s, 4 dummy scans, and total number of scans 32. Chemical shifts for protons were reported in parts per million (ppm) downfield from TSP at δ(1H) 0.00 ppm and referenced to that internal reference standard. 6.5.3 Statistical Analysis and Chemometrics Pre-processing spectral signal was performed by R language through the package specmine [8], available in CRAN. The remaining processing and data analyses were performed by using the Python language, and the libraries pandas, mumpy, ResearchPy, Scipy, Statsmodels, and Scikit-learn. 6.5.3.1 Spectral Data Processing and Data pre-processing The Peak alignment grouped proximal peaks together according to their median position across all samples, using a moving window of 0.03ppm. 2E4 was the baseline selection limit for the detection of minimum positive values. Typical resonance of the deuterium oxide (D2O)(4.80ppm) signal was removed from the dataset. No missing values and no outliers were detected. The processed dataset of the 56 honey samples was extracted into a file (comma separated values format .csv), where the first column indicates peak positions (ppm) and the second one represents the peak intensities containing a total of 36 peaks per sample, with an interval of peaks comprehending 1.34ppm to 5.41ppm. To note the the raw spectra without any treatment comprehends an interval of resonances between 0ppm to 13.78ppm. After spectral data processing, it was followed by z-score standardization, a selected treatment utilized to perform the data analyses. Two datasets were organized from the samples based on region of production and flora origin. To note that both datasets had the same samples but are organized according to class classification (production region or flora origin). All samples that had replicates were transformed by computing the respective spectra mean values. 118 6.5. MATERIALS AND METHODS 6.5.3.2 Data analysis Exploratory and statistical analysis A preliminary analysis included statistical univariate techniques and unsupervised multivariate statistical models to understand the samples relationships among the different classes according to the region of production and according to the flora origin. To understand if there are statistical differences between different combinations of must varieties an ANOVA with Tukey post hoc-test, between the each region pair and each flora pair, was carried to identify which specific resonances (ppm) could be utilized to discriminate between these classes. The unsupervised multivariate statistical models used in this study were Principal Component Analysis (PCA), Hierarchical Clustering Analysis (HCA) and T-Distributed Stochastic Neighbor Embedding (t-SNE), alternative approaches to interpret 1H-NMR data. PCA is widely used tool for interpreting and visualizing the main features, as for classifying samples and is this study was used to see if there is a clear discrimination between the samples with different classes in the latent space. HCA used an Euclidean distance between samples and the complete linkage method. This alternative was chosen because it tends to find compact clusters of approximately equal diameters and does not force clusters together due to single elements being close to each other, like in single linkage clustering. t-SNE is a method which applies a non-linear dimensionality reduction technique to keep the similar data points close in a lower-dimensional space. This method uses a student t-distribution to compute the similarity between two points in lower-dimensional space. Region and Flora classification and features analysis To classify each class a one-vs-all approach was followed, employing binary classification to separate a target class from any other set of honeys samples. The main objective is to classify a target class from any other source. Since this is an exploratory work with potential to add more classes (not regions but floras) in further studies, the binary classification is the best approach to apply into this scenario. The following models were used: Random Forest Classifier (RF), Gradient Boosting Classifier (GBC), Logistic Regression (LR), C-Support Vector Classification (SVC) and Partial Least Squares-Discriminant Analysis (PLS-DA). Several performance metrics were tested and compared, being the Receiver Operator Characteristic Under the Curve (ROC-AUC) the selected metric to evaluate the models performance. Due to the reduced number of samples, leave-one-out cross-validation was utilized to train and validate the model. Each model used hyperparameters is described in Table 24. 119 CHAPTER 6. MACHINE LEARNING APPROACHES TO CLASSIFY HONEYS BASED ON FLORA AND GEOGRAPHIC ORIGIN FROM SANTA CATARINA, BRAZIL USING 1H-NMR Table 24: Hyper-parameters used by the distinct models. Model Hyper-parameters Random Forest Classifier N of Estimators: 100; Max N of Features: sqrt; Max. of depth: None; Class weigth: balanced. SVC C-Support Vector Classification kernel: rbf; – Support Vector Machine Classification gamma: scale; C: 1; Class weigth: balanced Gradient Boosting Classifier N of Estimators: 100; Learning rate: 0.1; Max N of Features: None; Max. of depth: 3; Class weigth: balanced Logistic Regression C: 1.0; Partial Least Squares-Discriminant Analysis N of Components: 2 For each training, and validation set, the following preprocessing was carried: each model was calibrated using the treated raw spectra with the peak alignment, baseline selection and plus standardization using the aforementioned preprocessing pipeline. The calibration was performed with explicit feature selection. The explicit Feature Selection method consisted in dropping correlated features above a given threshold (in this work 90%). First, the correlation matrix of the spectral resonances was computed. Next, the correlation matrix was further processed by forming clusters of resonances with correlation values above the threshold. For each cluster, the resonance with the highest variance is kept and remaining ones were dropped. The calibrated models were evaluated in the test set, and the spectra treament was similar (using the mean and standard deviation values computed in the train dataset). Feature analysis To assess the feature importances, the best models from each class were chosen based on the best obtained performance in the aforementioned cross-validation schema for each target class (Region and Flora). Afterwards, each resonance influence on the target classes was inferred using the Shapley values. These values are the average marginal contribution of a feature value across all possible coalitions and can be applied for classification problems, as in our study [56], [58]. For binary classification, SHAP values for classes 0 and 1 are symmetrical. The selected feature contributes to a certain amount towards class 1, at the same time reduces the probability of being class 0 by the same amount. 120 7 Conclusions and Future Work The overall aim of this thesis was the development of pipeline approaches of statistical and machine learning models using metabolomics datasets to solve/improve different issues of food authentication problems. A variety of models were built and evaluated for three different authentication problems of two natural food products with a high-value market and PDO importance, mainly, wine, wine musts, and honey samples. Although the models developed in this work are still far from being used in food authentication processes, the work described here will be useful to guide the development of future computational methods for the design of novel and more accurate classification approaches. The developed work analyzed wine, wine musts and honey samples obtained from the Federal University of Santa Catarina and wine musts samples obtained from the Instituto Nacional de Investigação Agrária e Veterinária, Dois Polos, Portugal. Although the limited sample size of the metabolomics datasets posed some challenges, the results shed light on the potential of these tools to be use in the authentication of food products. Additionally, this study represents the first investigation focusing on honey samples from Santa Catarina Estate and the first to employ wine musts samples from Portuguese grape varieties grown in PDO regions of Portugal. Further studies are needed to expand on these preliminary findings. In the following sections, the main contributions of this thesis will be summarized, and suggestions for further improvement will be proposed. 7.1 Main contributions In Chapters 2 and 3, a review of the literature, regarding the state-of-the-art on metabolomics-based approaches for authentication of natural food products using statistical tools and machine learning algorithms, was performed. Multiple concepts across the metabolomics field were described and analyzed, including the main statistical and machine learning approaches used in food authentication processes. It became evident 121 CHAPTER 7. CONCLUSIONS AND FUTURE WORK that further development of novel model-based methods was needed to understand and overcome several hurdles still encountered in the reviewed studies. One of the main contributions of this work was the analysis of different methods tested on the classification performance of ML-based approaches to metabolomics data. A variety of preprocessing methods and modeling strategies were tested, in an attempt to understand which approaches produce the best results for each model classification. These contributions are present in the studies in chapters 4, 5, and 6. Chapter 4, performed the annotation of the metabolites present in wine samples with different ages of bottle storage, from different harvest years, and produced under similar climatic conditions. To classify each group of wine samples of the same age, different models were tested to select the one with higher performance with results between 80-92% of accuracy. To predict the influence of age on bottle storage and the climatic influences Linear Regression model was utilized. The used pipeline allows us to identify the resonances influenced by each of the climatic factors and age of bottle storage. This work is one of the few with a focus on the age of bottle storage using a dataset with a large time interval (years) and with the climatic influences over the production years. Our work proves the efficacy to combine metabolomics data combined with ML to classify problems related to harvest years and to predict the age of bottle and climatic influences in wine samples. In Chapter 5, ML was used to classify different Portuguese native wine musts varieties. During the development of a solution for this classification problem, it was necessary to deal with different wine must varieties samples resulting from the grape maturation process, from different but close production places. The implementation of feature selection methods combined with different models was assessed. This approach also included the test of different pre-processing treatments and different metrics evaluation, to understand which models were more suitable for the classification problem of the study. This pipeline is a possible approach to improve each variety model classification performance. This work also explored the interpretability of the model, being the first study in the authentication of wine musts samples, to use this method. With the use of the interpretation model, it was clear that based on our data that there are some resonance bands related to functional chemical groups which are associated with each variety classification. In Chapter 6, a similar analysis to chapter 5 was performed for the problem of botanical and geographic classification of honey produced in eleven agroecological regions from Santa Catarina, Brazil. Here the major difficulty was the number of regions with different sample sizes per region. Binary classification models were tested to deal with a large number of regions with different sample sizes. The implementation of feature selection methods combined with different models was well-suited to this type of classification task. And the use of different metrics evaluation, to understand which models were more adequate to the problem of the study was assessed. Here, it was selected the ROC-AUC metric, instead of MCC or PR-AUC, because it allows to understand the generalization of the model to each scenario. This pipeline was a good approach to improve each region’s model classification performance as to each botanical honey origin. As in the previous exploratory study, model interpretability, was explored. Based on the 122 7.2. FUTURE WORK available metabolomics data, it was possible to understand which resonances are associated with each region and botanical origin, and which resonances are related between regions and flora. These results can be used in future works to explore the authentication of honey based on both classifications and in examples of origin certifications. Until the moment of this dissertation submission, this is the first study of honey classification based on region and biological origin from Santa Catarina Brazil. The work of the three food products (wine, wine must, and honey) confirmed that ML algorithms are well-suited to the classification of authentication problems. The most important issues to pay attention are the quality of the data and what pre-treatments are more adequate for each metabolomics data. The use of feature selection methods, tested more than one model and evaluation metrics to get the most adequate model, able to interpret each authentication problem. Despite the ability of ML models to learn relevant features directly from raw data, the use of prior biological knowledge to select features is still useful to reduce the size of high-dimensional omics datasets, particularly the presence of solvents, water, and other elements present in IR and NMR datasets. The studies of food authentication have a reduced size of datasets and used similar and a reduced number of models to perform predictions or classifications. In this work, despite the reduced sample size for each study, several models were used. Distinct models may learn different aspects of the problem, even in cases with a reduced sample size. When tested more and different models combined with the proper validation model method, it may increase the performance of the final classification model. In this thesis, the issue of model interpretability was also addressed. Using the feature importance of the Shapley values method, we were able to demonstrate that the decisions made by our ML models are driven by relevant resonances related to functional chemical groups. The work undertaken in Chapter 4 contributed to extending and validating the Specmine software package, mainly giving support to metabolites annotation and metabolomics spectral data visualization developed by our research group. With these, reliable and reproducible analysis pipelines on a scientific level, the main objectives for this work were accomplished, although there are still many improvements that could be done in future works. 7.2 Future work The insights gained from this research can be used as a starting point for the development of better classification models for food authentication in the future. In this section, we suggest ways to extend and improve upon the work described in this thesis. A major limitation of our work was the reduced size of all datasets, with reduce number of samples not allowing to explore more issues of authentication process and limitation the validation of the models. 123 CHAPTER 7. CONCLUSIONS AND FUTURE WORK With this being said, the developed data analyses with the data available conditions allow us to perform reliable and reproducible analysis pipelines, on a scientific level. There are still, however, many aspects that could be improved in the web data analysis despite the size of dataset, including the improvement of already existing models implementation. The future work could encompass: • Combining multiple metabolomics datasets to have access to more training data, for example by using existing integrative databases; • Combining multiple models could be a good strategy to improve performance. Distinct models may learn different aspects of the problem, and thus combining them may increase the generalization capacity of the final model; • Preprocessing methods and evaluation metrics improve models performance in Chapters 4, 5 and 6. Thus, we recommend the use of more elaborate preprocessing and evaluation metrics and in the future food authentication studies. • Prior knowledge for feature selection proved to be beneficial in Chapter 5 and 6. In the future, other ways of leveraging this knowledge should be explored, including the use of alternative resonances for feature selection and methods that can incorporate chemical knowledge directly within the models. • Finally, the development of a more general and easy-to-use pipeline to food authentication processes with ML and Model interpretability could be applied to other wine varieties, different honeys and other food products in general. These pipelines could be give important contributions to the field of food authentication and would allow our work to be more widely used by the community. 124 Bibliography [1] D. I. Ellis, H. Muhamadali, S. A. Haughey, C. T. Elliott, and R. Goodacre, “Point-and-shoot: Rapid quantitative detection methods for on-site food fraud analysis–moving out of the laboratory and into the food supply chain”, Analytical Methods , vol. 7, no. 22, pp. 9401–9414, 2015. [2] D. J. Shaw, “World food security”, A History since , 1945. [3] M. Castro-Puyana, R. Pérez-Míguez, L. Montero, and M. Herrero, “Application of mass spectrometrybased metabolomics approaches for food safety, quality and traceability”, TrAC Trends in Analytical Chemistry , 2017. [4] A. Zhang, H. Sun, P. Wang, Y. Han, and X. Wang, “Modern analytical techniques in metabolomics analysis”, Analyst , vol. 137, no. 2, pp. 293–300, 2012. [5] O. Fiehn, “Metabolomics–the link between genotypes and phenotypes”, in Functional Genomics , Springer, 2002, pp. 155–171. [6] P. O. Larsen and M. Von Ins, “The rate of growth in scientific publication and the decline in coverage provided by science citation index”, Scientometrics , vol. 84, no. 3, pp. 575–603, 2010. [7] U. Roessner and J. Bowne, “What is metabolomics all about?”, Biotechniques , vol. 46, no. 5, p. 363, 2009. [8] C. Costa, M. Maraschin, and M. Rocha, “An r package for the integrated analysis of metabolomics and spectral data”, Computer methods and programs in biomedicine , vol. 129, pp. 117–124, 2016. [9] J. M. Cevallos-Cevallos, J. I. Reyes-De-Corcuera, E. Etxeberria, M. D. Danyluk, and G. E. Rodrick, “Metabolomic analysis in food science: A review”, Trends in Food Science & Technology , vol. 20, no. 11, pp. 557–566, 2009. [10] C. A. Georgiou and G. P. Danezis, Food authentication: Management, analysis and regulation . John Wiley & Sons, 2017. [11] S. G. Villas-Boas, J. Nielsen, J. Smedsgaard, M. A. Hansen, and U. Roessner-Tunali, Metabolome analysis: An introduction . John Wiley & Sons, 2007. 125 BIBLIOGRAPHY [12] K. O’Shea and B. B. Misra, “Software tools, databases and resources in metabolomics: Updates from 2018 to 2019”, Metabolomics , vol. 16, no. 3, pp. 1–23, 2020. [13] K. Böhme, P. Calo-Mata, J. Barros-Velázquez, and I. Ortea, “Recent applications of omics-based technologies to main topics in food authentication”, TrAC Trends in Analytical Chemistry , vol. 110, pp. 221–232, 2019. [14] R. Ramautar, A. Demirci, and G. J. de Jong, “Capillary electrophoresis in metabolomics”, TrAC Trends in Analytical Chemistry , vol. 25, no. 5, pp. 455–466, 2006. [15] D. S. Wishart, “Metabolomics: Applications to food science and nutrition research”, Trends in Food Science & Technology , vol. 19, no. 9, pp. 482–493, 2008. [16] S.-A. Sansone, T. Fan, R. Goodacre, J. L. Griffin, N. W. Hardy, et al. , “The metabolomics standards initiative”, Nature biotechnology , vol. 25, no. 8, pp. 846–848, 2007. [17] S. G. Villas-Bôas, U. Roessner, M. A. Hansen, J. Smedsgaard, and J. Nielsen, “Microbial metabolomics: Rapid sampling techniquesto investigate intracellular metabolite dynamics–an overview”, Metabolome Analysis: An Introduction , pp. 203–214, 2006. [18] A. Alonso, S. Marsal, and A. Julià, “Analytical methods in untargeted metabolomics: State of the art in 2015”, Frontiers in bioengineering and biotechnology , vol. 3, p. 23, 2015. [19] J. L. Markley, R. Brüschweiler, A. S. Edison, H. R. Eghbalnia, R. Powers, et al. , “The future of nmr-based metabolomics”, Current opinion in biotechnology , vol. 43, pp. 34–40, 2017. [20] A. Owen, “Fundamentals of uv-visible spectroscopy”, 1996. [21] E. Smith and G. Dent, Modern raman spectroscopy: A practical approach . John Wiley & Sons, 2013. [22] E. Acuna and C. Rodriguez, The treatment of missing values and its effect on classifier accuracy . Springer, 2004, pp. 639–647. [23] V. Hodge and J. Austin, “A survey of outlier detection methodologies”, Artificial intelligence review , vol. 22, no. 2, pp. 85–126, 2004. [24] J. L. Schafer, Analysis of incomplete multivariate data . CRC press, 1997. [25] A. R. Kapil, Methods of missing value treatment and their effect on the accuracy of classification models , 2018. [26] C. B. Costa, “Development of an integratedcomputational platform for metabolomics data analysis and knowledge extraction”, PhD thesis, 2014. [27] G. A. Morris, “Nmr data processing”, Encyclopedia of Spectroscopy and Spectrometry , pp. 125– 133, 2017. [28] O. M. Kvalheim, F. Brakstad, and Y. Liang, “Preprocessing of analytical profiles in the presence of homoscedastic or heteroscedastic noise”, Analytical Chemistry , vol. 66, no. 1, pp. 43–51, 1994. 126 BIBLIOGRAPHY [29] R. A. van den Berg, H. C. Hoefsloot, J. A. Westerhuis, A. K. Smilde, and M. J. van der Werf, “Centering, scaling, and transformations: Improving the biological information content of metabolomics data”, BMC genomics , vol. 7, no. 1, pp. 1–15, 2006. [30] K. H. Liland, “Multivariate methods in metabolomics–from pre-processing to dimension reduction and statistical analysis”, TrAC Trends in Analytical Chemistry , vol. 30, no. 6, pp. 827–841, 2011. [31] A. Savitzky and M. J. Golay, “Smoothing and differentiation of data by simplified least squares procedures.”, Analytical chemistry , vol. 36, no. 8, pp. 1627–1639, 1964. [32] K. Varmuza and P. Filzmoser, Introduction to multivariate statistical analysis in chemometrics . CRC press, 2016. [33] Å. Rinnan, F. Van Den Berg, and S. B. Engelsen, “Review of the most common pre-processing techniques for near-infrared spectra”, TrAC Trends in Analytical Chemistry , vol. 28, no. 10, pp. 1201– 1222, 2009. [34] S. Kus, Z. Marczenko, and N. Obarski, “Derivative uv-vis spectrophotometry in analytical chemistry”, Chem. Anal , vol. 41, no. 6, pp. 889–927, 1996. [35] J. Godzien, A. G. de la Fuente, A. Otero, and C. Barbas, Metabolite annotation and identification . Elsevier, 2018, vol. 82, pp. 415–445. [36] P. Patil, “What is exploratory data analysis”, Toward Data Science , 2018. [37] P. Filzmoser, K. Hron, and C. Reimann, “Univariate statistical analysis of environmental (compositional) data: Problems and possibilities”, Science of the Total Environment , vol. 407, no. 23, pp. 6100–6108, 2009. [38] T. W. MacFarland and J. M. Yates, Introduction to nonparametric statistics for the biological sciences using r . Springer, 2016. [39] R. B. Darlington and A. F. Hayes, Regression analysis and linear models: Concepts, applications, and implementation . Guilford Publications, 2016. [40] L. Van der Maaten and G. Hinton, “Visualizing data using t-sne.”, Journal of machine learning research , vol. 9, no. 11, 2008. [41] M. Rocha, P. Cortez, and J. M. Neves, Análise inteligente de dados: Algoritmos e implementação em java . 2008. [42] Y. Saeys, I. Inza, and P. Larranaga, “A review of feature selection techniques in bioinformatics”, bioinformatics , vol. 23, no. 19, pp. 2507–2517, 2007. [43] V. Kumar and S. Minz, “Feature selection”, SmartCR , vol. 4, no. 3, pp. 211–229, 2014. [44] R. Kohavi, Special issue on applications of machine learning and the knowledge discovery process . Kluwer, 1998. 127