2.1.

Patient Selection and Sample Collection

This study curated 149 patients diagnosed with PM from the University of Chicago Medicine (UCM) under a Health Insurance Portability and Accountability Act–compliant, Institutional Review Board–approved protocol from April 2016 to June 2022. Informed consent was obtained from all participants. The protocol allowed for the collection and biobanking of peripheral blood, saliva, and tumor samples. Tumor deoxyribonucleic acid was extracted from fresh frozen, paraffin-embedded tumor tissue blocks. Somatic mutations were identified using the UCM OncoPlus next-generation sequencing panel.²¹ Patients with confirmed somatic BAP1 mutations only (BAP1+, $N=68$) were included in the study. The remaining 81 patients presented with the wild-type allele (BAP1−). Immunohistochemical analysis of the BAP1 protein was conducted in a Clinical Laboratory Improvement Amendments–certified laboratory at UCM using the Santa Cruz C4 monoclonal antibody. Table 1 includes further details about the patients of this study.

Table 1

Patient demographics categorizing patient sex and age characteristics.

2.2.

Image Data Curation and Segmentation

Axial images from unenhanced chest CT scans of the patients were retrospectively collected (Table 2).²² Scans were acquired with the assistance of the University of Chicago’s Human Imaging Research Office,^23,24 which provided de-identified, compliant images for evaluation. For each patient, the CT section displaying the largest area of tumor was selected by a radiologist. This section and the immediate superior and inferior sections were used to create a three-dimensional (3D) volume for analysis. A Visual Geometry Group 16 (VGG16)/U-Net deep convolutional neural network (CNN) architecture was utilized to segment the tumor within this volume.²⁵ The two-dimensional (2D) architecture employed downsampling and upsampling paths. The downsampling path utilized a VGG16 model pre-trained on ImageNet with scale-jittering, applying $2\times 2$ max pooling with stride 2. Dropout layers of 0.5 probability were used to prevent model overfitting. The upsampling path employed a 2D operation with nearest-neighbor interpolation on the feature maps. The network generated $512\times 512\text{-}\mathrm{pixel}$ probability maps, which matched the input image size. Rectified linear unit and sigmoid activation functions were applied after the convolutional layers and the final layer, respectively. Lastly, the model was trained with a binary cross-entropy loss function using the Adam optimizer with a learning rate of ${10}^{-5}$. More details regarding the architecture of the model and its training scheme can be found in Gudmundsson et al.²⁵ For the present study, tumor contours were automatically generated and evaluated with no additional training or validation of the model.

Table 2

Image acquisition characteristics for the patient cohort analyzed in this study.

The resulting probability maps output by the network were thresholded at a value of 0.2; this threshold was determined to have maximal overlap with human contours using the Dice similarity coefficient (DSC) from prior work.^26,27 The radiologist adjusted the resulting segmentations to ensure the segmentations were highly specific to tumor pixels. The finalized contours were defined as the region of interest (ROI) and used for feature extraction.

2.3.

Image Resampling and Gray-Level Discretization

To mitigate the impact of different image acquisition parameters, all images were resampled to the mean resolution of all scans, with a pixel spacing of $0.75\times 0.75\mathrm{mm}$ and a slice thickness of 3 mm (see Table 2 and Fig. 2 for more details). Prior to texture feature extraction, gray-level discretization was applied using a fixed bin number of 32 gray levels, as small or large gray-level quantizations have been shown to impact texture feature values due to the reduction of information that can be extracted from an image.^28,29

Fig. 2

Histogram of the (a) pixel spacing and (b) slice thickness of CT sections of the original 149 scans. The red vertical line depicts the mean value in each of the distributions to which resampling was performed.

2.4.

Feature Extraction

Eighteen intensity-based and 123 texture features (111 second-order, 6 Laws’ texture energy, 2 Fourier, and 4 fractal dimensions) were extracted from the original ROIs. The 123 texture features were also extracted from the ROIs after applying seven different filtering operations on the images: two Laplacian of Gaussian (LoG) filters ($\sigma =0.75\mathrm{mm}$, 1.5 mm), four multi-channel wavelet decompositions [low low (LL), low high (LH), high low (HL), and high high (HH)], and a local binary pattern operator (radius = 0.75 mm). With 18 intensity-based features and 123 texture features extracted from the ROIs before the filtering operations and the 123 features extracted from the ROIs after the seven filtering operations, a total of 1002 features were computed from each ROI (the finalized tumor contours). Intensity-based features were obtained from the 3D volume.³⁰ All other features were computed by averaging the 2D feature values over the three CT sections. Features were calculated using the Python packages PyRadiomics,³¹ PyFeats, and Nyxus.

2.5.

Data Imbalance

Due to the imbalance of BAP1 mutation status among patients, a hybrid approach using the Synthetic Minority Over-sampling Technique (SMOTE) coupled with the removal of Tomek links was employed to over-sample the minority class and under-sample all classes,³² respectively, prior to the feature selection. The SMOTE algorithm generates artificial data in the feature space near the existing feature values of cases from the minority class. Tomek links are a pair of nearest neighbors of opposite classes with minimal distance between them compared with other neighboring data. Removal of Tomek links decreases noisy data or eliminates data near the decision boundary. Implementation of SMOTE–Tomek resulted in equal mutation prevalence, per fold, during training.

2.6.

Machine Learning Model and Feature Selection

The performance of 18 separate³³ calibrated machine learning models (Table 3) was evaluated using leave-one-out cross-validation (LOOCV), resulting in 149 iterations. Calibration was performed using the “sigmoid” method, which corresponds to fitting a logistic regression model to the scores of a classifier (Platt’s scaling). Although “isotonic” calibration, which fits a non-parametric isotonic regressor, could be performed, such calibration is recommended only for large datasets as overfitting could result in too few samples (i.e., fewer than 1000 cases).^34,35

Table 3

Types of models evaluated in the BAP1 classification task.

Feature selection was performed on the training set of each iteration of the LOOCV in the following order (with empirically determined parameters):³⁶ (1) features with variance less than 0.01 were discarded, (2) features were Z-score normalized, and (3) features with a Pearson correlation coefficient of 0.75 or greater with other features were removed (to assess linear independence among the features).³⁷ Lastly, the top four features were selected using the calculated F-statistic of the analysis of variance (ANOVA) test between the feature and the BAP1 mutation status.³⁸ These four features were then extracted from the left-out test case, per iteration, for the classification task.

Other training schemes were assessed. In particular, different-sized folds for repeated k-fold cross-validation were implemented as well as changing the number of top features selected.³⁸ Preliminary work was also performed to study the impact random state seeds had on the classification task.

2.7.

Evaluation Metric and Statistical Analysis

The receiver operating characteristic area under the curve (ROC AUC) was used as the figure of merit to assess the classification performances of the models to differentiate among BAP1+/− patients. The Wilcoxon rank-sum test was used to assess the differences in tumor volume, and age distributions between patients in the two classes and the DeLong³⁹ and Wilcoxon signed-rank tests were used to evaluate the differences in AUC values among the models. To assess the impact of human modifications on the segmentation of the PM tumor, DSC values were calculated between the CNN segmentations and radiologist-modified masks to determine the overlap between the two. Further, the classification task was performed employing the same models (Table 3) and using feature values extracted from the unmodified CNN probability maps thresholded at 0.2. Using the DeLong test, the AUC values computed from the unmodified segmentations were compared with the AUC values achieved from the modified segmentations. Due to the hypothesis-generating nature of this work, statistical significance was obtained at $p=0.05$.

3.1.

Tumor Volume

Figure 3 shows the distributions of the tumor volume contoured across the three sections selected per patient in the dataset; the median (range) volume of tumor contoured was $\mathrm{13,109}{\mathrm{mm}}^3$ (1630 to $\mathrm{108,331}{\mathrm{mm}}^3$) across all patients, $\mathrm{11,615}{\mathrm{mm}}^3$ (1630 to $\mathrm{108,331}{\mathrm{mm}}^3$) for BAP1+ patients, and $\mathrm{15,949}{\mathrm{mm}}^3$ (1688 to $\mathrm{92,352}{\mathrm{mm}}^3$) for BAP1− patients. The difference in median volume among the BAP1+/− patients failed to achieve statistical significance ($p=0.15$), which mitigated the impact of tumor size as a confounding factor for the classification task. BAP1+ patients had the larger range of tumor volumes, whereas BAP1− patients had the larger median. Differences in age among the patient cohorts failed to achieve statistical significance.

Fig. 3

Histogram of the tumor volume categorized by BAP1 mutation status. The difference in tumor volume between wild-type and mutated tumors failed to achieve statistical significance.

3.2.

Classification Performance

In the task of differentiating between BAP1+ and BAP1− patients, the top three models (sorted by AUC values) were decision tree, Gaussian process, and SVM with a radial basis function kernel (Table 3). Figure 4(a) shows the ROC curves obtained from the three models, along with their AUC values and the 95% confidence intervals (CIs). The AUC values and 95% CIs were constructed from 2000 bootstrapped samples of the prediction values during LOOCV. The decision tree classifier yielded an AUC value of 0.69 (95% CI: 0.60, 0.77). Figure 4(b) displays the distribution of scores obtained during the cross-validation for the top-performing model, the decision tree. No scores were less than 0.32 or greater than 0.70 for either class. The DeLong test failed to achieve a statistically significant difference in AUC values among the top three models as shown in Table 4.

Fig. 4

(a) ROC curves depicting the true-positive and false-positive fractions of the top three performing classifiers in the task of differentiating somatic BAP1 mutation status using feature values extracted from segmented regions. ROC curves were fitted using software created by Metz and Pan.⁴⁰ (b) Distributions of the decision tree classifier prediction scores across all cases. The histograms were normalized to have an equal area of 1.

Table 4

Comparisons of the three best-performing classification models: decision tree, Gaussian process, and support vector. The p-values comparing the differences in AUC values were calculated using the DeLong test, with their corresponding CIs. Significance levels (α) and widths of the CIs were adjusted for multiple comparisons. None of the comparisons achieved statistical significance after correcting for multiple comparisons using Bonferroni–Holm corrections.

The four features selected most frequently through the 149 iterations of the cross-validation are presented in Table 5. All the features were second-order [e.g., gray-level co-occurrence (GLCM) or gray-level size zone matrices (GLSZM)] and were extracted from LoG-filtered or wavelet-decomposed images.

Table 5

Four texture features most often selected during the 149 LOOCV iterations and the frequency each feature was chosen, i.e., the number of iterations in which a feature was selected.

3.3.

Change of k-Fold and Number of Features

Table 6 displays the AUC values achieved from the different number of folds used for the repeated k-fold cross-validation and the different number of features selected by the final ANOVA feature selection step: 200 repetitions were performed to ensure robust statistics for the calculation of the 95% CI. As reported in Sec. 3.2, the decision tree classifier resulted in the highest overall AUC value of 0.69 [0.60, 0.77]; however, this AUC value failed to achieve a significant difference ($p=0.1$) from the AUC value of the SGD classifier (0.63 [0.54, 0.72]) obtained when selecting the top 10 features, as determined using the DeLong test for correlated ROC comparison and setting the alternative hypothesis to “greater.”

Table 6

Model performance using various cross-validation approaches. ROC AUC values in the task of differentiating between BAP1+ and BAP1− patients and 95% CIs for the LOOCV were obtained using 2000 bootstrapped samples. For the 10-fold and 5-fold cross-validation, AUC values were acquired by averaging the AUC values per repeat of the cross-validation approach, and 95% CIs were obtained by calculating the 2.5% and 97.5% percentile of the distribution of AUC values.

A selection of four features yielded a different distribution of AUC values than the distribution of AUC values calculated with a selection of 10 features ($p<0.05$ as determined by the Wilcoxon signed-rank test). There was a significant difference between 10- and 5-fold cross-validation results when selecting four features ($p<0.05$); however, this trend did not occur for a selection of 10 features as there was a failure to achieve significance ($p=0.14$). Interestingly, the most-selected feature was the same across all cross-validation approaches: GLCM cluster prominence with an LoG filter applied of size $\sigma =1.5\hspace{0.17em}\mathrm{mm}$ (LoG_sigma = 2.0). The top-performing models encompassed different types, including ensemble, naive Bayes, discriminant analysis, neighbor, tree, and linear. Therefore, the classification schemes included all but the SVMs and Gaussian processes.

To assess the impact of the random state seed on the performance of a model, AUC values were recorded for 100 seeds of the decision tree classifier, resulting in a median AUC value of 0.66 [0.64, 0.68], with the 95% CI calculated using the percentiles for 2.5% and 97.5% of the distribution of the 100 AUC values calculated; the reported value of 0.69 obtained during LOOCV of the decision tree classifier was outside these boundaries of the CI constructed from the AUC values calculated for the 100 random seeds.

3.4.

DSC and Classification Performance of Unmodified Segmentations

When comparing the CNN segmentations to the radiologist-modified segmentations, an average DSC value of 0.79 with an interquartile range of 0.21 was achieved. The same feature extraction and selection were performed on the unmodified segmentations of tumor contours. The CNN failed to predict tumors for one case; therefore, that case was discarded from the analysis. Using LOOCV, the highest AUC value achieved across the 18 models was 0.61. The decision tree classifier, the highest-performing model as aforementioned, yielded an AUC value of 0.45 (0.36, 0.56), which was significantly different than 0.69 ($p<0.001$) as determined by the DeLong test.