A metabolomic data fusion approach to support gliomas grading

Magnetic resonance imaging (MRI) is the current gold standard for the diagnosis of brain tumors. However, despite the development of MRI techniques, the differential diagnosis of central nervous system (CNS) primary pathologies, such as lymphoma and glioblastoma or tumor‐like brain lesions and glioma, is often challenging. MRI can be supported by in vivo magnetic resonance spectroscopy (MRS) to enhance its diagnostic power and multiproject‐multicenter evaluations of classification of brain tumors have shown that an accuracy around 90% can be achieved for most of the pairwise discrimination problems. However, the survival rate for patients affected by gliomas is still low. The High‐Resolution Magic‐Angle‐Spinning Nuclear Magnetic Resonance (HR‐MAS NMR) metabolomics studies may be helpful for the discrimination of gliomas grades and the development of new strategies for clinical intervention. Here, we propose to use T2‐filtered, diffusion‐filtered and conventional water‐presaturated spectra to try to extract as much information as possible, fusing the data gathered by these different NMR experiments and applying a chemometric approach based on Multivariate Curve Resolution (MCR). Biomarkers important for glioma's discrimination were found. In particular, we focused our attention on cystathionine (Cyst) that shows promise as a biomarker for the better prognosis of glioma tumors.

Magnetic resonance imaging (MRI) is the current gold standard for the diagnosis of brain tumors. However, despite the development of MRI techniques, the differential diagnosis of central nervous system (CNS) primary pathologies, such as lymphoma and glioblastoma or tumor-like brain lesions and glioma, is often challenging. MRI can be supported by in vivo magnetic resonance spectroscopy (MRS) to enhance its diagnostic power and multiproject-multicenter evaluations of classification of brain tumors have shown that an accuracy around 90% can be achieved for most of the pairwise discrimination problems. However, the survival rate for patients affected by gliomas is still low. The High-Resolution Magic-Angle-Spinning Nuclear Magnetic Resonance (HR-MAS NMR) metabolomics studies may be helpful for the discrimination of gliomas grades and the development of new strategies for clinical intervention. Here, we propose to use T 2 -filtered, diffusion-filtered and conventional waterpresaturated spectra to try to extract as much information as possible, fusing the data gathered by these different NMR experiments and applying a chemometric approach based on Multivariate Curve Resolution (MCR). Biomarkers important for glioma's discrimination were found. In particular, we focused our attention on cystathionine (Cyst) that shows promise as a biomarker for the better prognosis of glioma tumors. age. 1 Primary brain tumors are mainly gliomas, with a number of histology subtypes: glioblastoma is the glioma of grade IV, which is the most invasive and frequent; anaplastic astrocytoma and anaplastic oligodendroglioma are gliomas of grade III and low-grade astrocytoma and low-grade oligodendroglioma are diffuse gliomas of grade II. Gliomas are classified into grades I to IV by the World Health Organization (WHO) using histopathological and clinical criteria. The new WHO 2016 classification prognosis is influenced by grading and biomolecular assessment in IDH mutation and MGMT methylation. 2,3 About 60-70% of the primary brain tumors are glioblastomas and their prognosis is still poor, with a median of 12-14 months of overall survival despite the new surgical, oncological and radiotherapeutic techniques.
Magnetic resonance imaging (MRI) is the current gold standard for the diagnosis of brain tumors. However, despite the development of MRI techniques, the differential diagnosis between CNS primary pathologies such as lymphoma and glioblastoma or tumour-like brain lesions and glioma or recurrent tumors and radionecrosis is often difficult. 4 Further studies in neuro-oncology for the diagnosis and care of brain tumors are thus needed. MRI can be flanked by in vivo magnetic resonance spectroscopy (MRS) in order to enhance its diagnostic power and multiprojectmulticenter assessments of classification of brain tumors have shown that an accuracy around 90% can be achieved for most of the pairwise discrimination problems. 5 In vivo MRS provides, in principle, important metabolic information on tumors, detecting signals such as N-acetylaspartate (NAA), choline containing compounds (ChoCC) and creatine (Cr) and, in good quality spectra, also glutamate (Glu), glutamine (Gln), lactate (Lac), alanine (Ala) and other more overlapped signals from myo-inositol (Myo), glycine (Gly) and taurine (Tau). These metabolites are often aspecific and do not represent clear-cut markers for the different neuropathologies, even though recently Myo signals have been reported to be significantly higher in high-grade gliomas with respect to CNS lymphomas and lipid signals to be higher in gliomas with respect to other brain tumors. 6 A complementary approach can be represented by ex vivo High-Resolution Magic-Angle-Spinning Nuclear Magnetic Resonance Spectroscopy (HR-MAS NMR) which enriches and strengthens the interpretation bases of in vivo spectra. The number of identifiable metabolites found in ex vivo spectra using HR-MAS are about 10-fold those recognizable in in vivo spectra, due to higher resolution. For instance, HR-MAS NMR clearly shows, using two-dimensional techniques (COSY, TOCSY and HSQC), that the in vivo ChoCC peak receives contributions not only from choline (Cho), phosphocholine (PC) and glycerophosphocholine (GPC), but also from phosphoethanolamine (PE), (Myo, ß-glucose (Glc), Tau and arginine (Arg). 7,8 The ex vivo HR-MAS NMR spectra most amenable for multivariate analysis are those T 2 -filtered, obtained applying a Carr-Purcell-Meiboom-Gill (CPMG) sequence. They retain signals from small metabolites undergoing fast motion, and hence characterized by narrow linewidths, whereas signals from slowly tumbling species (usually lipids and macromolecules) that give broad spectral components, partially observed as baseline distortions, are filtered off. Detailed studies on different grade gliomas based on statistics applied to CPMG spectra have been reported. 9,10 Nevertheless, T 2 -filtered signals represent only a part of the spectroscopic information, the part related to slowly diffusing molecules can be found in diffusion-filtered experiments and all the NMR visible metabolites contribute to the basic water-presaturated pulse-and-acquire spectra.
In this work, we evaluate not only T 2 -filtered spectra but we also try to extract as much information as possible from diffusion-filtered and conventional water-presaturated spectra fusing the data gathered by these different NMR experiments and then applying the Multivariate Curve Resolution (MCR) approach on each experiment individually. MCR coupled with an aligning method is particularly important on tissues, since the chemical shifts of the signals of the same metabolite can slightly change in the different samples, for instance due to slight pH variations. As these changes cannot be controlled by using buffers, as it is usually done in NMR on fluids, the alignment step can be used to avoid this potential issue.
Aims of this study are i) to investigate the different glioma grades and the biomarkers most important for their discrimination, taking advantage of the details that can be obtained by ex vivo HR-MAS NMR, in order to support the interpretation of in vivo MRS; ii) to gain further insight into the biochemistry of the tumor subtypes and iii) to help in the identification of metabolic alterations in biofluids, to be implemented in the screening for early diagnosis and for therapy planning.

| Patient population
This study was approved by the local ethics committee (CE 131/10, approved in 07/09/2010). All patients received detailed information regarding the procedure and written informed consent was obtained. Based on pre-surgical MRI parameters, image-guided tissue samples (n = 44) were collected from 35 enrolled subjects by biopsy or surgical resection, as summarized in Table 1.
The pre-surgery MRI (data not reported in this paper) were acquired and used to plan the tissue sample locations prior to surgery. Samples were collected from a specific tumoral area identified by MR images (the same analyzed by in vivo MRS) and transferred to the neuronavigation system (BrainLAB Inc., USA). Each sample was split into two parts; one was signed and used for standard histology, the other for HR-MAS NMR measurements. The last one was snap-frozen in liquid nitrogen and stored at -80 C. The diagnosis of grade was assessed by a single pathologist on the basis of standard histological criteria. At the time of sample collection, the mutations in isocitrate dehydrogenase 1 and 2 (IDH1 and IDH2) were not assessed. Mutations in IDH1 and IDH2 and complete codeletion of chromosome arms 1p and 19q have been recognized as favorable prognostic molecular markers. As these markers have represented a major breakthrough in the diagnosis of brain tumors, they have been integrated into the 2016 WHO classification of gliomas. 3,11 The IDH mutation analysis was performed only for some patients with recurrent cancer.

| HR-MAS Nuclear Magnetic Resonance Spectroscopy
HR-MAS NMR spectra were recorded with a Bruker Avance400 spectrometer (Bruker BioSpin, Karlsruhe, Germany) operating at a frequency of 400.13 MHz for 1 H and 100.61 for 13 C. The instrument was equipped with a 1 H, 13 C HR-MAS probe whose temperature was controlled by a Bruker Cooling Unit. The samples, 4-70 mg, were placed into a MAS zirconia rotor (4 mm outer d), added with 5-20 μL of D 2 O to provide deuterium for the lock system, and the volume was adjusted to 50 μL (or 12 μL in the case of very small samples) using the proper insert and then transferred into the probe cooled to 5 C to prevent tissue degradation processes. The set up takes about 20 min. Experiments were performed at a temperature of 5 C, spinning the samples at 4 kHz and three different types of one-dimensional proton spectra were acquired 12 using the sequences implemented in the Bruker software: (i) a composite pulse sequence (zgcppr, hereafter zg) with 2 s water presaturation during the relaxation delay, 8 kHz spectral width, 32k data points, 32 scans, (ii) a water-suppressed spin-echo CPMG sequence (cpmgpr, hereafter cpmg) with 2.5 s water presaturation during the relaxation delay, 1 ms echo time (τ) and 360 ms total spin-spin relaxation delay (2nτ), 8 kHz spectral width, 32k data points, 64 scans; and (iii) a sequence for diffusion measurements based on stimulated echo and bipolar-gradient pulses with longitudinal eddy current delay (ledbpgp2s1d, hereafter led) with big delta 200 ms, eddy current delay Te 5 ms, little delta 2*2 ms, sine-shaped gradient with 32 G/cm followed by a 200 μs delay for gradient recovery, 8 kHz spectral width, 8k data points, 256 scans. Assignments of 1 H signals were based on previous findings, 7 checked with standard two-dimensional homonuclear and heteronuclear techniques (COSY, TOCSY and HSQC) and compared with data from open-access databases, mainly HMDB (http://www.hmdb.ca/ version 4.0) and BRMB (http://www.bmrb.wisc.edu/).

| Multivariate Data Analysis
Three data sets, corresponding to the three types of NMR spectra which were acquired, were considered, namely: led (4096 data points), cpmg (15000 data points) and zg (7500 data points). The data analysis workflow consisted of the following steps: (i) preprocessing of each single data set; (ii) features extraction applying MCR; (iii) exploratory data analysis of each single data set (not reported for the sake of brevity) and on the fused data set by Principal Component Analysis (PCA); (iv) classification of the different glioma grades by two different approaches, namely classmodeling applying Soft Independent Modelling of Class Analogies (SIMCA) and Discriminant Analysis applying Partial Least Squares-Discriminant Analysis (PLS-DA). Given the limited number of samples, a double cross-validation scheme was applied for model validation. The salience of metabolites in differentiating the tumor grades were assessed using the discriminant power 13 in SIMCA and Variables Important in Projection (VIPs), 14 regression vectors and selectivity ratio 15,16 in PLS-DA.
Each step is described in detail in the following sub-sections.

| NMR data preprocessing
The three different groups of spectra, namely zg, cpmg and led, were processed using the Bruker software, applying line broadening of 0.5, 0.5 and 5 Hz, respectively, prior to Fourier transformation, and then phase-corrected. The zg and cpmg spectra were referenced to the chemical shifts of the CH 3 doublet of alanine at δ 1.48 ppm, whereas the led spectra were referenced to the CH 3 peak of fatty acid (FA) chain at δ 0.89 ppm.
The zg NMR data points were reduced from 16k to 8k. The selected spectral width spanned from δ 9 ppm to δ 0.05 ppm and the final number of data points was 7500. For the cpmg spectra, the NMR data points were reduced from 32k to 16k. The selected spectral width spanned from δ 9 ppm to δ 0.05 ppm and the final number of data points was 15000. Finally, for the led spectra, the NMR data points were reduced from 16k to T A B L E 1 Tissue samples and patient population 8k. The selected spectral width spanned from δ 9 ppm to δ 0.05 ppm and the number of final data points was 4k. All NMR data were exported in ASCII format and imported in Matlab using a routine written by us.
Preprocessing of NMR data was then performed on each spectral data set individually. The procedure consisted of three steps: denoising, baseline correction and peak alignment. Denoising was performed in wavelet (WT) domain 17 applying a global threshold criterion 18,19 to the detail coefficients obtained by WT decomposition of the NMR spectra using a Daubechies 4 wavelet filter and 5 as decomposition level.
Baseline correction was then performed on the reconstructed denoised spectra by means of weighted asymmetric least-squares, 20 fitting a fifth order polynomial. Next, the denoised and baseline-corrected spectra were aligned, first globally (i.e. rigidly shifting the whole spectrum) and then operating on a set of carefully defined intervals, which were individually aligned.
To this aim, the icoshift algorithm 21 was used. The option 'average' 22 was selected both for the global and interval alignments: this setting performs the alignment procedure two times in a row, leading to better aligned peaks.

| Features extraction and data fusion
MCR 23 was successfully applied to extract features from the NMR spectra, allowing to obtain a total of 95 resolved components from the three data blocks, of which 49 were resolved from the zg spectral dataset, 27 from cpmg and 19 from led. A brief description of MCR is given hereinafter, while a detailed description can be found in the literature. 24 MCR is a decomposition method based on Beer's Law which allows to obtain a pure spectral profile S i and a vector of relative concentrations C i (one concentration value for each sample) for each ith-extracted component. The pure spectral profiles carry the information about the signal's shape and position and can be easily used for matching signals with metabolites by comparison with literature, spectral library data and previous knowledge.
MCR can recover pure contributions from overlapping signals applying the decomposition equation reported in Figure 1 and solving it by alternating least squares, starting from an initial guess of either matrix C (concentrations) or S (spectral profiles). 23 Given the complexity of the NMR spectrum, it has been demonstrated that working on intervals instead of the whole spectral width can improve interpretability and model performances. Differences due to baseline noise, signal intensity and metabolites can be more easily handled, and meaningful chemical quantifications can be obtained from each region. 25 An initial spectral profiles S matrix was estimated by the SIMPLISMA algorithm, 26 , 27 with 10% maximum noise allowed. The non-negativity constraint was applied to both the concentrations and the profiles, using in both cases the fast non-negative least squares algorithm. 28 Finally, the S matrix was normalized using the Euclidean norm. The maximum number of iterations was set to 500 and, whenever convergence was not reached within this limit, the evolution of the spectral profile of each resolved component was inspected, to confirm that the components of interest (those related to NMR signals) were stable enough at the last iteration.
For each defined interval a distinct local MCR model was developed: the S profiles of the resolved components were inspected to select meaningful spectral profiles, while suspicious and baseline-like profiles were discarded. The identification of the resolved components was based on literature assignments, the digital libraries of HMDB and BMRB and our knowledge about the shape and position of the signals. 7 F I G U R E 1 Extraction of chemical features from NMR spectra through interval-based MCR modeling For each dataset, the concentration vectors corresponding to the identified components were organized into a new data matrix (MCR extracted features). The features matrices were then normalized by the sample's mass, to account for differences in signals' intensity that could affect the correct sample characterization.
Once each data set (zg, cpmg and led) was resolved, the extracted features, named after the identified metabolites, were combined in a new dataset. This approach of combining features extracted from different sources is referred to as mid-level data fusion 29,30 and all the results shown hereafter were obtained from the data-fused matrix (represented in blue in Figure 1).

| Exploratory analysis
Exploratory data analysis was performed by PCA 31,32 on the autoscaled (i.e. the mean value was subtracted from each dataset's column, which was then divided by its standard deviation) mid-level fused data set containing glioma samples of grade IV and III. Results are presented in the Results and Discussion section.

| Classification analysis
Taking into account both the clinical relevance and the limited number of samples, especially of grade III, two distinct models were derived: (i) a one-class model for assessing if a glioma is of grade IV using SIMCA and (ii) a PLS-DA model to discriminate between grade IV and grade III gliomas. The basis of the two methods are briefly summarized below.

| SIMCA
SIMCA 13,33 belongs to the class-modelling methods, in which each class is individually modeled, independently from the presence of other classes.
The focus of class-modeling is on intra-class similarities, i.e. similarities among samples belonging to the same class, and the classification rule consists of establishing if a sample belongs to the modeled class or not. To capture characteristic class variability, SIMCA builds a distinct PCA model for each class being modeled (preprocessing is also applied distinctly to each class) and the number of PCs is established independently for each class and may differ for each of them.
The classification rule is then based on verification of whether a sample is accepted or rejected by the given class model in terms of distance of the sample to the class model. Two distances are defined: (i) the scores distance (T 2 ), which indicates how far a sample is from the training objects of the class in the PC space, and (ii) the orthogonal distance, which measures the distance of a sample to the PC space of the class and is given by the sum of squared residuals (Q). The implementations of SIMCA 34 may differ with respect to how these two distances are combined in an overall distance measure and on how the statistical acceptance limit, for this measure, is defined. In this work we used the implementation of the PLS_Toolbox (see Software section) that considers distinct statistics for T 2 and Q and bases the classification rule on the reduced distance D r : In this implementation, a sample is accepted if its D r from a given class is smaller or equal to the square root of 2 or rejected if it is higher. The classification performance of the model was evaluated in terms of class sensitivity, i.e. percentage of samples belonging to the class, which are correctly accepted, and specificity, i.e. percentage of samples not belonging to the class, which are correctly rejected.
The class dimensionality was established according to the minimum classification error in leave-one-out cross-validation, which corresponds to a one-component model. Validation of the model was carried out using Double cross-validation, 35 as described below.
Interpretation of the SIMCA model in terms of salient metabolites for discriminating the samples with respect to the modeled classes can be accomplished calculating the variable discriminatory power. 13 This parameter is calculated comparing the square residual standard deviations of a variable calculated for the objects when projected on the model of the class they do not belong to with respect to that of the class they belong to. The concept is that we expect the residuals of samples on a given variable being higher when fitted to a different class model with respect to their own, if this variable is important in discrimination. Therefore, this variable will have a higher value of discriminatory power. Since the discriminatory power is not upper bounded, we have used both the 95 percentile and the average of the calibration samples, as thresholds to assess the significance of each variable. Only to calculate the discriminatory power a SIMCA model for the grade III class was also evaluated.

| PLS-DA
PLS-DA is a discriminant method 36 in which discriminant boundaries among classes are computed. As opposed to SIMCA, this method focuses on the differences between the classes, so the classes are not individually modeled and therefore at least two classes have to be defined. Moreover, samples will necessarily be assigned to one of the modeled classes. PLS-DA is based on PLS regression 14 where the dependent variables are dummy variables (one for each modeled class) taking values of 0 if the sample does not belong to the class and 1 when it belongs to it. Implementation of PLS-DA may differ on the basis of how the classification rule is defined. Here we use the rule to assign a sample to the class for which the predicted y-value is the highest (i.e., Y predicted values are continuous and not dummy as they were codified).
The number of PLS-DA components was established according to the minimum classification error in leave-one-out cross validation, which corresponds to a two-components model. The performance of the model was established in terms of non-error classification rate, i.e. percentage of samples correctly assigned to the respective classes.

| Model validation: Double cross-validation
Double cross-validation (double CV) is an iterative procedure based on two cross-validation loops, one "external" and one "internal". The internal loop is used to optimize the complexity of the model, while the external one is used to obtain prediction values for each one of the samples. The dataset is initially split into balanced segments, and it is processed according to these steps:1 one segment j is excluded; 2 the remaining j-1 segments are used to build a series of models using a leave-one-out CV scheme (internal loop), to compute the optimal complexity of each jth model; 3 once the optimal complexity is obtained, the jth model is built and the segment excluded in step 1 is predicted (external loop); 4 a set of prediction errors is stored.
Using this procedure, each sample is excluded and predicted one time, so one prediction value for each sample is obtained. These prediction errors are used to estimate the predictive capability of the model. A final calibration model is then built whose complexity (number of components) was established taking the median of the optimal complexity (established by the internal CV loop) of the j models.
Considering the presence of different number of replicates for each sample, here the internal and external segments have been customized, taking always replicates in the same segment. Moreover, in the case of PLS-DA representativeness of both classes in each segment was ensured.
Summarizing: the grade IV and III samples were divided into four balanced splits for the external loop, and for SIMCA modeling, only grade IV samples were considered in the modeling step.

| Software
The whole data analysis process was carried out on MATLAB 2017b (Mathworks, MA, USA). PCA analysis was performed using the PLS_Toolbox 8.6 (Eigenvector Research Inc. WA, USA) and WT denoising using the Matlab Wavelet Toolbox 19 . NMR spectral alignment was operated using icoshift, 21

| RESULTS AND DISCUSSION
A total of 35 subjects were evaluated, of which 3 patients with glioma II, 8 with glioma III, 19 affected by glioma IV and 5 patients with preliminary diagnosis of glioma and then reclassified after histological analysis as lymphomas. More than 70% of patients died between the first and second year after surgery. Despite numerous treatment strategies, these data show that the survival rate is still low and the NMR metabolomics studies may be helpful for discrimination of gliomas grades and the development of new strategies for clinical intervention. Three one-dimensional experiments were recorded for each sample (Figure 2). The first one is a conventional 1 H HR-MAS NMR spectrum with water presaturation (zg), which allows to detect the whole NMR visible metabolome of intact tissues, formed by both rapidly and slowly tumbling molecules that give narrow and broad lines, respectively. This is the real and complete fingerprint of a tissue. The second one is a T 2 -filtered spectrum, obtained through a cpmg sequence that filters off broad resonances, characterized by short spin-spin relaxation times (i.e. T 2 ) and highlighting signals from small metabolites (mainly monosaccharides and polyols, aminoacids, organic acids and osmolites). The last one is a diffusion-edited led spectrum that retains the broad signals coming from slowly moving species (macromolecules and lipids) at expense of the narrow ones due to rapidly diffusing metabolites.
The majority of HR-MAS NMR studies reported in the literature relies on cpmg spectra, as they show resolved signals and flat baselines.
Despite the predictable difficulties implicit in the processing of water-presaturated zg and diffusion-edited led spectra, we believe that restricting the statistical approach to cpmg data could limit the spectroscopic information that can be in principle gathered from tissues. We then tried to derive metabolic and grading knowledge from the combined statistical analysis of the three types of spectra. To this aim, zg, cpmg and led spectra were processed for each sample as described in the experimental section.

| Exploratory analysis
A PCA model to explore the data was built prior to the classification step. Figure 3a shows a clear separation between grade III and IV samples, with some important exceptions. Three grade IV samples were found to be either extreme (04), extreme-borderline (28) or lying among the grade III samples (17). Two grade III samples also had interesting positions: sample 33-33_2 (replicates of the same specimen) was found in borderline position, while sample 22 was lying among the grade IV samples. This can be better seen in Figure 3b.
The separation evidenced by the dash-dotted red curve (a guide for the eye) in Figure 3 correlates with neuroradiological characteristics of the two groups of tumors, in particular with the presence of necrosis and marked contrast enhancement in grade IV and with the absence of necrosis and the reduced contrast enhancement in grade III. Samples 17 and 28 were classified by histological analysis as grade IV, but possess neuroradiological characteristics similar to those of grade III tumors, i.e. the absence of necrosis and the reduced contrast enhancement typical of grade III.
Grade II gliomas and lymphomas were projected on this PCA model and they are also shown in Figure 3a. Lymphoma samples seem to be most similar to grade IV samples, while grade II gliomas are more similar to grade III gliomas.
It is evident that sample 04 is extreme, also compared to the other grade IV samples. Inspecting the corresponding T 2 -contribution plot ( Figure S1 reported in the Supplementary Materials), which shows the contribution of the original variables in determining the scores values of the samples, it emerges that its extreme position is mainly due to the significantly higher contributions from FAs. Sample 04 was obtained from a tumor lesion characterized by a cystic necrotic area, which could induce the high presence of FAs. Moreover, this sample, at the histological level, has features of atypical cells such as gemistocytic cells that define a different type of astrocytoma, that relapses and has a poor prognosis.
Based on this analysis, three grade IV samples (04, 17,28) and one grade III sample (22) were not considered in the calibration set for the classification step. A summary of spectra considered in classification models is presented on Table 2.

| Classification analysis: SIMCA
Lymphomas and samples belonging to grade II were not considered in the classification modeling step, together with three samples of grade IV and one sample of grade III, as reported on Table 2. The total number of samples included for classification modelling was 23 (32 spectra), F I G U R E 2 1 H NMR spectra of a grade IV sample obtained with a zg sequence (A), a cpmg sequence (B) and a led sequence (C). The first one displays both narrow and broad signals, the second one retains the narrow signals from fast tumbling molecules, whereas the third one only the broad resonances from slowly tumbling species. Major metabolites are labelled: alanine (Ala), 2-aminoadipic acid (2-aa), adenine, adenosine, ascorbate This allowed us to estimate sensitivity and specificity of the grade IV model in prediction using the classification assignments of the external loop. These results are reported in Table 3 and Figure 4.
The model wrongly rejected only three grade IV samples out of 21 (03_1, 16,19). It is important to clarify that samples 03_1 and 16 lie within the acceptance limit in Figure 4, because this figure refers to their estimation when they are present in the calibration set. However, the above two samples were rejected during the external cross-validation assessment, which is the information to be used when referring to predictive capability for grade IV samples.
The model wrongly accepted only one grade III sample (33) out of 11. It is interesting to notice that sample 33 was found in PCA (Figure 3b) very close to grade IV samples, together with its replicate 33_2, which was correctly rejected by the SIMCA model, but can be found very close to the acceptance limit in Figure 4. Samples 33 and 33_2 derive from the same cancer tissue, and their difference is probably due to the heterogeneity of the tissues. The grade II glioma and lymphoma samples listed in Table 3 were also predicted by the SIMCA model, and the prediction results are reported in Table 4. All the grade II samples were correctly rejected, but two lymphoma samples (23,34) were wrongly accepted.
A comparison of the metabolites' relative content (peak areas obtained by MCR) of lymphoma samples 23 and 34 with respect to the correctly rejected ones reveals that the two rejected samples have significantly lower content of all metabolites. Moreover, sample 05 was correctly rejected but it can be found close to the acceptance limit (Figure 4), and its profile resulted more similar to the two wrongly accepted samples (23,34) than to the other two rejected ones (24,31). Recently, a paper has been published about lymphoma tumor. 38 This paper focuses on the discrimination between malignant lymphomas and gliomas, and the authors report that it is usually difficult to preoperatively distinguish between the two. The use of MRS can be useful for preoperative diagnoses, and quantitative analysis is considered to be a valuable method for distinguishing between gliomas and malignant lymphomas. The limitation of the cited 38 and of our study is the number of samples.
As expected for their position in exploratory PCA model, the three glioma grade IV samples 04, 17, 28 ("excluded" in Table 2) were wrongly rejected and the glioma grade III sample 22 wrongly accepted. This last finding can be explained by the final diagnosis of 22 as xanthoastrocytoma.

| Classification analysis: PLS-DA
Only samples belonging to grades III and IV were considered in this classification modeling step. As opposed to the SIMCA classification, it is not possible to predict samples belonging to classes different from grade III and IV using the PLS-DA model, as explained in the Multivariate data F I G U R E 4 Reduced Orthogonal distance (Q/Q lim ) vs. reduced Scores distance (T 2 /T 2 lim ) of the SIMCA model built for the grade IV gliomas. The red circle highlights the acceptance area, which is delimited by the reduced distance threshold, computed using the formula reported in the paragraph about SIMCA, under the Experimental Section. Symbols are reported in the legend; the symbols with lighter colors represent the twelve projected samples i.e. 3 gliomas of grade II, 5 lymphomas, 1 glioma of grade III and 3 of grade IV, listed in the last column of analysis section about PLS-DA. As with the SIMCA classification, instead of dividing the dataset into a training and a test set, the double CV approach was used. The model performance results are reported in Table 4.
As expected, the samples' distribution in the scores space of the PLS-DA model ( Figure 5) resembles the one found in PCA (Figure 3), with a clear separation between the two classes. It is important to clarify the reason why the model wrongly predicted three samples, while Figure 5 shows a perfect separation between the two classes. As for the SIMCA model, the information obtained from the double CV is the one to be used when referring to the model's predictive capability. Figure 5 refers instead to the model built using all the samples altogether, which generally performs better than the single models built in cross-validation.
It is interesting to inspect the position of the wrongly predicted samples. These samples result somehow extreme within their own class ( Figure 5, labelled samples).
The PLS-DA model wrongly assigned sample 33_2 to grade IV, as opposed to the SIMCA model, which wrongly assigned sample 33 (a replicate of sample 33_2) to grade IV. These two replicates can be found close to each other in Figure 5, but also close to the grade IV samples. As already underlined, the difference found for these adjacent samples is probably due to the heterogeneity of the glioma sample. Furthermore, subject 33 was initially diagnosed with an oligodendroglioma that became an anaplastic tumor after a few years.

| Significant metabolites
In this section the SIMCA and PLS-DA sets of discriminant metabolites are compared. The shared metabolites are reported on Figure 6.
Twenty-eight discriminant signals were recovered by SIMCA, thirty-six were recovered by PLS-DA and twenty-one were found in common between the two classification methods. Six signals related to FAs were found, and the main source for them were the led NMR spectra.
Maleschlijski and co-workers 9 report that grade IV compared to grade III tumors have a higher content of FAs, and the same trend was found in our study, as it can be seen in Figure 6. In the same study, Cr and Myo are reported to be less abundant in grade IV compared to grade III tumors: this trend was also found in our study and it can be seen in Figure 6 as well. 9 Other metabolites such as NAA, 39 Cyst and acetate (Ac) follow the same descending trend ( Figure 6).
These trends are also confirmed by Figure 7, where the discriminant metabolites important for grade IV (bars in blue) are mainly FAs, which are more abundant for this grade, while the other discriminant metabolites are important for grade III (bars in green) and are also more abundant for this grade.
The role of FAs in the distinction between lymphomas and grade IV gliomas is still an open question. It has been suggested that their signals can be useful to distinguish between malignant lymphomas and gliomas, selecting regions without necrosis, for malignant lymphomas and glioblastomas have different mechanisms for the generation of FAs. 38 In the case of gliomas, a higher FA peak is due to cystic or necrotic components, 40  For this reason, we report in Figure 8 and 9 the fingerprints of Cyst and in Table S1 its chemical shifts (together with those of some selected metabolites). In details, the methylene protons of Cyst at 3.12 ppm correlate with a carbon at about 35 ppm, and not at 44 ppm as in ethanolamine, and the methylene protons of Cyst at 2.73 ppm correlate with a carbon at 30 ppm, and not at 58 ppm, as in the case of HTau. In some samples, both HTau and Cyst can be found and distinguished. In our samples instead, ethanolamine was present in small amounts and Cyst signals usually dominated 1 H spectra in the region around 3.1 ppm, as can be seen from Figure 8. The T 2 filter acts more on Cyst than on HTau, the signals of which become more evident in the cpmg spectrum.
One recently published study has reported the first measurement of Cyst by in vivo MRS. 44 The identification of Cyst was confirmed comparing in vivo spectra acquired gliomas (with IDH mutations) with the Cyst spectrum measured in a phantom. In the present study we support these findings and the direct assignment of Cyst in vivo MRS.
Finally, we want to stress the importance of verifying 1D 1 H NMR assignments through selected 2D experiments, when possible. A signal (triplet) around 2.25 ppm can be attributed both to 2-aminoadipate and to 2-hydroxyglutarate (2OHG) 10 and the assignment can be done on the F I G U R E 8 1 H NMR water-presaturated spectra of a GIII sample obtained with a zg sequence (A) and a cpmg sequence (B). Signals from Cyst are marked by red squares. Signals from HTau are marked by blue squares. Those in broken line are highly overlapped to other signals F I G U R E 9 2D NMR COSY (A) and HSQC (B) spectra of a GIII sample with Cyst correlations marked basis of TOCSY spectra (see Figure S2 and Table S1). TOCSY correlations are very helpful in this respect and allowed us to assign it to 2-aminoadipate in seven samples and to 2OHG only in one.
Cyst and 2OHG are considered important metabolites in cancer pathologies. Their role is related to IDH mutations. IDH mutations can be found in gliomas and are characterized by a specific cellular metabolism, causing the accumulation of 2OHG in tumor cells. 45,46 Higher expression of cystathionine-β-synthase (CBS), the first enzyme of the transsulfuration pathway, has been associated with better prognosis in In IDH-mutated 1p/19q codeleted gliomas. 47 Cyst derives from the condensation of homocysteine with serine catalyzed by CBS, which is the initial and rate-limiting step in the transsulfuration pathway. Cyst is subsequently cleaved by the enzyme cystathionine gamma-lyase (CTH) to form cysteine, a precursor of GSH. Moreover, CBS participates in the desulfuration reactions that contribute to endogenous hydrogen sulfide production. Deregulation of CBS and the associated alterations in homocysteine and/or hydrogen sulfide levels leads to a wide range of pathological disturbances, and CBS activity also plays an important but complex role in cancer biology. 48 Reduced serine biosynthesis may lead to increased reliance on the CBS/CTH pathway as a critical response to increased oxidative stress. 47,49 In particular, high CBS expression has been shown to confer better prognosis in IDH-mutated 1p/19q codeleted gliomas 47 in line with a previous study showing that decreased expression of CBS promotes glioma tumorigenesis in tumor xenografts. 50 Correlation between higher CBS expression and survival in IDH-mutated 1p/19q codeleted gliomas has been confirmed also from the POLA public dataset. 44

| CONCLUSIONS
This study highlights metabolic differences between anaplastic astrocytomas (grade III gliomas) and glioblastomas (grade IV gliomas), as identified at the time of the first diagnosis. The application of MCR allowed to efficiently extract important and meaningful metabolic details from the brain tumour samples obtained using HR-MAS NMR. This approach is particularly crucial on tissues, since slight changes in the chemical shifts of the signals of the same metabolites are commonly observed in the different specimens, and they cannot be controlled using buffers, as usually done in NMR on fluid samples.
The three one-dimensional spectra datasets under examination (zg, cpmg and led) were analysed and modeled together, within the framework of data fusion.
Two classification models of different nature (SIMCA, class-modelling vs PLS-DA, discriminant analysis) were built, leading in both cases to good classification rates as reported in Tables 3 and 4. Moreover, the two approaches provided valuable information on metabolites that can be utilized for the distinction between grade III and grade IV gliomas. Among these metabolites, Cyst surely is the most interesting result. Cyst and other related amino acids such as 2OHG identified in this study are relatively new in the MRS and HR-MAS NMR panorama and they seem to be good candidates as markers to monitor brain cancer prognosis and treatments.
Once again, we would like to underline that the use of 2D NMR spectra is essential for the correct assignment of metabolites with similar chemical shifts and signal shapes.

AKNOWLEDGEMENTS
Centro Interdipartimentale Grandi Strumenti of Università di Modena e Reggio Emilia is greatly acknowledged for the use of Bruker Avance400 spectrometer.