Discovering Glioma Tissue through Its Biomarkers’ Detection in Blood by Raman Spectroscopy and Machine Learning

The most commonly occurring malignant brain tumors are gliomas, and among them is glioblastoma multiforme. The main idea of the paper is to estimate dependency between glioma tissue and blood serum biomarkers using Raman spectroscopy. We used the most common model of human glioma when continuous cell lines, such as U87, derived from primary human tumor cells, are transplanted intracranially into the mouse brain. We studied the separability of the experimental and control groups by machine learning methods and discovered the most informative Raman spectral bands. During the glioblastoma development, an increase in the contribution of lactate, tryptophan, fatty acids, and lipids in dried blood serum Raman spectra were observed. This overlaps with analogous results of glioma tissues from direct Raman spectroscopy studies. A non-linear relationship between specific Raman spectral lines and tumor size was discovered. Therefore, the analysis of blood serum can track the change in the state of brain tissues during the glioma development.


Introduction
The most commonly occurring malignant brain tumors are gliomas, and among them is glioblastoma multiforme (GBM), which accounts for 14.3% of all tumors and 49.1% of malignant tumors [1][2][3]. GBM is the most aggressive, invasive, and undifferentiated type of tumor, and has been designated grade IV by the World Health Organization [4,5]. GBM is among the deadliest neoplasms.
One reason for the poor outcome of glioblastoma is a late-stage diagnosis, since most of the existing methods of noninvasive diagnostics, such as magnetic resonance imaging (MRI) [6][7][8] and computer tomography [9], are ineffective for diagnosing small-size tumors. Optical methods, such as fluorescence imaging [10], multiphoton microscopy [11], photoacoustic imaging [12], optical absorption spectroscopy [13][14][15], Raman spectroscopy and imaging [16,17], and terahertz (THz) imaging [18][19][20] are widely used for direct detection of glioma tissues spectral fingerprints, including glioma's molecular biomarkers, under cancer tissue surgery resection. Similar methods are very sensitive to tissue chemical content, but have too low a penetration depth to in vivo noninvasive glioma detection. Tissue biopsies can be taken to study cellular and molecular composition. It is an invasive and traumatic procedure.
Mouse blood serum was used for spectroscopic measurements. The study was conducted on 60 male mice of the immunodeficient SCID line at the age 6-7 weeks. The model of orthotopic xenotransplantation of the U87 human glioblastoma cells into mice was used [49,50]. A total of 5 µL of the U87 MG cell suspension (500,000 cells per one injection) was introduced in the subcortical brain structure through a hole in the animal's cranium. Animals from the control group were injected in a similar manner with 5 µL of the culture medium. The tumor size was measured using a horizontal 11.7 Tesla MRI tomograph (Biospec 117/16; Bruker, Billerica, Waltham, MA, USA) [49]. Each experimental group had a corresponding control group. There were 10 mice in each group. The animals were removed from the experiment by decapitation on days 7, 14, and 21 after injection. Blood samples were collected, serum was separated, and frozen at −80 • C. The study was carried out in accordance with the EU Directive 2010/63/EU and the ARRIVE 2.0 guidelines, and approved by the Inter-Institutional Commission on Biological Ethics at the Institute of Cytology and Genetics, Siberian Branch of the Russian Academy of Sciences (Permission #78, 16 April 2021).

Raman Spectroscopy
DXR Raman microscope (Thermo Fisher Scientific, Waltham, MA USA) with a magnification of 10×, excitation wavelengths of 532 nm, and range 80-4000 cm -1 was used. Each sample of blood serum was a droplet with a volume of 10 µL placed on a special aluminum plate [51,52]. This plate had identical holes in the form of a funnel with a diameter of 5 and a depth of 2 mm (see Figure S1). This form of the plate leads to a uniform distribution of the sample upon drying. We can assume that the distribution of samples in each well and the film thickness of the sample after drying will be the same. The measurements were carried out after complete drying of the sample, which took place for 20 min at room temperature. The Raman spectrum of dried blood serum was recorded at three points in the center of the well and averaged (see Figures S2 and S3). Each spectrum was averaged over 300 scans to remove cosmic ray noise.

Machine Learning Pipeline
The implemented machine learning pipeline is shown in Figure 1. At preprocessing step, Raman spectra were normalized by maximal value, in order to be able to compare spectral intensities. Also, we applied adaptive iteratively reweighed penalized least squares method to remove background noise [53]. A Python library BaselineRemoval contains several methods to remove the noise, caused by fluorescence in Raman spectra, and according to our previous findings [54], the airPLS algorithm [51] has a good performance. Principal component analysis (PCA) was used to extract informative features [55,56]. Support vector machine (SVM) [57], random forest (RF) [58], and XGboost [59] were used to build prediction data models. The ten-fold cross-validation was used [60]. The At preprocessing step, Raman spectra were normalized by maximal value, in order to be able to compare spectral intensities. Also, we applied adaptive iteratively reweighed penalized least squares method to remove background noise [53]. A Python library Base-lineRemoval contains several methods to remove the noise, caused by fluorescence in Raman spectra, and according to our previous findings [54], the airPLS algorithm [51] has a good performance. Principal component analysis (PCA) was used to extract informative features [55,56]. Support vector machine (SVM) [57], random forest (RF) [58], and XGboost [59] were used to build prediction data models. The ten-fold cross-validation was used [60]. The prediction data model performance was estimated using a confusion matrix (see Table 1 for reference), where true positive is the number of correctly diagnosed U87 mice, true negative is the number of control mice labeled as healthy, false negative and positive values are respective errors of classification. The matrix members can be combined using receiver operating characteristic-area under curve (ROC AUC) based on true positive rate (TPR) versus false positive rate (FPR) values [61]. Here, TPR = TP/(FP + TN), and FPR = FP/(FP + TN). The ROC curve is a plot of TPR (vertical axis) versus FPR (horizontal axis). Whenever the area under the ROC curve is close to 1 or 0, the data model has a perfect quality. If it is close to 0.5, then the performance is equal to random guessing. Combination of k-fold cross-validation and ROC AUC analysis can visualize mean and standard deviation of ROC curve [60].

Characteristics of Experimental Animals
In the dynamics of the U87 glioblastoma development, the tumor size was 2.6 ± 0.4, 10.6 ± 1.8, and 89.6 ± 11.5 mm 3 (means ± standard error) on days 7, 14, and 21 after tumor cells' injection, respectively. Figure 2 shows that there is a non-linear increase in tumor size from day 7 to day 21 of the experiment.

Raman Spectra Analysis
Below, in Figure 3, we present preprocessed Raman spectra of U87 and control groups for all weeks of the experiment.

Raman Spectra Analysis
Below, in Figure 3, we present preprocessed Raman spectra of U87 and control groups for all weeks of the experiment.   As depicted in the Figure 3, standard deviation of Raman spectra decreases over time, and this can be related to the tumor growth [50,52]. The major changes in brain tissue metabolism are associated with tumor growth in U87-3 groups, hence, we expect to see regular, rather than random, peaks on the Raman spectra (see Figure 3c).
The groups' separability was visually estimated by PCA method. Thorough analysis is required to plot all possible combinations of more than 3000 principal components (PCs) (see Figure S4 in Supplementary Materials). Therefore, we focus only on 10 principal components, which cover more that 75% of variance. The visualization of the best group's separability by PCA is presented in Figure 4. As depicted in the Figure 3, standard deviation of Raman spectra decreases over time, and this can be related to the tumor growth [50,52]. The major changes in brain tissue metabolism are associated with tumor growth in U87-3 groups, hence, we expect to see regular, rather than random, peaks on the Raman spectra (see Figure 3c).
The groups' separability was visually estimated by PCA method. Thorough analysis is required to plot all possible combinations of more than 3000 principal components (PCs) (see Figure S4 in Supplementary Materials). Therefore, we focus only on 10 principal components, which cover more that 75% of variance. The visualization of the best group's separability by PCA is presented in Figure 4.
As shown in Figure 4, the tumor and control groups differ, yet have intersections. Interestingly, Raman data groups are separated in principal components planes with low explained variance. That points to the minor difference in relative intensity. Also, for the second week of the experiment, there is a good separability, yet some samples are closer to the control group. We believe that the progression of the glioma is lower for them.
Next, we tested separability of the tumor group versus control group for each week by: SVM, RF, and XGBoost with and without application of PCA. The results for the third week of the experiment are presented in Figures 5-7. The results for the first and second weeks of the experiment are presented in Figures S5-S10 of Supplementary Materials. Table 2 shows the total data for all groups and methods. regular, rather than random, peaks on the Raman spectra (see Figure 3c).
The groups' separability was visually estimated by PCA method. Thorough analysi is required to plot all possible combinations of more than 3000 principal components (PCs (see Figure S4 in Supplementary Materials). Therefore, we focus only on 10 principa components, which cover more that 75% of variance. The visualization of the best group's separability by PCA is presented in Figure 4.  As shown in Figure 4, the tumor and control groups differ, yet have intersections. Interestingly, Raman data groups are separated in principal components planes with low explained variance. That points to the minor difference in relative intensity. Also, for the second week of the experiment, there is a good separability, yet some samples are closer to the control group. We believe that the progression of the glioma is lower for them. to the control group. We believe that the progression of the glioma is lower for them.
Next, we tested separability of the tumor group versus control group for each week by: SVM, RF, and XGBoost with and without application of PCA. The results for the third week of the experiment are presented in Figures 5-7. The results for the first and second weeks of the experiment are presented in Figures S5-S10 of Supplementary Materials Table 2 shows the total data for all groups and methods.     We can see from Table 2 that SVM shows the best value for the U87-1 and U87-3 groups with the application of PCA, and for the U87-2 group without the application of PCA. For RF, we can see opposite patterns for the U87-1 and U87-2 groups compared to SVM. Application of PCA does not affect ROC AUC analysis for XGBoost.

Informative Feature Selection
The results of the informative feature selection for the third week of the experiment are shown in Figure 8. The same results for the first and second weeks can be found in Supplementary Materials (Figures S11 and S12).
All methods point to the Raman shifts near 1000 cm −1 and another one near 3000 cm −1 , yet RF and XGBoost marks 600 cm −1 and 2600 cm −1 , respectively. The variation can be explained by a high dimensionality of the data, so it should be verified by the respective mean spectra.
It should be noted that molecule or molecular band identification requires both central line position and band height/width parameters, and a good important feature selection algorithm should choose both. From this perspective, SVM feature ranking is reasonable, because Raman shifts near the central line are also significant. There is a special point about feature selection in random forests and extreme gradient boosting methods [62]. If variables are highly correlated, which is the case for the central line and surrounding band, the algorithm either chooses randomly between them (RF) or focuses on of them (XGBoost). What we can see in Figure 8 is not incorrect results, but specifics of the method, which should be taken into account, which is that is a variation of the central line position is possible. Additionally, we have filtered out features with low information gain for better representation.
We also measured tumor volume growth during the experiment (see Figure 2). As we can see, the change in volume is a non-linear function, so we applied a non-linear regression model based on RF (Python module sklearn.ensemble.RandomForestRegressor with default parameters except regularization parameter ccp_alpha was set to 0.015 to avoid overfitting), with Raman spectra as independent and tumor volume as dependent variable. The quality of RF regressor was estimated by mean and standard deviation of absolute error (MAE) and determination coefficient (R 2 ). Validation was performed by repeating ten-fold CV 20 times. This method also allows selecting the most informative features (see Figure 9). To avoid plotting more the 3000 features with importance values, we select only those with an absolute value greater than 0.02. Such a small value of importance is due to the large number of features. The results are the following: MAE = 25.3 ± 7.4, R 2 = 0.919. Informative Raman shifts are: 417, 420, 818, 1746, 2807, 3039, 3212, and 3563 cm −1 . Table 3 shows the correspondence between the Raman shifts for the constructed regression model, molecular vibrations, and the corresponding substances [44,[63][64][65][66][67][68][69][70][71][72][73][74][75]. The results of the informative feature selection for the third week of the experiment are shown in Figure 8. The same results for the first and second weeks can be found in Supplementary Materials (Figures S11 and S12). also allows selecting the most informative features (see Figure 9). To avoid plotting more the 3000 features with importance values, we select only those with an absolute value greater than 0.02. Such a small value of importance is due to the large number of features. The results are the following: MAE = 25.3 ± 7.4, R 2 = 0.919. Informative Raman shifts are: 417, 420, 818, 1746, 2807, 3039, 3212, and 3563 cm −1 . Table 3 shows the correspondence between the Raman shifts for the constructed regression model, molecular vibrations, and the corresponding substances [44,[63][64][65][66][67][68][69][70][71][72][73][74][75].

Discussion
The main idea of this paper is to estimate dependency between glioma tissue and biomarkers in blood serum, revealed by Raman spectroscopy. We used the most common model of human glioma, when continuous cell lines such as U87, derived from primary human tumor cells, are transplanted intracranially into the mouse brain [49,76]. This model is characterized by rapid glioma growth following the injection of tumor cells into the subcortical structures of the brain. We have a controlled growth of the tumor and can trace the change in blood composition in the dynamics of tumor development.
The modeling of glioblastoma with orthotopic transplantation of glioblastoma cells in the brains of laboratory mice is associated with many traits of traumatic brain injury syndrome [77,78]. It is especially clear in the first week of the experiment. The U87 cell suspension was introduced in the subcortical brain structure through a hole in the animal cranium. Animals from the control group were injected in a similar manner with 5 µL of the culture medium. When injected, brain injury occurs and an inflammation reaction develops, both in control and experimental mice. Therefore, in the first week of the experiment, the development of the inflammation reaction prevails over the development of the tumor [79,80]. The intensity of the characteristic peaks in the fingerprint range of 400-1700 cm −1 is higher in control animals and it is the same in high-frequency range of 2700-2950 cm −1 (see Figure 2a). Opposite trends are seen for the second and third weeks of the experiment (see Figure 2b,c): the intensity of the Raman peaks is higher in U87 mice.
We study the separability of the experimental and control groups for different glioma stages after the introduction of tumor cells by machine learning methods and try to discover most informative Raman spectral bands for this separation. SVM, RF, and XGboost methods were chosen because they have good generalization ability and additionally allow the selection of informative features. Figures 5-7 and S5-S10, as well as Table 2, show that SVM has a better generalization ability, so removing noise with PCA has a positive role for SVM. For RF and XGBoost, on the contrary, the fewer the number of features, the worse their learning. For group U87-3, quite good results are obtained with and without PCA, because there are obvious differences and PCA does not have an influence on them.
We would like to focus on the principal component (PC) selection for visualization. PCA sorts PCs in a descending order, so the first PCs correspond to the largest amount of variance in data. It is not necessarily related to the groups' separability [81]. In our case, the first several PCs do not affect group separability, so we do not describe them thoroughly.
At the next stage, we identified informative spectral features, according to which the experimental and control groups were divided at each stage of the experiment. The use of three methods makes it possible to establish differences with the highest degree of reliability (see Figure 8, Figures S11 and S12). As can be seen from the frequency distribution, Raman spectroscopy makes it possible to follow the change in the contribution of certain metabolites at each stage of the experiment with an increase in tumor size. It is shown earlier that attenuated total reflectance Fourier-transform infrared (ATR-FTIR) spectroscopy coupled with supervised learning methods allows for identifying glioma patients with tumor volumes of 0.2 cm 3 [40]. In our study, tumor volumes vary from 0.002 to 0.089 cm 3 at different stages of the experiment, that is, several orders of magnitude lower than in [40].
We applied a non-linear regression model to select the most informative Raman spectra features associated with tumor growth (see Figure 9 and Table 3). As can be seen from Table 3, the frequencies found by us are in good agreement with the Raman spectra of glioma tissues presented by other authors [44,[63][64][65][66][67][68][69][70][71][72][73][74][75]. It should be understood that the Raman spectrum of a sample is a superposition of vibrations of all substances present in the sample [15]. One vibration may be due to the contribution of several substances present in the sample. Vice versa, one substance can correspond to a complete set of peaks [82]. The discrepancy in the identification of informative Raman shifts can be explained by two reasons. First, it is computational. Machine learning uses optimization methods, hence, an approximate solution is sought and there can be several solutions. The stability of the approximate solution is improved by averaging. Such averaging can produce bias in informative Raman shifts [83]. Second, it is physical. Different physical conditions [84] and concentrations of substances, such as salt, can lead to changes in the position of Raman shifts [85]. It has also been shown that, depending on the glioma grade, the tissue microenvironment affects the shift of Raman bands by up to 38% [86]. According to Figure 9 and Table 3, with the progression of the tumor, an increase in the contribution of lactate, tryptophan, fatty acids, lipids, and water molecules is observed. Using the method of THz spectroscopy, which is known to be sensitive to the state of water in biological samples [87], we have previously shown that in the U87-3 group, the state of water changes, and the contribution of bound water molecules increases [79]. The same patterns are seen in the most informative Raman spectral features for each stage of the experiment. Thus, we show that the analysis of blood serum can track the change in the state of brain tissues during the development of glioblastoma.
The presented results confirm the possibility of glioblastoma detection by Raman spectra analysis of blood serum. This study has several limitations. First of all, the number of samples in each group is not large (10 per group). Though a proper validation of constructed machine learning models is made, a comparison to a reference method is required, namely, MRS, to determine metabolites in tumor and blood tests in order to find the quantity of these metabolites in blood. Second, in our previous work [50], we show that THz spectroscopy can be useful for the same task, so the combination of these modalities by ensemble machine learning techniques can be performed. Third, interesting data separation into groups of the mixed tumor and control samples can be seen on PCA visualization figures. This can be caused either by physiological processes of tumor development or by external factors, such as inflammation, due to the injection of the culture medium. The reason for such behavior will be investigated in further research.
The strong sides of this study are the following. The conducted experiments with a small animal glioblastoma model were performed in strict accordance with the research protocol and can be reproduced. State-of-the-art machine learning methods were used to verify the separability of the classes for each week of the experiment. The validity of the obtained results is tested by k-fold CV techniques and ROC AUC analysis. Raman spectra informative ranges were selected, which allowed us to speed-up further analysis. We can focus only on narrow bands, which reduces spectra acquisition time. A non-linear relationship between specific Raman spectral lines and tumor size was discovered. This information coupled with THz spectroscopy data allows the possibility of pointing to a metabolic profile associated with glioblastoma development.

Conclusions
The major goal of the work was discovering glioma through the detection of its biomarkers in the blood by Raman spectroscopy. In the work, the dried blood serum of mice was studied during the U87 glioblastoma development. It is a commonly accepted laboratory model of human glioma implemented under controlled conditions. Machine learning methods were applied to analyze the Raman spectra and identify the most informative frequencies. Thus, we show that the analysis of blood serum can track the change in the state of brain tissues during the glioma development. Raman spectroscopy combined with machine learning allows us to carry out label-free and real-time analysis of brain tissue and blood serum, as well as the differentiation of glioma grades and rapid monitoring of treatment.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.