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 techniques to other approaches, such as MRI [4].

In the last two decades, RS demonstrated to be 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 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 known as 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 improved 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, managing the dataset originally collected in [6]. Our analysis was focused either on the basic problem of distinguishing CS from EC (or on assigning a correct grade to the CS cases), either on the crucial question to identify the most relevant Raman bands of any malignat degree. About the first question, the results demonstrated the effectiveness of ML and DL models in identifying the aforementioned CS tumor, with classification accuracy reaching 99.7%. But much more relevant, the method for the identification of the RBI, i.e., the Raman bands evidencing the degree of malignancy of tissues under investigation, basically inspired by Permutation Feature Importance [20], allowed to identify Raman spectral components, more relevant to the classification process. Such results represent the most original contribution of our paper. This is because the analysis of RBI 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 μ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 [21] 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× objective and the out-of-focus components of the detected light were suppressed through an aperture (pinhole) of diameter 25 μ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 ~2 cm⁻¹ and they were defined within the spectral interval between 400 and 1800 cm⁻¹. In Figure 2, we showed images of H&E stained sections associated with the four types of tissue examined in this work.

Figure 1. Schematic representation of a Confocal Raman Microscope, equipped with a Rayleigh rejection filter, a pinhole, a high-resolution diffraction grating, and an electron multiplied CDD camera as a detector.

Figure 1. Schematic representation of a Confocal Raman Microscope, equipped with a Rayleigh rejection filter, a pinhole, a high-resolution diffraction grating, and an electron multiplied CDD camera as a detector.

Figure 2. (a) image of a stained histological section of tissue corresponding to G1; (b) image of a stained histological section of tissue corresponding to G2; (c) image of a stained histological section of tissue corresponding to G3; (d) image of a stained histological section of a tissue corresponding to EC. The images of the histological sections were obtained with a magnification of 20×.

Figure 2. (a) image of a stained histological section of tissue corresponding to G1; (b) image of a stained histological section of tissue corresponding to G2; (c) image of a stained histological section of tissue corresponding to G3; (d) image of a stained histological section of a tissue corresponding to EC. The images of the histological sections were obtained with a magnification of 20×.

The spectra were collected in grids of pixel resolution 4 μm × 4 μ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 a single Raman spectrum as a vector x of dimension $N_p$. In the supervised learning framework, a discrete variable y, called label, assuming, is associated with x. As stated before, we will have a set $\left\{(x_i,y_i)\right\}$, with $i=1,...,N_s$, representing all the available instances. The instances $\left\{x_i\right\}$ can be conveniently conceived as the rows 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). In this case, y can assume the values EC and CS;
• The problem of distinguishing G1, G2 and G3 (G1-G2-G3). In this case, y can assume the values G1, G2 and G3;
• The problem of distinguishing EC, G1, G2 and G3 (EC-G1-G2-G3). In this case, y can assume the values EC, G1, G2 and G3.

To this aim, we employed three ML protocols:

• Support Vector Machine (SVM): this ML was originally conceived to find the hyperplane that correctly separates the clouds in the space of dimension $N_p$ corresponding to the different values of y. Let $w·x+b=0$ the equation of the hyperplane, where w is an $N_p$-dimensional vector normal to the hyperplane, b a real number, and · denotes the scalar product. In the version of SVM implemented in this work, the parameters related to w and b were optimized by solving the following bound minimum problem:

  $$\left\{\begin{array}{c}\underset{w,b,{ϵ}_i}{min}{\left|\right|w\left|\right|}^2+C\sum _{i=1}^{N_s}{ϵ}_i\hfill \\ y_i(w·x_i+b)\ge 1-{ϵ}_i,\hfill & i=1,...,N_s\hfill \\ {ϵ}_i\ge 0\hfill \end{array}\right.$$

  (1)
  In this equation $\left|\right|w\left|\right|$ denotes the Euclidean norm of w. In addition, C, here fixed to 0.4, represents a regularization parameter accounting for the penalty related to the presence of outliers, i.e. instances. Finally, the parameters $\left\{{ϵ}_i\right\}$, here fixed to 0.8, are variables aimed at relaxing the constraints [22]. Despite more advanced versions of SVM, conceived to separate the classes of interest through hyper curves more complex than a hyperplane, we adopted the classical linear kernel as representative of a linear model and to compare it with more complex ML classifiers;
• Random Forest Classifier (RFC): RFC is a “bagging” ML routine, building the classifier as an ensemble, or forest, of decisional trees. A single tree is built from an ensemble of Raman spectra obtained through a procedure called bootstrap sampling. This selection procedure is aimed at selecting $N_s$ spectra, with the possibility of repeated spectra within the resulting dataset. Then, k spectral components, with $1<k\le N_p$ are randomly selected to build each node of a single tree. Among the previously selected spectral components, the best one is selected according to its capability to discriminate the values of the label y. In this paper, we adopted Gini’s index as the parameter quantifying this ability. The growth of each tree was stopped when a perfect separation of the spectra belonging to the training dataset was obtained. The final RFC classifier contained 4000 trees, considered a good solution to obtain small levels of overfitting and small out-of-bag error [23]. The response of the entire classifier is thus given by the majority of the prediction of its trees.
• Multi-Layer Perceptron (MLPC): in this investigation, this simple DL algorithm consisted of a single hidden layer with 400 neurons and zero bias. The hidden layer is followed by an activation layer employing a ReLU activation function. We carried out the training either with the ADAM [24] or with the L-BFGS-B [25] 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;

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. IN particular, in the following, we will refer to S, $SP$, and A as quantities averaged over the values of the label y and over the folds of the Cross-Validation.

One of the purposes of this work was the search for the Raman components associated with molecules relevant to the classification process. By definition, these molecules can be considered biomarkers of the diseases under investigation. In particular, the search of cancer biomarkers could provide interesting information about the biochemical processes behind the occurrence of such diseases. Furthermore, this analysis could help to reduce the amplitude of the spectral region needed to reach a reliable diagnosis. This perspective could have important consequences in terms of the future implementation of this technology in simple and low-cost engineered devices. To this aim, we adopted the so-called Raman Band Identification (RBI), inspired by the so-called Permutation Feature Importance (PFI) [20,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 RBI 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 RBI associated to the j-th spectral component is defined as

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

(2)

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 $RB{I}_j$, the larger the importance of the j-th spectral component. Despite RFC providing an intrinsic definition of Feature Importance, based on the ability of a single feature to increase the “purity” of the classification domains [28], the use of RBI allowed to compare the results of different ML models directly 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 RBI 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 3a, 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 (~1003 cm⁻¹) 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 cm⁻¹. 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 cm⁻¹. We attributed the two narrow peaks at ~729 cm⁻¹ to carbonates that, alongside bioapatites (604, 849, 1035, 1057 and 1098 cm⁻¹), 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 cm⁻¹, 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 3a resulted in small differences between the averaged Raman spectra of G1, G2 and G3, as shown in Figure 3b. 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 3a 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 RBI 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 RBI through the MinMax rule [59], to highlight the spectral components corresponding to the maximum RBI. We emphasize that the RBI provides a measure of the relevance of a spectral component in detecting the values of the label. In this sense, the behavior of RBI for MLPC (Figure 4a,c,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 RBI of SVM and RFC, represented in Figure 4b,d,f, showed a small number of well-definite peaks in comparison to MLPC. Therefore, despite SVM and RFC showing 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 4b, we detected a peak of RBI at ~602 cm⁻¹, 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 RBI for the G1-G2-G3 classification problem, we mention also the peaks related to SVM and located at ~ 730 and ~1449 cm⁻¹. 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 (~1003 cm⁻¹) observed either for SVM or for RFC in Figure 4d 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 4f) 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, as well as the identification of the most relevant Raman bands associated to any specific malignat degree.

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%. Hence, 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, in comparison to MLPC, SVM and RFC achieved worse classification performances. However, the analysis of RBI revealed that SVM and RFC rely on a smaller number of spectral components in comparison to MLPC to reach the aforementioned classification performances. This feature results in two important consequences: first of all, the spectral components maximizing the RBI of SVM and RFC potentially indicated the presence of biomarkers (proline, DNA, bioapatites) for the diseases under investigation. 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. Secondly, these findings suggest the possibility to restrict the spectral interval on the bands maximizing RBI, with negligible loss of the classification performances. In view of future implementations of this technology in engineered devices, the perspective of restricting the spectral interval needed for accurate and robust diagnosis undoubtedly represents an advantage, either in terms of acquisition speed or of cost of the hardware.

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. In general, the employment of datasets where the number of features of each instance is similar or larger than the number of instances available is traditionally associated with the occurrence of overfitting. Despite in this work we observed a negligible level of overfitting, further investigations have to be carried out in order to further monitor the robustness of the models. For this reason, the results generated by our analysis have to be considered an exploratory investigation, representing the basis for further studies relying on more statistically significant and larger RS-based datasets.

In general, further efforts have to be devoted, distributed on multiple routes: primarily, to test ML or DL routines on specific Raman bands, such as the Raman bands corresponding to the peaks of RBI, to assess the ratio between the computational time needed to obtain the response and the classification performances; therefore, 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.