Multi-level multi-modality (PET and CT) fusion radiomics: prognostic modeling for non-small cell lung carcinoma

Our study was conducted on NSCLC datasets obtained from the Cancer Imaging Archive (TCIA) (Prior et al 2013, Bakr et al 2018, Shiri et al 2020b), including 211 histologically proven NSCLC patients from two independent institutions namely Palo Alto Veterans Affairs Healthcare System, and Stanford University School of Medicine (referred to as VA and Stanford in the rest of this work). Two hundred and one patients underwent ¹⁸F-FDG PET/CT scans (10 patients were excluded owing to unavailability of PET images), prior to surgical treatment between April 2008 and September 2012. Our exclusion criteria were: visual observation of images with high levels of noise and/or artifacts, miss-segmentation, and miss-registration errors. Overall, 182 patients (87 from VA and 95 from Stanford) were included. Imaging data extracted from the DICOM headers were as follow. A GE Discovery PET/CT scanner was used to acquire images in VA, while Stanford used GE Discovery D690, Philips Healthcare and Siemens Healthcare PET/CT scanners. Key acquisition parameters are reported separately for VA and Stanford datasets in supplementary table S1 (available online at stacks.iop.org/PMB/66/205017/mmedia). Imaging data were also paired with clinical information such as histopathological grade, TNM stage, age, sex, smoking status and survival time outcomes. However, TNM stage and histopathological grade were not available for 28 subjects in the Stanford dataset. Characteristics of patients categorized in VA and Stanford datasets are described in table 1.

For radiomics analysis on multi-modal PET/CT imaging data, PET and CT data underwent fusion at feature- and image- levels (figure 1). For image-level fusion, features were extracted from a single image, obtained from the merging of PET and CT images via wavelet fusion strategy. For feature-level fusion, features were extracted from separate scans and then two different approaches of feature concatenation and feature averaging were pursued. Individual PET and CT models were also constructed for comparison purposes. Moreover, all models were considered with and without clinical features in order to assess the additive prognostic value of clinical parameters to each model.

Zoom In Zoom Out Reset image size

All malignant lesions in PET images were delineated and segmented manually using OSIRIX^®. 3D-slicer was utilized to segment tumor volumes from CT images using automatic region growing. Both PET and CT segmentations were edited and verified by an experienced radiologist. Subsequently, CT and PET segmentations were merged into a single mask to reduce segmentation errors and ensure that the features are extracted within the same volume of interest (VOI) over all single-modality and fusion models. In the merged segmentation, we generated a binary mask, where each voxel was assigned a label of 1, if the corresponding voxels in PET or CT masks were equal to 1 (supplementary figure S1). In addition, ∼10% of the patients had more than one lesion. Supplementary figure S2 illustrates the number of lesions for each patient and the procedure of analyzing multiple lesions together.

For image-level fusion, a wavelet-based technique was adopted. In order to integrate the spatial and frequency characteristics of CT and PET imaging modalities, the 3D discrete wavelet transform (DWT) was utilized. However, prior to image fusion, the volumes need to be registered to ensure that all voxels from one volume are mapped to their corresponding voxels in the other volume. Since PET and CT scans were captured simultaneously on a hybrid PET/CT scanner with no major motion artefacts, they were reasonably well aligned. However, in this dataset, not only PET and CT images had different in plain resolutions, but also, in many patients the size, and the position of in plain field-of-view was different (images' origins had to be matched). To address this issue, 'Image Position' in DICOM headers was utilized to crop and/or pad scans to the same size to get a registered field-of-view. Subsequently, PET and CT images were resampled to an equal resolution with isotropic voxel spacing of 2 × 2 × 2 mm³ using cubic interpolation and an antialiasing kernel. Pixel spacings of the original CT and PET images were within ranges of [0.88–1.37] mm and [3.65–5.47] mm, respectively. Hence, downsapmling CT to PET image resolution, or upsampling PET to CT image resolution would significantly affect the intensity values of PET or CT images. As such, we used 2 × 2 × 2 mm³ isotropic voxel size to keep both images as intact as possible.

Subsequently, after preparing images for fusion, PET and CT images were decomposed up to one level by applying 3D DWT using wavelet basis function symlet8. Subsequently, the spatially corresponding coefficients of the wavelet domains of images were weight averaged to obtain a single set of fused wavelet coefficients. Finally, the set of fused wavelet coefficients underwent 3D inverse DWT to obtain a fused image.

We investigated the influence of two wavelet fusion parameters on the performance of the models constructed from the consequent fused image. One adjustment modified the weights of CT and PET corresponding wavelet coefficients in the averaging procedure. CT weights (W) were set to 1/4, 1/2, and 3/4; and the corresponding PET weights were (1—CT weights). The other adjustment was applied through wavelet band-pass filtering (WBPF). This was accomplished by allocating different weights to band-pass sub-bands (HHL, HLH, LHH, HLL, LHL, LLH) of the wavelet domain, in comparison to the low- and high-frequency sub-bands (HHH, LLL). The ratio (R) of the weight assigned to band-pass sub-bands, over the weight assigned to the other sub-bands, was set to 2/3, 1, and 3/2. This band-pass filtering procedure was proposed by Carrier-Vallières (Carrier-Vallières 2013), suggesting that the spatial definition of textures in the fused images can be enhanced by modification of corresponding frequency properties of the images in the wavelet domain. Nine different fused images were constructed from all possible combinations of three CT weights and three wavelet band-pass ratios (a specific fused image with CT weight of 1/2 and band-pass ratio of 3/2 is denoted as WF_R3/2_W1/2) (figure 2). In the case of feature-level fusion, two different strategies were considered. The first method involved concatenation of features from different modalities, which is a common feature-level fusion strategy (ConFea). For the second strategy, features from the two modalities were averaged, and then the model was constructed from the above-mentioned averaged feature set (AvgFea). Thirteen different radiomic models from three different strategies, namely single-modality strategy (individual CT and PET), feature-level fusion strategy (AvgFea and ConFea), and image-level fusion strategy (nine wavelet fusion models) were analyzed in this study.

Zoom In Zoom Out Reset image size

In addition to radiomics models, a clinical model was constructed from two clinical parameters, namely TNM staging and histopathologic grade. Other clinical parameters such as sex, age, treatment regimen and smoking status had very limited univariate performance (C-indices ∼0.5). By comparison, TNM staging and histopathological grade had higher performance in our univariate model (C-indices of 0.63 and 0.59, respectively) and were selected for clinical modeling. Moreover, all 13 radiomic models were also extended to contain both radiomics and clinical features. To ensure the presence of clinical features in all radiomics + clinical models for thorough comparisons, radiomics features first went through the feature selection procedure and then the model was constructed from concatenation of the selected radiomics features and the two clinical parameters. Altogether, we investigated 26 + 1 different models including 13 radiomics, 13 radiomics + clinical strategies, plus one clinical model based on clinical-only (TNM staging and grade) parameters.

To standardize images with different voxel sizes and obtain rotationally invariant texture features, images were interpolated to isotropic voxel spacing of 2 × 2 × 2 mm³. Then, prior to computation of texture features, intensity levels inside the tumor volume were discretized into 16, 32, and 64 bins, so each radiomics model was extended to three models resulting in a total of 79 different models (3 × 13 radiomic models, and 3 × 13 radiomics + clinical models as mentioned earlier, plus one clinical-only model) to be analyzed. From each VOI, 225 radiomic features were extracted. The feature-set included 79 first-order features (statistical, histogram, intensity-histogram and morphological features), 136 three-dimensional textural features extracted from gray-level co-occurrence matrix (GLCM), gray-level run length matrix (GLRLM), gray-level size zone matrix (GLSZM), gray-level distance zone matrix (GLDZM), neighborhood gray tone difference matrix (NGTDM), and neighborhood gray-level dependence matrix (NGLDM), and finally 10 moment invariant features. All features were extracted using The Standardized Environment for Radiomics Analysis (SERA) Package ¹⁰ (Ashrafinia 2019), a Matlab^®-based framework, in agreement with guidelines from the Image Biomarker Standardization Initiative (IBSI) (Zwanenburg et al 2016) and evaluated in multi-center standardization studies (McNitt-Gray et al 2020, Zwanenburg et al 2020), for improved features reproducibility.

Since two independent datasets from two different institutions were involved in this study, in order to validate our results more comprehensively, we investigated two train/test arrangements (or partitions); namely, using VA for training and Stanford for testing cohorts (referred as VAS), and alternatively switching datasets (using Stanford for training and VA as test cohort, referred to as SVA). In order to minimize redundancy within feature subsets, one of each features pair with Spearman's rank correlation coefficient higher than 0.9 were omitted. For each feature in each model, univariate Cox proportional hazard (Cox PH) regression was conducted by performing 100 repetitions, using bootstrap resampling. Then the prognostic performance of each feature was measured using Harrell's concordance index (C-index), and top 10 most relevant features with the highest C-indices were selected for each model to use in further multivariate analysis. For each strategy, we tested the performance of multivariate Cox models constructed from all possible combinations of candidate features (3 up to 10 features) in the training cohorts, and the combination with the highest C-index was selected as the optimum. The number of all possible combinations of features in this study containing 10 candidates was $\left(\begin{array}{c}10\\ 3\end{array}\right)+\left(\begin{array}{c}10\\ 4\end{array}\right)+\left(\begin{array}{c}10\\ 5\end{array}\right)+\left(\begin{array}{c}10\\ 6\end{array}\right)+\left(\begin{array}{c}10\\ 7\end{array}\right)+\left(\begin{array}{c}10\\ 8\end{array}\right)+\left(\begin{array}{c}10\\ 9\end{array}\right)+\left(\begin{array}{c}10\\ 10\end{array}\right)=967.$ Testing all possible combinations of candidate features is time-consuming, but aids in finding the optimum feature signature. For radiomics + clinical models, the two clinical parameters were added to the optimum feature set, and then fed to the multivariate Cox model. The optimum model of each strategy was applied to training and testing cohorts by 1000 repetitions using bootstrap resampling. After model was trained, in each boot, model was applied to 100 subjects selected (by replacement) from training cohort and then tested on 100 subjects selected (by replacement) from the testing dataset. The final C-index of each model was obtained by averaging 1000 C-indices captured by bootstraps. Significant difference between each pair of models was obtained by applying Student's t-test, to compare the mean values of C-indices of each model captured from bootstrapping on the testing cohort. In order to set strict acceptance, only t-test p-values lower than 0.001 were considered significant. We used the median prognostic score of each strategy achieved from the training cohorts, untouched as a threshold to classify patients into low- and high-risk groups. Then log-rank test was applied to find the difference between Kaplan–Meier curves of the subsets. In order to average and combine results from both VAS and SVA groups, C-index values were averaged, and log-rank p-values were combined using Fisher's method (Fisher 1992). P-values lower than 0.05 were considered statistically significant. All statistical analyses were performed on R (R-3.5.3).

The top 10 features with the highest performance (ranked by univariate Cox regression) were collected in each radiomics model. Supplementary figure S3 demonstrates the feature popularity. Large distance high gray level emphasis from GLDZM, busyness from NGTDM, dependence count variance from NGLDM, and global intensity peak were among the top ranked features of more than 20 out of a total of 39 radiomic models (including single- and multi-modality models). Moreover, the final selected features for each model in VAS partition are listed in supplementary table S2. The number of selected features ranges from 4 to 10 (mostly within the range [6–8], where only one model was trained with 10 features).

Figure 3 depicts a visual example of textures displayed on PET, CT and fused images, as well as feature maps of one of the most repeated features, NGTDM-busyness, for both a low- and a high-risk patient.

Zoom In Zoom Out Reset image size

The comparison of features between different image-level fusion models and PET/CT single-modalities was performed in the VAS partition for varying number of gray-level intensity discretization (16, 32, 64). However, since morphological features were constant through all models, the comparison was applied to 182 intensity and texture features. Supplementary table S3 shows the correlation between corresponding features in PET/CT single-modality models and different image-level fusion models for 16, 32, and 64 gray-level discretization, respectively.

Here we report the performance of each model for OS prediction with C-index (regression task) and log-rank p-value (classification task) in the test cohort (independent and unseen dataset for external validation). Table 2 summarizes the C-index, its standard deviation, and p-value for each radiomics, radiomics + clinical, and clinical-only model, by averaging the VAS and SVA results. Moreover, the results of all models can be found in a more detailed way, separately for VAS and SVA arrangements, in supplementary table S4. Supplementary figure S4 illustrates the box-plots of C-indices for all radiomics (supplementary figure S4(a)) and radiomics + clinical models (supplementary figure S4(b)). Figure 4 depicts Kaplan–Meier curves (for VAS training/testing arrangement) of models with the best performance in classifying patients into low- versus high-risk groups, among single-modalities, feature-level and image-level fusion strategies, as well as the clinical-only model. As shown in table 2, considering the regression task, the best performance was achieved by WF_R3/2_W1/2_16 (C-index = 0.689) in radiomics group, and WF_R2/3_W1/4_64 (C-index = 0.683) in radiomics + clinical group. The clinical model also reached a C-index of 0.656. Note that in specific partition VAS (supplementary Table S4), several image-level fusion models showed significantly higher performance. WF_R3/2_W1/2_16 (C-index = 0.735) from radiomics and WF_R1_W1/2_16 (C-index = 0.762) from radiomics + clinical groups had the best results.

Table 2. The C-index, its standard deviation, and p-value of different models in testing cohort averaged for both VAS and SVA partitions. C-indices higher than 0.65 are marked as red, and p-values lower than 0.01 are color marked as blue.

Zoom In Zoom Out Reset image size

In the case of classification competency, in radiomics models, CT-16 from individual strategy, AvgFea-16 from feature-level fusion strategy and 10 models from image-level fusion strategy such as WF_R2/3_W3/4_16, WF_R1_W1/4_32 and WF_R3/2_W1/2_64 reached p-value <0.001. In radiomics + clinical group, WF_R1_W3/4_16 (p-value: 0.003) had the best performance. Moreover, clinical model was able to classify patients in low- versus high-risk groups with an acceptable p-value of 0.001. Further results are aimed to provide a deeper insight into our findings by comparing different models from different strategies (clinical, single-modality, feature- and image-level fusion) in radiomics and radiomics + clinical groups, separately. The corresponding p-values of the t-tests applied to all possible pairs of models can be found in supplementary figure S5 for radiomics versus radiomics supplementary (figure S5(a)), radiomics versus radiomics + clinical (figure S5(b)), and radiomics + clinical versus radiomics + clinical (figure S5(c)). Figure 5 compares the best models of each strategy for radiomics (figure 5(a)) and radiomics + clinical (figure 5(b)) groups. Moreover, a comparison of strategies (comparison of average c-index of models from single-modality strategy, feature-level fusion strategy, image-level fusion strategy, and clinical model) is presented in supplementary figure S6 for radiomics and radiomics + clinical, respectively.

Zoom In Zoom Out Reset image size

In radiomics group (figure 5(a)), the best models in single-modality strategy were CT-64 and PET-16 (C-indices = 0.631 and 0.635, respectively). Regarding feature-level fusion strategy, the best model was AvgFea-64 (C-index = 0.629); although the difference with CT-64 was not significant (p-value = 0.308), it was predominated by PET-16 (p-value < 0.001). Regarding image-level fusion strategy, WF_R3/2_W1/2_16 (C-index = 0.689) significantly outperformed CT-64, PET-16, and AvgFea-64 (p-values <0.001).

Amongst radiomics + clinical models (figure 5(b)), the best results of single-modality strategy were again achieved by CT-64 and PET-16 (C-indices = 0.617 and 0.605, respectively). Regarding feature-level fusion strategy, the best model was AvgFea-64 (C-index: 0.612) which was outperformed by CT-64, but outperformed PET-16 (both p-values <0.001). Considering image-level fusion strategy, best model was WF_R2/3_1/4_64 (C-index = 0.683), which significantly outperformed CT-64, PET-16, and AvgFea-64 (p-values <0.001).

The clinical model (C-index = 0.656) outperformed all single-modality and feature-level fusion models but was outperformed by some image-level fusion models, such as WF_R3/2_W1/2_16 and WF_R2/3_1/4_64 from radiomics and radiomics + clinical groups, respectively.

Strategies were also evaluated and compared by applying Student's t-test on the average C-index of models in different strategies (supplementary figure S6). While the results of feature-level fusion strategy and single-modality strategy were not significantly different (p-values = 0.474, and 0.690), image-level fusion strategy significantly outperformed single-modality (p-values <0.01) and feature-level fusion strategies (p-value <0.001 and p-value <0.005) when comparing within radiomics and radiomics + clinical groups, respectively. Moreover, the result of the clinical model was significantly higher than single-modality, and feature-level fusion strategies, but was comparable with image-level fusion strategy (p-values = 0 .749 and 0.367 for radiomics and radiomics + clinical group, respectively).

Figure 6 illustrates the scatter plots of results from radiomics (figure 6(a)) and radiomics + clinical (figure 6(b)) groups, based on their strategies (single-modality, feature- and image-level fusion), with Spearman's rank correlation coefficients and the corresponding p-values indicating statistically significant differences. Spearman's rank correlation test was performed between ranks of different strategies (single-modality, feature- and image-level fusion ranked 1, 2, and 3, respectively) versus corresponding C-indices. The correlation coefficients between single-modality and feature-level fusion ranks versus their C-indices were 0.19 and zero for radiomics and radiomics + clinical, respectively (p-values for significant correlation were 0.55 and 1.0), the correlation coefficients between individual and image-level fusion ranks versus their C-indices were 0.47 and 0.45 for radiomics and radiomics + clinical groups (p-values < 0.01 for both), and finally the correlation coefficients between feature- and image-level fusion ranks versus their C-indices were 0.44 and 0.42 for radiomics and radiomics + clinical, respectively (p-values = 0.01 for both).

Zoom In Zoom Out Reset image size

Statistically significant differences between the performance of each radiomics model and its peer within radiomics + clinical models were determined in the VAS partition, and the corresponding p-values are presented in supplementary table S4. Out of a total of 39 radiomics models, in 14 models, the performance was significantly improved; in 13 other models, the performance was reduced, whereas the performance of the remaining 12 models did not differ significantly.

In order to assess the effect of wavelet-fusion parameters (WBPF ratio, and CT-weight) and number of discretized gray-levels on the performance of the image-level fusion models, we performed Spearman's correlation between values of a specific parameter (for example 2/3, 1, 3/2 for WBPF ratios) and C-indices of the resulting models in order to see if there is any increasing or decreasing trend. The scatter plot of the results obtained from image-level fusion models along with trend-lines fitted on them is displayed as a function of different WBPF ratios (supplementary figure S7(a)), different CT-weights (supplementary figure S7(b)), and different number of discretized gray-levels (supplementary figure S7(c)). The correlation coefficients for C-indices versus WBPF ratios R, CT-weights W, and number of discretized gray-levels were −0.256, −0.122, and −0.099, respectively (p-values of significant difference with zero were equal to 0.197, 0.543, and 0.623, respectively).