2.1.

Raman Spectroscopy System

Measurements were made using an intra-operative hand-held single-point Raman spectroscopy probe system manufactured by the company Reveal Surgical (Montreal, Canada). The fiber-optics probe design of similar systems was presented elsewhere.^30,36 Briefly, it integrates ten fiber-optics of $100\text{-}\mu \mathrm{m}$ diameter core. Of these fibers, nine were used for light detection and a central fiber was used for tissue excitation. A band-pass filter covered the tip of the excitation fiber, allowing illumination centered around 785 nm, and a high-pass filter was used for the collection fibers. The probe was connected to a 785-nm laser and a spectrometer with high sensitivity at wavelengths ranging from 800 to 900 nm, with an average resolution of $1.8{\mathrm{cm}}^{-1}$ across the spectral detection domain. A converging lens at the tip of the probe ensured contact measurements interrogated a spot size of diameter $\sim 0.5\mathrm{mm}$. The system was controlled by a proprietary software (Reveal Surgical, Montreal, Canada) that allowed acquisition parameters to be set by the user, including laser power, exposure time per spectrum, and number of repeated measurements (i.e., accumulations) at each point.

2.2.

Patients Selection

About 20 patients were recruited for this study. Of those, 19 underwent breast surgery (lumpectomy or mastectomy) following a diagnosis of invasive breast cancer (lobular or ductal) and one patient underwent breast reduction surgery. For that patient, spectroscopic measurements were made from a breast in which no tumor was radiologically detectable. However, the patient was diagnosed with breast cancer associated with a tumor detected in the contralateral breast. All patients recruited in the study were undergoing breast surgery for the first time, did not have neoadjuvant therapy, had a cancer grade inferior to 4, and had a tumor larger than 1 cm.

Patient specimens were utilized to build an ex vivo dataset of Raman spectroscopy measurements combined with histopathological and clinical data. Informed consent was obtained before the patient underwent surgery (McGill University Health Center Ethics Committees, approval number 2021-5997). Clinical data available included age, tumor type, tumor size, and Nottingham histologic score. Patient demographic details are provided in Table 1 (see also Table 3 for demographic details on each patients). The breast reduction surgery patient was recruited with the intention to acquire more measurements in healthy tissue to ensure a more balanced dataset.

Table 1

Clinical and pathological characteristics of all patients undergoing breast surgery (std: standard deviation; NA: not available; and IQR: interquartile range).

2.3.

Specimen Handling and ex vivo Spectroscopic Measurements

Twenty breast specimens (one per patient) were extracted by the surgeons (either S.M. or F.T.). The specimens were weighed, measured, and marked with ink, following institutional standards. They were then cut in serial slices of $\sim 5\text{-}\mathrm{mm}$ thickness. For each specimen, two smaller samples were cut from two different slices. Those samples were selected by the pathologist (A.O.) based on visual inspection with the objective of having one sample containing cancer cells and one sample mostly composed of healthy tissue. The samples were fixed on a cardboard support using two pins and a photograph was taken. The support was placed on a grid template and the Raman spectroscopy probe was fixed above the center of the grid [Fig. 2(a)]. Before each measurement, the cardboard support was moved to a different location on the template. The $xy$ position with respect to the lower right corner of the cardboard support was recorded for each measurement.

Fig. 2

Spatial registration of Raman spectroscopy measurements with histopathology: (a) Representation of the technique where a fresh specimen was initially pinned to a cardboard support placed on a guiding template in the form of a grid. The circle superposed with the specimen represents the position of the probe above the grid. (b) End result of the methodology leading to the spatial registration of all spectroscopic measurements with an H&E-stained section. Different colors were used for some of the circles to improve visualization in cases when they overlapped.

The laser power of the system was set at a value that resulted in 100 mW being delivered to the tissue surface. The number of repeat measurements per point (accumulations) was fixed at $N=10$ and the exposure time per spectrum ranged from 0.1 to 4 s, including dark noise (background) measurements made with the laser off. The exposure time was adjusted to ensure a raw signal intensity above 60% of the sensor dynamical range while avoiding saturation. All measurements were acquired within a one-hour timeframe after surgical excision and the size of the samples was limited to $\sim 2{\mathrm{cm}}^2$ to preserve the integrity of the specimen and the tumor for margin evaluation, as part of the regular post-surgery clinical diagnostic workflow. The size limitation and the time constrain limited the number of independent Raman measurements to—on average—10 per sample (i.e., 20 per patient). After spectroscopic measurements, the samples were fixed in formalin, embedded in paraffin, and sectioned in 3 to $10\mu \mathrm{m}$ histological sections. They were then stained with H&E as per institutional standards, resulting in one stained image for each sample.

2.4.

Registration with Optical Measurements and Histology LABELS

A four-step methodology was developed leading to spatial registration of all Raman spectroscopy measurements with corresponding locations on the H&E images using the Inkscape software (Inkscape’s Contributors). First, macroscopic photographs of the breast sample on the cardboard support [Fig. 2(a)] were taken. Second, the H&E images were superposed to the photograph and rotated to ensure the contours of the stained tissue sections matched as much as possible the contours seen on the photograph. Irregularities in the specimen shape were used as geometrical landmarks guiding this spatial registration process. Visible structures, including blood vessels and regions with high adipose content, were also used as tissue-based fiducial markers to ensure an accurate superposition of the photographs with the H&E images. Third, a numerical version of the registration grid was superposed to the H&E image, and it was used to pinpoint the location of every spectroscopic measurement. Fourth, each measurement location was annotated with a 3-mm diameter circle. The size of that circle was set to be two times larger when compared to the tip of the probe, accounting for potential spatial registration inaccuracies associated with sample deformations. As exemplified in Fig. 2(b), this led to images where colored circles indicated measurement locations superposed with the H&E images. Different colors were used for some of the circles to improve visualization in cases when they overlapped.

The resulting images were analyzed by pathologists (D.T., F.A.) to determine which types of cells were contained within each circled region. Cells were classified into three categories: (1) cancer (tumor cells, tumor stroma, or necrosis); (2) normal: either normal breast (connective tissue, stroma, fibroblast, and collagen) or breast parenchyma (ducts and lobules); (3) fat (adipose cells). The number of cells from each category was counted and the percentage of cancer, normal, or fat cells within each circle were computed. A measurement was labeled as cancer if the percentage of cancer cells was equal or superior to 80%, normal if the percentage of normal cells was equal or superior to 80%, and fat if the percentage of adipose cells was equal or superior to 70%. All measurements that did not fit into one of the three classes were excluded.

2.5.

Data Processing and Labeling

The data processing steps preceding the production of machine learning models led to the extraction of the inelastic scattering signature, for each measurement, from background contributions, including intrinsic tissue fluorescence.³⁷ Spectral pre-processing included the following steps: (1) averaging of the $N=10$ repeat spectra to increase overall signal-to-noise ratio, (2) averaging and subtraction of background spectra acquired with the laser turned off between measurements, (3) normalization with a NIST Raman standard (SRM 2214) to correct for the instrument response, (4) $x$-axis (wavenumber shift) calibration based on a Raman spectrum acquired on polycarbonate, (5) removal of background signals using the custom background removal algorithm BubbleFill,³⁸ and (6) standard normal variate (SNV) normalization.

A quantitative quality factor metric (between 0 and 1) was then computed from each resulting spectrum. It provided a statistical assessment of the likelihood the SNV-normalized signal is associated with tissue Raman peaks or stochastic noise.³⁹ Spectra with a quality factor metric inferior to 0.6 were excluded as they did not always clearly show all the ubiquitous Raman peaks encountered in biological tissue (e.g., amide bands for protein-rich tissue). All spectra were then visually inspected individually to qualitatively identify any glaring labeling errors here, mostly based on the fact that adipose tissue Raman spectra have spectral features that are dramatically different when compared to protein-rich tissue. For example, pure fat does not contain phenylalanine and as a result does not present a visually detectable peak at $1004{\mathrm{cm}}^{-1}$. Spectra labeled as cancer or normal that were associated with clearly identifiable fat spectra were removed from the dataset, as well as spectra identified as fat that did not show the expected biomolecular features associated with proteins. All remaining spectra were then plotted along with the standard deviation for each of the following groups: normal, cancer, and fat. The average spectrum associated with the normal breast measurements of the breast reduction surgery patient was also plotted against all other normal spectra and the cancer spectra (Fig. 7).

2.7.

Machine Learning Workflow

A machine learning workflow was applied to develop three different models: cancer versus normal + fat (model A), cancer versus normal (model B), and cancer versus fat (model C). Due to the small size of the dataset and to minimize the risk of over-fitting the data, the first step of the machine learning workflow consisted of dimensionally reducing the number of features per spectra to $<10$ for each model. This dimensional reduction procedure was done using a linear support vector machine (SVM) with a regularization by Lasso regression, in which individual statistical weights were assigned to each feature. This statistical weight acted as a surrogate for feature relevance, in effect quantifying the ability of each feature to capture inter-class variations and allowed less important bands to be discarded.^40,41 Features were then ranked according to their weight and boxplots were produced to show the statistical distribution of spectral intensity for all features that were retained to produce the machine learning models (Table 2 and Fig. 5).

Table 2

Raman peaks and main corresponding vibrational bonds associated with the Raman-predicted molecular tissue content associated with the spectra in Fig. 3. The bands associated with spectral features used as inputs to the classification models (Model A, Model B, and Model C) are identified in the last column. A tentative molecular assignment (specific molecules and families of biomolecules) based on literature findings is shown.42

The classification models training process involved SVM with a linear kernel and a regularization parameter $C$. This hyperparameter, which varied between ${10}^{-3}$ and 1, controlled the penalty for errors in the training process and reduced the risk of overfitting the data. The SVM algorithm also considered the class imbalance for all datasets by setting a higher cost $\gamma$ to misclassified spectra from the least populated class.⁴¹ A grid search method was used to optimize both hyperparameters $C$ and $\gamma$, and a five-fold cross-validation procedure was performed to assess classification performances. The five-fold cross validation took the training dataset and randomly split it in five equal-size groups. Then, five models were created with the same hyperparameters. Each model was trained on four groups and tested on the remaining one. The average of the performance of the five classifiers performance was used to assess predictive performance.⁴³ During the cross-validation phase, data from the same patient were never used in both the training and testing sets.

The technique returned, for each measurement, a posterior probability $0\le p\le 1$ that a measurement was classified within one of two classes. A receiver-operating-characteristic (ROC) curve was computed by comparing this posterior probability $p$ with a parameter $\lambda (0\le \lambda \le 1)$. All observations associated by the classifier were assigned the label 1 when $p\ge \lambda$ or assigned the label 0 if otherwise. Different values of the parameter $\lambda$ corresponded to different points of the ROC curve. Sensitivity and specificity were computed by comparing the label assigned by the model to the label given by the histopathology analysis for every measurement. Results were reported for the optimized hyperparameters only, yielding the highest area-under-curve (AUC) value. The final optimized model, with reported accuracy, sensitivity, and specificity, corresponded to the ROC curve point with the smallest distance to the upper left corner of the curve ($x$-axis: 1 − specificity, $y$-axis: sensitivity).

3.1.

Spectroscopic Measurements

Application of the Raman spectroscopy measurement protocol resulted in 388 spectroscopic measurements with co-located histopathology analyses. Of those, 58 spectra were excluded because of visibly low spectral quality, non-tissue-related artifacts (e.g., from residual ambient light, Raman peaks from the cardboard substrate) or a quality factor lower than 0.6. Moreover, 58 measurements were excluded because they did not qualify for association with either the normal, cancer, or fat categories. Further, 34 spectra were excluded due to suspected labeling errors based on visual inspection comparing measured spectra with the known Raman signature of adipose tissue. Of those, 28 that were labeled as normal or cancer were excluded because the spectra were identical to fat spectra. The remaining six spectra, which were labeled as fat, were excluded because they showed clearly visible spectral features associated with protein-rich tissue. The final dataset was composed of 238 measurements: 58 were labeled as normal, 87 as cancer, and 93 as fat. The raw and processed average measurements for each tissue category, along with the signal variance associated with each spectral bin, are shown in Figs. 3(a) and 3(b), respectively. Selected bands of interest are identified in Fig. 3(b) and their biochemical interpretation is provided in Table 2.

Fig. 3

Spectra associated with normal, cancer, and fat measurements in breast tissue: (a) average spectroscopic measurements (raw data) for each tissue category along with the standard deviation for each spectral bin; (b) average Raman spectra and standard deviation with selected bands identified with their wavenumber value. All identified bands were used for the biochemical interpretation of the data (Table 2), and a subset of those bands was used for machine learning model development.

3.2.

Machine Learning Models and Biomolecular Predictions

A classification model (Model A) was trained to distinguish between the cancer and non-cancer (fat + normal) spectra, from 238 measurements (87 cancer spectra, 58 normal spectra, and 93 fat spectra). This yielded an accuracy of 94%, a sensitivity of 93% and a specificity of 95%. The ROC curve for this model is shown in Fig. 4(a) and it is associated with an AUC of 0.98. These results were achieved using a model trained using only two spectral features, namely the peaks mostly associated with C–C stretching of proteins around $940{\mathrm{cm}}^{-1}$ and the symmetric aromatic ring breathing at $1004{\mathrm{cm}}^{-1}$ that is associated with the presence of phenylalanine. A classification model (Model B) was also trained to distinguish cancer from normal spectra from 145 measurements (87 cancer spectra and 58 normal spectra). This model had an accuracy of 91%, a sensitivity of 92%, a specificity of 90%, and an AUC of 0.96 based on four spectral features [Fig. 4(b)]. Those features corresponded to the same features as Model A, plus a peak associated with C–C skeletal stretching from proteins, lipids, and carbohydrates ($1129{\mathrm{cm}}^{-1}$) and a band around $1157{\mathrm{cm}}^{-1}$ associated with C-C/C-N skeletal stretching from proteins. Model C was developed to detect cancer from fat based on 180 Raman measurements (87 cancer spectra and 93 fat spectra). It performed with an accuracy of 98%, a sensitivity of 97%, and a specificity of 98%, with an AUC of almost 1 [Fig. 4(c)]. This model required two spectral features, namely, the $1004{\mathrm{cm}}^{-1}$ peak and the ${\mathrm{CH}}_2/{\mathrm{CH}}_3$ twisting deformation band at 1301 to $1304{\mathrm{cm}}^{-1}$.

Fig. 4

Average Raman spectra and spectral features (dotted lines) used by the classification models. Classification results are expressed in term of accuracy, sensitivity, and specificity. Confusion matrices and ROC curves are also shown for the three models: (a) Model A: cancer versus non-cancer (normal + fat); (b) Model B: cancer versus normal, and (c) Model C: cancer versus fat.

The spectral features that were used by the machine learning models are given in Table 2 along with their associated vibrational bonds and a tentative assignment with biomolecules based on the literature. Figure 5 shows the statistical distribution of these bands in the form of boxplots. The median (second quartile, Q2), first quartile (Q1), and third quartile (Q3) were clearly distinct for the two features used in model A (940 and $1004{\mathrm{cm}}^{-1}$). Q1, Q2, and Q3 were also well separated for cancer and normal tissue features associated with model B (940, 1004, 1129, and $1155{\mathrm{cm}}^{-1}$) as well as for cancer and fat for the features used in model C (1004 and $1304{\mathrm{cm}}^{-1}$).

Fig. 5

Statistical distribution of the intensity of the five features used by classification models for fat, normal and cancer spectra. Model A (cancer versus normal + fat) used features associated with bands at 940 and $1004{\mathrm{cm}}^{-1}$; Model B (cancer versus normal) used bands at 940, 1004, 1129, and $1155{\mathrm{cm}}^{-1}$; and Model C (cancer versus fat) used bands at 1004 and $1304{\mathrm{cm}}^{-1}$. Outliers are represented by black diamonds.

A fourth model was trained to distinguish fat from all other tissue types (normal + cancer) from 238 spectra (Fig. 6), resulting in an almost perfect AUC of 1 with specificity and sensitivity $>98\%$. This model was used to test whether the exclusion criteria leading to the rejection of suspected mis-labeled samples (28 labeled as normal or cancer, 6 spectra labeled as fat) was justified. The model was directly applied to those 34 spectra and predicted, in all cases, that the histology labels were inconsistent with the measured Raman signature. Specifically, the model predicted that all 28 normal or cancer spectra were actually fat, and that all 6 fat spectra were actually normal or cancer.

Fig. 6

A Fat detection classifier was trained from 238 spectra. It was applied on the 34 spectra that were excluded from the study because of suspected mis-labeling likely associated with spatial registration errors of Raman spectroscopy measurements with the H&E images: the 28 spectra that were labeled as normal or cancer were all classified as fat (false negatives) and the six spectra that were labeled as fat were classified as normal or cancer (false positives).

In this study, more than 30% of all normal measurements were associated with the breast conserving surgery patient, which could have led to classification biases. However, Fig. 7 qualitatively shows that the average normal spectrum from the breast conserving surgery patient ($n=18$ spectra) was similar to the average spectrum associated to all other normal spectra in the dataset ($n=40$). More importantly, the average spectra associated with the normal category matched perfectly for the Raman bands used to train Models A, B, and C (Table 2), providing strong evidence no bias resulted from including a larger proportion of normal measurements from the same patient.

Fig. 7

Raman spectra associated with the normal category from the breast reduction surgery patient, compared to all other normal measurements. Cancer spectra are also shown for visual comparison. The four features used in the classification models are identified with dotted lines.