Chondrogenic Cancer Grading by combining Machine and Deep Learning with Raman Spectra of Histopathological Tissues

1. Introduction

Chondrosarcoma (CS) turns out to be the second most common form of primary bone tumor worldwide, with almost seven cases per million people registered every year [1]. In the continuous effort to introduce effective and minimally invasive therapies, the ability to grade CS and to distinguish this form of bone cancer from healthy tissues and/or benign tumors, such as Enchondroma (EC), represents a crucial factor. The common protocol for the diagnosis of CS is based on an initial examination of radiological and/or Magnetic Resonance (MRI) images, followed by the histopathological analysis of a tissue biopsy. One of the drawbacks of this procedure is represented by the frequent discrepancies between the response of pre- and post-surgical histopathological analysis, mainly due to the presence of tissue inhomogeneities within the same tumor mass. This last feature also determines the presence of not well-defined tumor margins, undoubtedly affecting the outcome of the surgical therapies. As previously demonstrated [2], the presence of not removed tumor margins represents a concrete risk factor, increasing the probability of recurrence and thus decreasing the 5-year survival rate. Setups for image-guided surgery, e.g., Computed Tomography (CT) and MRI, were introduced to help the surgeon in maximizing the size of the excised tumor mass by preserving the surrounding healthy tissues. Despite the progress of these techniques, they still show weaknesses. In particular, CT is based on the use of ionizing radiation, with possible side effects for the patient and/or the operator. Furthermore, MRI employs intense magnetic fields, not suitable for patients with metallic devices implanted within the body. Among the imaging techniques, we also mention the Ultrasonographic (US) analysis. However, this technique is not yet considered a reliable solution, since it suffers from problems in the detection of soft tissues characterizing bone metastasis [3]. Finally, nuclear medicine, such as Positron Emission Tomography (PET), represents an interesting solution, since it potentially provides additional information about metabolism and drug activity. Despite the large values of sensitivity and specificity in detecting tumor masses, nuclear medicine is characterized by small spatial resolution. This feature pushed the researchers to couple such technique to other approaches, such as MRI [4].

Raman Spectroscopy (RS) represents a promising, easy-to-use, and low-cost solution to respond to the aforementioned needs [5,6,7]. RS is based on the measurement of the so-called Raman effect [8], i. e. a light scattering phenomenon produced by the interaction between a beam of photons and virtual molecular energy levels [9]. The difference between the wavelengths of the scattered and the incident photons is strongly correlated to the chemical properties of the molecular target that determined the scattering. Therefore, by performing a spectral analysis of the radiation scattered by a sample of interest, it is possible to retrieve detailed information about the chemical composition of the target under investigation [10]. Since the Raman effect doesn’t involve the interaction between photons and well-defined electronic energy levels, the technique is intrinsically not subject to non-radiative relaxation phenomena, possible source of sample cooling and consequent degradation. This feature makes RS particularly suitable to in vivo applications, e.g., for intra-operative cancer detection, aimed at estimating the tumor margins. Unlike the conventional histopathological analysis, requiring staining procedures to bring out the characteristic features of the tissues under investigation, RS is a label-free approach, i. e. applicable directly to the tissue with no preliminary treatments.

Despite the aforementioned strengths, RS presents two main limitations for its application in cancer diagnosis: first of all, the low acquisition speed could represent a relevant obstacle if diagnostic times of fractions of seconds are needed, e.g., in intraoperative cancer detection. In addition, the presence of a large amount of information characterizing a single Raman spectrum dramatically increases the difficulties in the visual interpretation of Raman spectra. This last issue assumes particular significance in the study of biological tissues, often characterized by complex chemical composition. This class of samples is characterized by frequent cases of subtle molecular differences between tissues corresponding to different diseases affecting the same type of organ, often associated with partially overlapped Raman bands.

In this sense, the employment of Machine Learning (ML) algorithms to interpret Raman-derived experimental data could represent a powerful solution to retrieve the information of interest in times compatible with the final use requirements [11,12,13]. In particular, Several works [14,15,16,17] showed how RS in combination with ML allows to detect malignant tissues in timescales of minutes or even seconds, being a potential help in the intraoperative estimation of the tumor margins.

Analogously, recently, RS has been applied to chondrogenic tumor classification with excellent classification performances, by employing either classical ML or algorithms belonging to the sub-field of ML called Deep Learning (DL) [6,18,19]. Such investigations demonstrated the effectiveness of such algorithms in systematically interpreting Raman-derived experimental data for the diagnosis of bone tumors. This finding represents an advantage over the analysis of stained tissue sections, requiring expert pathologists to express an ultimate diagnosis. In addition, the high level of detail of such an approach could potentially detect sub-types or sub-stages of a disease, intermediate to the conventionally accepted ones, e.g., the three different grades of CS. This outcome could represent a breakthrough toward an accurate oncological diagnosis and, therefore, the planning of personalized therapies.

In this paper, we explored the ability of ML and DL models to perform an accurate diagnosis of CS and EC from Raman-based data, collected from ex-vivo tissues. Our analysis was focused either on the problem of distinguishing CS from EC or on assigning a correct grade to the CS cases. The results demonstrated the effectiveness of ML and Deep Learning (DL) models in identifying the aforementioned CS tumor, with classification accuracy reaching 99.7%. This analysis of the Permutation Feature Importance (PFI) allowed to identify Raman spectral components, relevant for the classification process. In particular, the analysis of PFI potentially represents an important tool in the attempt to restrict the spectral range, to reduce the acquisition time and the computational cost of the classification phase.

2. Materials and Methods

2.1. Samples

The study involved ten patients, treated in 2018 in the Azienda Ospedaliero Universitaria Pisana, Pisa, Italy. Among the examined patients, three of which corresponded to the diagnosis of EC, three to CS of grade 1 (G1), two to CS of grade 2 (G2) and two to CS of grade 3 (G3).

Further details about the patients’ cohort are available in reference [6]. The resulting bone-excised tissues were subjected to formalin-fixing and paraffin-embedding, without decalcification. Two tissue sections of thickness 5-$\mu \mathrm{m}$ were obtained from the resulting excised masses for each patient. After being deposited on a glass slide, commonly employed in microscopy, the paraffin was removed through immersion in two baths of xylene for 10 min. Then, the residuals of formalin were removed by rinsing in Polybutylene succinate (PBS). One of the aforementioned two sections was stained with Hematoxylin and Eosin (H&E) for the following histopathological examination. This procedure allowed to assign a label to the resulting Raman spectra and therefore to employ ML in a supervised fashion. The other tissue section was not subjected to staining and it was employed for the acquisition of Raman spectra. The employment of not stained tissues for the acquisition of Raman spectra allowed to minimize the contribution to the resulting signal attributable to the sample fluorescence.

2.2. Raman Apparatus and Data Pre-Processing

Here, we will provide only a brief description of the Raman apparatus employed in this paper and of the pre-processing operations carried out on the resulting spectra. Further details are available in [6]. The Raman setup employed in this work was a Thermo Fisher Scientific DXR2xi [20] Confocal Raman Microscope, whose setup is schematically represented in Figure 1. The Raman setup was equipped with a 532 nm laser, whose power was set within the range between 5-10 mW, representing a right compromise between the need to avoid tissue damage and to achieve an adequate signal-to-noise ratio.

The back-scattered photons were collected with a $100\times$ objective and the out-of-focus components of the detected light were suppressed through an aperture (pinhole) of diameter $25\mu \mathrm{m}$. A 1200 gr/mm diffraction grating allowed to separate the spectral components from the detected light. We reported a schematic representation of the Raman setup in Figure 1. Each spectrum was collected with 10 accumulations and an acquisition time of 0.2 s for each accumulation. The resulting Raman spectra had a spectral resolution of $\sim 2{\mathrm{cm}}^{-1}$ and they were defined within the spectral interval between 400 and $1800{\mathrm{cm}}^{-1}$. In Figure 2, we showed images of H&E stained sections associated to the four types of tissue examined in this work.

The spectra were collected in grids of pixel resolution $4\mu \mathrm{m}\times 4\mu \mathrm{m}$, within the portions of unstained samples. The final dataset included $N_s=400$ spectra, 80 spectra for EC, 80 for G1, 84 for G2 and 93 for G3. We applied a baseline subtraction algorithm to the raw Raman spectra through the tool Omnicxi, integrated into the application controlling the CRM. We carried out this operation by subtracting a 5th-order polynomial. Then, we referred the spectra to the minimum signal. We didn’t normalize the spectra, to exploit the additional information of the absolute signal intensity for the following ML and DL classification routines.

2.3. Data Analysis

In the following, we will represent the Raman data as couples $\left\{(x_i,y_i)\right\}$, $i=1,...,N_s$, where $x_i\in {\mathbb{R}}^{N_p}$ is a vector storing the Raman intensities within a single spectrum, while $y_i$ represents the value of the label associated to the i-th spectrum. The spectrum $x_i$ can be also viewed as the i-th row of a matrix X of size $N_s\times N_p$.

As stated in the Introduction section, the main objective of this work was to employ RS in combination with ML to classify bone tissue. In particular, we focused our attention on three distinct classification problems:

• The problem of distinguishing EC and CS (EC-CS);
• The problem of distinguishing G1, G2 and G3 (G1-G2-G3);
• The problem of distinguishing EC, G1, G2 and G3 (EC-G1-G2-G3).

To this aim, we employed three ML protocols:

• Random Forest Classifier (RFC): this non-linear ML algorithm is based on building decision trees by training them on datasets obtained by bootstrap-sampling spectra from the initial dataset. The result of this procedure is a “forest”, whose prediction is based on the majority of the responses of the trees belonging to it. In our analysis, we adopted a forest of 4000 trees, to minimize the so-called out-of-bag error [21];
• Multi-Layer Perceptron (MLPC): in this investigation, this simple DL algorithm consisted of a single hidden layer with 900 neurons. We introduce non-linearity through a ReLU activation function. We carried out the training either with the ADAM [22] or with the L-BFGS-B [23] solvers, with an upper limit of 600 iterations. In the following, we will refer to the aforementioned DL models as MLPC(ADAM) and MLPC(L-BFGS-B), respectively;
• Support Vector Machine (SVM): this ML algorithm is aimed at determining the so-called maximum-margin hyperplane, separating the vectors $\left\{x_i\right\}$ corresponding to different values of the label [24]. The general equation of a hyperplane can be written as

  $$w^{T}x-b=0,$$

  (1)

  where w is a vector normal to the hyperplane and b is a real constant. According to the linear version of SVM, the classification is performed by solving the following minimum problem:

  $$\underset{w,b}{min}{w}^{T}w$$

  (2)

  with the constraint $y_i(w^T{x}_i+b)\ge 1$$\forall i=1,\dots ,N_s$. Despite the introduction of non-linear and more advanced versions of this algorithm, we adopted the original linear version as representative of a linear ML model, intending to compare the resulting performances with the aforementioned non-linear ML routines;

We assessed the performances of the ML models mentioned before through a 5-fold Cross Validation, in terms of Sensitivity (S), Specificity ($SP$) and Accuracy (A), widely employed in medicine.

One of the purposes of this work was the search for the Raman bands associated with biomolecules (collagen components, DNA, lipids, hydroxyapatite, etc.). These findings could allow to understand the biochemical processes behind the occurrence of such diseases and, on the other hand, to reduce the amplitude of the spectral region employed for the acquisition. In the framework of ML and DL, this problem is referred to as the assessment of the Feature Importance (FI). In turn, we adopted the so-called Permutation Feature Importance (PFI) [25,26,27]. Consider an ML or DL model M, previously trained on the dataset represented by the matrix X. Let s be a generic score, assessing the performance of M on a test dataset $X^{\prime }$. Suppose to calculate the parameter $s_j$, representing the aforementioned score, on the dataset $X^{\prime }$ after a permutation of the elements of the j-th column. The basic idea behind the definition of PFI is that, if $s_j\sim s$, the spectral component corresponding to the j-th column of $X^{\prime }$ can be considered irrelevant for the classification. The PFI associated to the j-th spectral component is defined as

$$PF{I}_j=s-\frac{1}{K}\sum _l^K{s}_{l,j}.$$

(3)

in this definition, $X^{\prime }$ is subjected to K permutations of the j-th column, and $s_{l,j}$ represents the score obtained at the l-th permutation. By definition, the larger $PF{I}_j$, the larger the importance of the j-th spectral component. Despite RFC provides an intrinsic definition of FI, based on the ability of a single feature to increase the “purity” of the classification domains [28], the use of PFI allowed to compare the FI of different ML models and to reinforce our considerations about the relevant chemical compounds in the biochemical processes at the origin of the diseases under interest. Hence, we calculated the PFI associated with the prediction accuracy ($s\equiv A$), on a $X^{\prime }$ dataset obtained by randomly selecting spectra from X, with a proportion of 20% with respect to the total number of spectra. We chose this proportion in line with the 5-fold cross-validation adopted. In addition, we chose $K=5$.

All the calculations were carried out through a hybrid script, obtained by customizing an Orange Data Mining script with Python codes [26].

3. Results and Discussion

In Figure 3 (a), we showed the averaged spectra associated with CS and EC. The shaded areas represent the related standard deviation. The corresponding peak assignments are reported in Table 1. The general behavior of the averaged spectra highlighted a decrease in the Raman signal of CS with respect to EC. In particular, we qualitatively observed the largest differences between CS and EC in the Raman intensities in the following spectral bands: the most intense peak ($\sim 1003{\mathrm{cm}}^{-1}$) is attributed to Phenylalanine (Phe), recognized as a strongly Raman-active molecule [29]. Phe is a precursor of several amino acids, such as Tyrosine, represented by the peaks at 815, 1172 and $1207{\mathrm{cm}}^{-1}$. Tyrosine was previously recognized as an amino acid regulating the production and activity of osteoclasts [30] The presence of Phe is testified also by the peak at $1609{\mathrm{cm}}^{-1}$. We attributed the two narrow peaks at $\sim 729{\mathrm{cm}}^{-1}$ to carbonates that, alongside bioapatites (604, 849, 1035, 1057 and $1098{\mathrm{cm}}^{-1}$), represent one of the most abundant inorganic components of the bone tissue. Other relevant spectral components can be found in the peaks at 830 and $1453{\mathrm{cm}}^{-1}$, assigned to proline (collagen) and ${\mathrm{CH}}_2$ wagging, respectively.

The small standard deviation in comparison to the average Raman spectra observed for the CS class in Figure 3 (a) resulted in small differences between the averaged Raman spectra of G1, G2 and G3, as shown in Figure 3 (b). However, the peak of Phe qualitatively highlighted the largest differences between the three grades of CS. This last feature makes the peak of Phe a promising candidate for the grading of CS.

In view of the employment of the collected Raman spectra to classify the analyzed tissues, we emphasize that the apparent macroscopic difference in the signal intensity of CS and EC observed in Figure 3 (a) is not necessarily sufficient to perform an accurate classification based on the sole visual interpretation of the Raman spectra. We stress that the shaded areas associated with the mean spectra represent the standard deviation. For this reason, CS and EC could have a superposition region in the intensity domain, due to Raman spectra corresponding to the tails of the related intensity distributions, possible source of misclassification problems. In this case, the employment of ML algorithms allowed to rely on additional information, provided by the spectral components, in a systematic and efficient fashion.

In Table 2, we resumed the performances of the ML and DL models for the examined classification problems. In particular, we averaged the Accuracy A, the Sensitivity S and the Specificity $SP$ over the 5 folds and over the values of the label. The most evident result highlighted is the substantial difference between the performances of the linear SVM with respect to the other models, for all the classification problems considered. This feature revealed how non-linearities in the definition of the models resulted in enhanced classification performances, with relative uncertainties of $\sim {10}^{-3}\%$, indicating a negligible overfitting. As expected, the simplest binary classification problem EC-CS led to the maximum performances, i. e. $A=(99.7\pm 0.1)\%$, $S=(99.7\pm 0.1)\%$ and $SP=(99.0\pm 0.1)\%$, reached with MLPC(ADAM). This result is in line with previous works about the Raman-based detection of CS and EC [6,18,19]. On the other hand, the classification problem EC-G1-G2-G3 resulted in the minimum classification accuracy, with the maximum value of $A=(97.6\pm 0.1)\%$ obtained with MLPC(ADAM). In conclusion, the tested routines showed excellent classification performances, making the technique a promising candidate for future applications in the diagnosis and grading of CS.

Figure 4 reports the PFI associated with the prediction accuracy for the models examined as a function of wavenumber, for all the classification problems of interest. In particular, we normalized PFI through the MinMax rule [59], to highlight the spectral components corresponding to the maximum PFI. We emphasize that the PFI provides a measure of the relevance of a spectral component in detecting the values of the label. In this sense, the behavior of PFI for MLPC (Figure 4 (a), (c) and (e)), characterized by the presence of a large number of peaks of similar intensities, indicated that the DL models studied rely on a large number of spectral components to reach the aforementioned classification performances. This trend is probably the result of the strongly non-linear nature of MLPC. On the other hand, the PFI of SVM and RFC, represented in Figure 4 (b), (d) and (f), showed a small number of well-definite peaks in comparison to MLPC. Therefore, despite SVM and RFC showed worse classification performances than the DL models, they are potentially capable of providing more information about the biochemical mechanisms on the basis of the malignant degree under interest. In particular, by observing the Figure 4 (b), we detected a peak of PFI at $\sim 602{\mathrm{cm}}^{-1}$, exhibiting large feature importance either for SVM or for RFC. As shown in Table 1, this peak is attributable to phosphate groups (P-O bending) associated with hydroxyapatite. This result is in line with previous investigations, attributing the growing presence of bioapatites to the production of calcifications [60]. To conclude the analysis of PFI for the G1-G2-G3 classification problem, we mention also the peaks related to SVM and located at $\sim 730$ and $\sim 1449{\mathrm{cm}}^{-1}$. These spectral components can be attributed to DNA and ${\mathrm{CH}}_2$ bending, respectively. These findings suggested differences in EC and CS attributable to cell proliferation. Analogously, the peak of Phe ($\sim 1003{\mathrm{cm}}^{-1}$) observed either for SVM or for RFC in Figure 4 (d) suggested relevant differences between different grades of CS in the content of this molecule. Finally, as expected, in the classification problem EC-G1-G2-G3, the behavior of the PFI (Figure 4 (f)) appeared to be intermediate between the PFI in the EC-CS and the G1-G2-G3 classification problems.

4. Conclusions and Future Perspectives

Chondrosarcoma and the classification of its various degrees of malignancy has recently become a case study for the application of Raman spectroscopy in oncology. In this investigation, we tested the ability of different ML and DL routines to accurately discriminate between CS and EC from a Raman-based experimental dataset, on which several multivariate statistical approaches have already achieved good classification performances [6], or DL algorithms based on wavelet transform [18], as well as topology machine learning [19].

The classifiers used in this paper allowed to achieve two main objectives. First of all, we observed how the MLPC(ADAM) classifier was able to distinguish CS and EC with an accuracy reaching 99.7% and to distinguish the specific grades of CS with an accuracy reaching 99.2%. In addition, the aforementioned DL algorithm was able to distinguish the three grades of CS and EC within a single-stage process, with an accuracy reaching 97.6%. This outcome highlighted the technique as a promising route for future applications in the intra-operative diagnosis of bone cancer, to correctly estimate the tumor margins. Furthermore, the analysis of PFI highlighted that, despite the worse classification performances of RFC and SVM in comparison to MLPC, they offer the capability to detect characteristic tissue components (proline, DNA, bioapatites), particularly relevant for the grading. In particular, according to this analysis, Raman bands associated with bioapatites turned out to be the most relevant in the distinction between CS and EC. On the other hand, Phenylalanine (Phe) turned out to be the most relevant amino acid in grading CS.

The main limitation associated with our analysis is represented by the employment of a relatively small statistical sample, either in terms of the number of patients or the number of collected Raman spectra. Therefore, the results generated by our analysis have to be considered an exploratory investigation, representing the basis for further analysis relying on more statistically significant and larger RS-based datasets. Furthermore, as already described in the Introduction section, the richness of information within a single Raman spectrum complicates the data interpretation, requiring the employment of ML or DL to retrieve the information of interest. This operation often requires high computational cost and/or long computational times. This issue may represent a problem, especially for the implementation of such technology on engineered probes for real-time intra-operative diagnosis, where diagnostic times of seconds or fractions of seconds are strongly required.

To achieve such challenging objectives, further efforts have to be devoted, distributed on multiple routes: first of all, to test ML or DL routines on specific Raman bands, such as the Raman bands corresponding to the peaks of PFI, to assess the ratio between the computational time needed to obtain the response and the classification performances; in turn, to realize cost-effective and powerful hardware, able to speed up the calculation; to operate on the software, by implementing the classification routines in high-performance computing codes.