2.1.

Patient Samples and Imaging Method

The dataset consisted of 483 PC patients from three different institutions: Centre hospitalier de l’Université de Montréal (CHUM), University Health Network (UHN), and Centre hospitalier universitaire de Québec-Université Laval (CHUQc-UL). The characteristics of the dataset are shown in Table 1. Tumor stages defines where cancer is present in prostate tissue; there are three tumors stages: pT2 (organ confined), pT3a (extraprostatic extension or microscopic invasion of bladder neck), and pT3b (seminal vesicle muscle wall invasion). The Gleason score details the arrangement and pattern of cancer cells. There are five patterns, and the Gleason score is equal to the sum of the most and second most common pattern present prostate tissues. Confocal $\mathrm{R}\mu \mathrm{S}$ measurements on formalin-fixed paraffin embedded tissue microarrays (TMAs) of PC and IDC-P tissue were acquired as described in our previous work.² Briefly, sections from the TMAs were transferred onto aluminum slides with low Raman activity (Miro5011, Anomet, Brampton, Ontario) and dewaxed according to the CHUM standard clinical dewaxing protocol. Raman acquisitions were made using a confocal Raman microscope (Renishaw, Gloucestershire) equipped with a 785-nm line focus laser and a grating of 1200 lines/mm allowing acquisition of Raman shifts between 602 and $1726{\mathrm{cm}}^{-1}$. For each TMA, five accumulations of 10 s each at a 150-mW laser output power were acquired using a $50\times$ objective with a 0.75 numerical aperture. There were always four spectra taken per TMA core at four different locations, whereas the number of TMA cores per patient varied between institutions. The diameter of TMA cores ranged between 0.6 to 1.2 mm and the probed region corresponded to a rectangular area of $24\mu {\mathrm{m}}^2$² ($8\mu \mathrm{m}\times 3\mu \mathrm{m}$, approximately corresponding to single-cell analysis). IDC-P identification was done by two independent pathologists. All Raman spectra files are available in the Dryad Digital Repository database.¹⁶

Table 1

Patient clinicopathological characteristics separated by institution: CHUM, Centre hospitalier de l’Université de Montréal; CHUQc-UL, Centre hospitalier universitaire de Québec–Université Laval; IDC-P, intraductal carcinoma of the prostate; IQR, interquartile range; PSA, prostate-specific antigen; TMA, tissue microarray; UHN, University Health Network.

2.2.

Data Preprocessing

The following data preprocessing procedures are applied to obtain each spectrum: (1) summed the five accumulations of 10 s each, (2) remove cosmic rays, (3) remove background signal produced by the aluminum slides and tissue fluorescence using the rolling ball algorithm,¹⁷ (4) apply standard normal variate (SNV) normalization, (5) average the four spectra of each TMA core resulting in between 1 and 3 spectra per patient. During background removal, the rolling ball algorithm distorts, due to border effects, Raman signals by creating peaks that are not present in the raw spectra. To avoid these artifacts, regions between 602 and $669{\mathrm{cm}}^{-1}$, as well as between 1707 and $1726{\mathrm{cm}}^{-1}$ were removed.

A limitation of the rolling ball algorithm comes from the presence of a finite size structuring element within the algorithm. At the beginning and end of a spectrum, half of this structuring element extends past the range of spectra and causes distortions in the extracted Raman signals due to this border effect. The magnitude of this effect depends on the smoothness of the signal baseline. In this particular dataset, the baseline shows a step non-smooth gradient at the beginning of the signal, which increases the distortion, and a smooth baseline at higher wavelengths. It was determined, empirically, that to remove any residue due to this effect, signal regions between 602 and $669{\mathrm{cm}}^{-1}$ as well as between 1707 and $1726{\mathrm{cm}}^{-1}$ had to be removed, i.e., a portion of the beginning and end of each spectrum.

The next preprocessing step is to fit Gaussian functions over the most common Raman peaks present in the dataset. Since the Gaussian function is strictly positive, the minimum value of SNV normalized spectra was subtracted from the entire dataset to obtain only positive values. The peak fitting algorithm, which is an improvement of the algorithm first used in Ref. 15, considers a peak only when it appears in at least 50% of spectra at the same location or within $\pm 2{\mathrm{cm}}^{-1}$ and if its height is above a threshold height previously determined empirically. For every peak, both inflection points are determined as well as the height of the maximum and its position. The distance between inflection points is used to compute an approximate standard deviation ($\sigma$) and to determine the fitting domain for a nonlinear least-squares regression algorithm.¹⁸ This standard deviation, and the maximum and position of the peak determine the three starting points given to the regression algorithm to fit a Gaussian over the wavelength domain of a peak to obtain fitted heights and widths (standard deviation, $\sigma$). These values are then SNV normalized and used as features.

2.3.

Statistical Analysis

A machine learning model with several hyperparameters was trained on the CHUM dataset and tested on two independent hospital cohorts: CHUQc-UL and UHN. Patient clinicopathological characteristics and number of spectra of each classification task are shown in Table 1. The model uses two types of features: individual features; intensity set, and peak features; peak set. To determine the features in each set, two separate feature selection algorithms are used. The selected features are then combined in a single feature vector and passed to the RBF kernel SVM classifier. Before training a model on the CHUM dataset, cross validation (CV) is required to determine the optimal value of the hyperparameters. For each classification task, the CV identifies the set of hyperparameters, through a grid search, that maximizes the area under the curve (AUC) of the receiver-operator-characteristic (ROC) curve of the CHUM dataset for three independent models: one model using only the intensity features, one using only the peak features, and a combined model that used both feature sets; $\text{intensity}+\text{peak}$. Figure 1 shows a schematic workflow of the statistical analysis.

Fig. 1

Data acquisition and machine learning workflow. Annotated stained TMAs are used to guide Raman microspectroscopy measurements. Data are split into training (CHUM) and testing sets (UHN & CHUQc-UL). A double nested five-fold CV is performed to determine the hyperparameters through a randomized grid search. Feature selections on peak and intensity features are done separately using linear kernel L1 SVM. Classification models are RBF kernel SVMs that use selected features. The hyperparameter set that yields the highest AUC on training data is used to build the final model (green), which is evaluated on predicting outcome on the testing set (orange).

2.4.

Feature Engineering and Feature Selection

Two different sets of SVN normalized features were used for classification. The first set, intensity, consists of individual intensity values of single Raman bands. To extract the most relevant intensities, an L1 regularization linear SVM with one hyperparameter (regularization parameter $C_{\text{intensity}}$) is utilized for feature ranking. The linear SVM assigns coefficients to individual intensity features, which are then used to rank them by importance. Features that are close to one another, which are most likely correlated and contain the same information, are removed using the hyperparameter $n_{\text{neighbor}}$. Every feature that is $\pm n_{\text{neighbor}}{\mathrm{cm}}^{-1}$ from the feature with the highest coefficient is removed. The feature set is then updated, and the same procedure is applied for the feature with the second-highest coefficient and so forth. For example, if the intensity features with the highest coefficient has a Raman shift of $780{\mathrm{cm}}^{-1}$ and $n_{\text{neighbor}}=5$, intensity features with a Raman shift between 775 and $785{\mathrm{cm}}^{-1}$ will be removed. This procedure also reduces the dimensionality of the feature set.

The second set, peak, consists of Raman peaks fitted with a Gaussian distribution. The height and the width of these Gaussian fitted peaks are used as features. A linear SVM with L1 regularization (regularization parameter $C_{\text{peak}}$) is used to rank these engineered features by importance.

To further reduce the number of features within each feature set, the $k_{\text{intensity}}$ and $k_{\text{peak}}$ features with the highest ranked coefficient are used for modeling and the others are discarded. Both $k_{\text{intensity}}$ and $k_{\text{peak}}$ are hyperparameters determined during CV.

2.5.

Modeling and Cross Validation

An RBF kernel SVM was used as a classifier. This algorithm has two hyperparameters: $C_{\mathrm{RBF}}$ and $\gamma$. A grid search was performed using CV to optimize the following hyperparameters: $C_{\mathrm{RBF}}$, $C_{\text{intensity}}$, $C_{\text{peak}}$, $\gamma$, $k_{\text{intensity}}$, $k_{\text{peak}}$, and $n_{\text{neighbor}}$. Both $k_{\text{intensity}}$ and $k_{\text{peak}}$ ranged from 5 to 100, $C_{\mathrm{RBF}}$ ranged from 0.001 to 1, $\gamma$ ranged from 0.001 to 0.01, $n_{\text{neighbor}}$ from 0 to 10, and $C_{\text{intensity}}$ and $C_{\text{peak}}$ from $1\times {10}^{-6}$ to $1\times {10}^{-4}$. The CV scheme performed was a five repeat five-fold CV. The combination of hyperparameters generating the model with the highest AUC value was chosen to train the model using the complete training set. The selected model was then applied to the two independent testing sets.

3.1.

Benign/Cancer Classification

The selected hyperparameters, optimized by cross-validated grid-search, for the intensity, peak, and $\text{intensity}+\text{peak}$ models are : $k_{\text{intensity}}=(15,NA,16)$, $C_{\text{intensity}}^{L1}=(1\times {10}^{-4},NA,1\times {10}^{-4})$, $n_{\text{neighbor}}=(5,NA,6)$, $k_{\text{peak}}=(NA,\mathrm{6,8})$ with $C_{\text{peak}}^{L1}=(NA,1\times {10}^{-6},5\times {10}^{-5})$, $C^{\mathrm{RBF}}=(0.07,0.2,0.06)$, and ${\gamma }^{\mathrm{RBF}}=(5\times {10}^{-3},5\times {10}^{-4},5\times {10}^{-3})$, where $NA$ stands for Not Applicable, e.g., the intensity model does not have a $k_{\text{peak}}$ hyperparameter and is thus NA. Combining the two feature types improved the classification metrics for the UHN and CHUQc-UL testing sets and, in both cases, the intensity model was superior to the peak model. The difference in classification performance [accuracy, sensitivity, specificity, AUC] between the intensity and $\text{intensity}+\text{peak}$ feature sets for both training sets was: $+1\%$, $+0\%$, $+4\%$, $+2$. Performance metrics obtained using the combined model compared to our previous study changed according to the following: $+0\%$ ($-5\%$), $+0\%$ ($-7\%$), $+2\%$ ($+1\%$), $+4$ ($+6$). While some classification metric decreased, the AUC increased for both testing sets. The results for the training set are similar to the results obtained on the testing sets, i.e., a slight decreased in accuracy ($-2\%$), sensitivity ($-1\%$), specificity ($-4\%$), and an increased in AUC ($+3$) while retaining similar uncertainties.

The average Raman spectra for both benign and cancer prostate tissue and the five highest ranked features of each feature set are shown in Fig. 2. The description of the five best features for each the peak and the intensity models are presented in Table 3 along with their respective molecular assignments.²³ For the peak height feature set, four Raman peaks ranging from 719 to $1032{\mathrm{cm}}^{-1}$ were increased in the average spectrum of benign tissue. DNA, RNA, proteins, and amino acids (i.e., tryptophan, proline, valine, and phenylalanine) were the main biochemical components associated with these peaks. The only feature associated with an increase in cancer tissue was the peak at $1294{\mathrm{cm}}^{-1}$, which was assigned to lipids. For the intensity feature set, a unique peak ($1250{\mathrm{cm}}^{-1}$) associated with the $\beta$-sheet secondary structure of proteins was more intense in the Raman spectrum of cancer tissue. There were also 4 Raman peaks increased in benign tissue, from 719 to $1003{\mathrm{cm}}^{-1}$ and their molecular assignments are mostly DNA, RNA, proteins, and amino acids (i.e., tyrosine, proline, and phenylalanine).

Fig. 2

Average Raman spectra of cancer (red, $n=272$) and benign (black, $n=99$) tissue annotated with the five highest ranked individual intensity features (purple) and the five highest ranked peak features (teal, for height features).

Table 3

Highest ranked features used for the classification of benign and malignant prostate tissue found using a linear SVM with L1 regularization, and their associated Raman peaks. Tentative molecular assignment of prostate Raman peaks based on the literature. The Raman shift of the features (feature column) and of its associated Raman peak (peak, center column) are reported.2,9–12,20–22

3.2.

Cancer/IDC-P Classification

The optimal hyperparameters to classify cancer from IDC-P tissue using the (intensity, peak, and $\text{intensity}+\text{peak}$) models are: $k_{\text{intensity}}=(77,NA,76)$ with $C_{\text{intensity}}^{L1}=(1\times {10}^{-3},NA,5\times {10}^{-6})$, $n_{\text{neighbor}}=(7,NA,8)$, $k_{\text{peak}}=(NA,74,90)$ with $C_{\text{peak}}^{L1}=(NA,1\times {10}^{-2},1\times {10}^{-2})$, $C^{\mathrm{RBF}}=(2,2,5)$, and ${\gamma }_{\text{peak}}^{\mathrm{RBF}}=(0.01,0.01,0.01)$ where NA stands for Not Applicable. The performance metrics of the combined model for the UHN (CHUQc-UL) testing set improved our previous results by: $+4\%$ ($+8\%$), $+7\%$ ($+9\%$), $+2\%$ ($+6\%$), $+9$ ($+9$). It also improved the results on the CHUM training set by: 2%, 1%, 3%, 4.7 while also decreasing the standard deviation by: $-1\%$, $-3\%$, $-1\%$, 2.7.

The description of the five highest ranked features for each the peak and the intensity models are shown in Fig. 3 and Table 4 along with their respective molecular assignments. Within the five highest ranked peak features, three Raman peaks were increased in the average spectrum of IDC-P compared to cancer tissue. These peaks were mostly associated with amino acids: $825{\mathrm{cm}}^{-1}$ (tyrosine), $1170{\mathrm{cm}}^{-1}$ (tyrosine), and $932{\mathrm{cm}}^{-1}$ (proline and valine). Raman peaks assigned to phenylalanine ($1001{\mathrm{cm}}^{-1}$), protein, lipids, DNA, and RNA ($1666{\mathrm{cm}}^{-1}$) were reduced in IDC-P. For the intensity feature set, the $1484{\mathrm{cm}}^{-1}$ Raman peak, assigned to the nucleobases adenine and guanine, was the single decreased peak in the average spectrum of IDC-P. The four other Raman peaks, ranging from 944 to $1266{\mathrm{cm}}^{-1}$, were identified mostly in IDC-P. The biochemical constituents associated with these peaks were predominately linked to proteins: $\alpha$-helix and $\beta$-sheet secondary structures, proline, valine, and phenylalanine.

Fig. 3

Average Raman spectra of cancer (red, $n=272$) and IDC-P (black, $n=112$) tissue along with the 5 highest ranked individual intensity features (purple) and the five highest ranked peak features (teal showing height features; blue showing width features).

Table 4

Highest ranked features used to classify intraductal carcinoma of the prostate and invasive prostate cancer tissue found using a linear SVM with L1 regularization, and their associated Raman peaks. Tentative molecular assignment of prostate Raman peaks based on the literature. The Raman shift of the features (feature column) and of its associated Raman peak (peak, center column) are reported.2,9–12,20–22

3.3.

HGPIN/IDC-P Classification

By using only the individual intensity features a classification with 100% AUC was obtained for both testing sets, and combining both feature sets also yielded the same result. The optimal hyperparameters for all three data sets were $k_{\text{peak}}=42$ with $C_{\text{peak}}^{L1}=1\times {10}^{-4}$, $k_{\text{intensity}}=49$ with $C_{\text{intensity}}^{L1}=1\times {10}^{-4}$, $C^{\mathrm{RBF}}=3$, ${\gamma }_{\text{peak}}^{\mathrm{RBF}}=0.003$, and $n_{\text{neighbor}}=4$. The improvements in accuracy, sensitivity, specificity, and AUC with respect to our previous study were: $+2.2\%$ ($+1.7\%$), $+4.5\%$ ($+3.6\%$), $+0\%$ ($+0\%$), $+2.3$ ($+0$) for the UHN (CHUQc-UL) testing set. Those results were obtained while also improving the result on the training set (CHUM) by: $+2.5\%$, $+0.8\%$, $+2.9\%$, and $+2.4\%$.

Following the classification of IDC-P and HGPIN, the five highest ranked features were assigned to specific molecules (Fig. 4 and Table 5). Only one feature was found to be a dominant contributor for identifying HGPIN using the peak height feature set. More specifically, the $937{\mathrm{cm}}^{-1}$ Raman peak was identified as $\alpha$-helix secondary structure of proteins and two amino acids, proline and valine. The other four features, from 1001 to $1238{\mathrm{cm}}^{-1}$, were all increased in the average spectrum of IDC-P tissue and were assigned to the $\beta$-sheet secondary structure of proteins, phenylalanine, tryptophan, and tyrosine. For the classification using the intensity feature set, two Raman peaks were increased in HGPIN tissue: phenylalanine ($1032{\mathrm{cm}}^{-1}$); $\alpha$-helix secondary structure of proteins, fatty acid, and the nucleobase thymine ($1667{\mathrm{cm}}^{-1}$). Raman peaks at $1003{\mathrm{cm}}^{-1}$ (phenylalanine) and at $1250{\mathrm{cm}}^{-1}$ ($\beta$-sheet) were mostly found in IDC-P tissue. The carotenoid biochemical component ($1152{\mathrm{cm}}^{-1}$) was a main contributor to the IDC-P classification in both feature sets.

Fig. 4

Average Raman spectra of HGPIN (red, $n=170$) and IDC-P (black, $n=112$) tissue along with the five highest ranked individual intensity features (purple) and the five highest ranked peak features (teal, for height features).

Table 5

Highest ranked features used to classify high-grade prostatic intraepithelial neoplasia and IDC-P found using a linear SVM with L1 regularization, and their associated Raman peaks. Tentative molecular assignment of prostate Raman peaks based on the literature. The Raman shift of the features (feature column) and of its associated Raman peak (peak, center column) are reported.2,9–12,20–22