Radiomics is the science of high-throughput mining of quantitative image information from medical images to reveal a spectrum of characteristics and information content in form of features that are not appreciable by the human eyes. ^{1 – 7} Many studies have shown that medical images contain abundant information about shape, intensity, and texture of the volume of interest that can reveal characterization of underlying pathology. ^{1 , 2 , 5 – 19} In modern medicine, such information can play an important role in detection of image-based biomarkers/signatures to enhance the predictive power of diagnostic and prognostic tools. Radiomics information constitutes a part of multi-omics, ²⁰ which typically consists of ample of features extracted from different information modalities. Radiomics features show promise toward serving as biomarkers of outcome of cancer and other human pathology. ^{1 – 20}

In the recent years, as many complex multi-omics datasets have been applied in different fields of healthcare, the importance of feature engineering process with capability of feature interpretation has increased. ²¹ Nowadays, feature engineering plays a key role in the field of multi-omics. Feature engineering enables mining of associations between input variables and relevant outcomes, including underlying pathology, to reveal feature- predictive biomarkers among different multi-omics data. ²¹

Feature engineering or multivariable feature analysis refers to the process of recruiting comprehensive statistical methods and domain knowledge to reveal the most relevant subset of features from the raw feature sets to be used in supervised or unsupervised predictive analytics processes. ^{22 – 28} Such analyses can be either performed in the spatial-domain (image feature space), or the frequency-domain to highlight specific information and reveal the most relevant and explanatory features towards the application of interest. The focus of the frequency-based feature engineering techniques is to decompose the raw information into different frequency bandwidths or information components in the frequency-domain which can be compared among different samples/cases and used for model development. ^{29 – 33} Multiresolution and filter-based radiomics analyses such as wavelet transform, ^{29 , 30 , 34 , 35} Fourier-based adaptive filter or transform, ^{36 – 39} Discrete Orthonormal Stock-well Transform (DOST), ^{40 , 41} Law’s textural information, ^{31 , 42} and Gabor filter banks ^{32 , 33} are examples of such frequency-based analysis approaches.

Wavelet decomposition (WD) techniques ^{43 , 44} form a more recent addition to the collection of multispectral signal analysis techniques. Unlike other frequency-based analysis methods such as Discrete Fourier Transform analysis, ^{38 , 39 , 45 – 48} Laplacian pyramids ⁴⁹ and etc., WD techniques provide a complete image representation by performing a full decomposition of the data according to scale and orientation, which are considered essential properties of feature description in the feature engineering process.

Our research group ^{50 – 53} among others ^{29 , 30 , 34 , 35 , 54 – 66} have shown that the frequency and wavelet analysis of multi-spectral features in the frequency domain can reveal more structured and relevant information details pertinent to the outcome of interest compared to the spatial domain. ^{29 , 30 , 36 , 48 , 67} We have recently demonstrated that radiomic features extracted from primary tumor gross tumor volumes (GTVs) delineated on contrast-enhanced (CE) CT images of oropharyngeal squamous cell carcinoma (OPSCC) patients can be used to construct a highly sensitive and specific classifier for characterization and prediction of human papilloma virus (HPV). ^{53 , 68 , 69} In a related study, we have also shown ^{36 , 50 – 52} that frequency analysis of these images can provide valuable information to enhance the predictive power of the radiomics-based model for classification of HPV+/- in OPSCC patients.

In this pilot study, we investigated the potential predictive value and interpretability of the three-dimensional (3D) Discrete Wavelet Decomposition (3D-WD) analysis for HPV characterization of GTVs of OPSCC patients. The primary objective of this research was to characterize and identify the spatial frequency sub-bands and their respective regions of primary tumor most likely responsible for characterization of HPV.

Workflow for Quadrature Mirror Filter (QMF) bank for analysis and synthesis components of the discrete WD decomposition and reconstruction is shown in subfigure 1A. Subfigure 1A shows the main structure of a two-channel QMF bank which consists of analysis and synthesis components. Subfigure 1B shows the ‘mother’ wavelet ( ψ) and its scaling functions (‘father’ wavelet, φ) in the time domain along with their four corresponding discrete analysis filters ( H( z)) and synthesis filters ( G( z)) generated for a Symlet function with two vanishing momentums. As shown in subfigure 1B, the Symlet function (with two momentums) is nearly symmetric and meets the orthogonality and biorthogonality properties.

Subfigure 2A shows one level of decomposition for the 3D discrete wavelet decomposition analysis used in this work. As shown in subfigures 2A and B, the 3D GTV information is successively decomposed using the basic operator block, the two-channel filter bank (QMF, see Subfigure 1A). At each decomposition level, seven Wavelet Coefficients (WCs, details) plus one approximation Wavelet Coefficient were generated.

The GTVs of the patients were decomposed at three levels. As shown in subfigure 2A and B, the 3D data is decomposed into a series of frequency-based details and approximation coefficients (corresponding to 22 WCs, 21 details plus one approximation). The first filtering operation at the end of each decomposition level, LLL, contains the low-pass component (highest approximation of the information) at that particular decomposition level. The rest of the filtering operations, ‘HLL’, ‘LHL’, ‘HHL’, ‘LLH’, ‘HLH’, ‘LHH’, and ‘HHH’ refer to the detailed coefficients (higher frequency sub-bands). The sequence of letters in each filtering operation gives the order in which the separable filtering operations are applied from left to right. For instance, ‘LHH’ means that the low-pass filter (scaling) with down-sampling (order of two) is applied to the rows (X dimension) of the GTV data, followed by the high-pass filter (wavelet) with down-sampling (order two) applied to the columns (Y dimension) of the GTV data. Finally, the high-pass filter with down-sampling (order of two) is applied to the 3 ^rd dimension (Z) of the GTV information. Subfigure 2B demonstrates the workflow and processing stages of the 3D-WD analysis, feature extraction, feature engineering, model development, and frequency-domain interpretation analysis used in the study. Among 22 LEE feature set corresponding to 22 different frequency sub-bands from three decomposition levels, only 13 LEE features showed high predictive power and were discriminant between the two patient cohorts (HPV- and HPV+).

Subfigure 3A shows a 10-fold nested cross validation with its inner and outer loop cores that was used in the model development and validation parts of the study. During each iteration, in the inner loop, 116 subjects were used as the training plus validation set, while 12 cases were used as the test set. Subfigure 3B illustrates average of the Receiver operating characteristic (ROC) curves corresponding to 10 folds nested CV for the inner loop, with each iteration related to a different division of the training plus validation sub-cohort for the Lasso-GLM analysis. Subfigure 3C demonstrates Receiver operating characteristic (ROC) curves corresponding to 500 iterations (using RPS), with each iteration related to a different division of training and testing for the optimal (λ=0.0088) Lasso-GLM classifier.

Table 1 illustrates the average values for AUC, PPV, NPV, and their corresponding values at 95% Confidence Intervals as well as their standard deviations and standard errors for the outer loop of the nested cross validation (NCV, N=10) process.

Subfigure 4A shows cross-validation analysis (K=10) results for finding the optimal L1-regularization parameter (λ at one standard error, Lambda 1SE). As shown in this figure, the average value of the λ at 1SE over ten different folds was 0.009. Subfigure 4B illustrates the Lasso-GLM feature ranking of the selected sub-band features in the Wavelet feature space at the optimal L ₁-regularization (λ=0.009). Subfigure 4C shows the trace plot of Lasso-GLM selected coefficients for the LEE based Wavelet features after L ₁-regularization. Trace plots with negative and positive coefficients are associated with HPV- and + classes respectively. As shown in this figure, the LEEs of the two WCs (HHL1 and HLL1 corresponding to details information) at decomposition level of one were associated with positive Lasso coefficients (shown in yellow and purple colors respectively, corresponding to HPV+ cohorts). Similarly, the LEEs of the two WCs (one approximation, LLL3, and one detail information, LHL3) at decomposition level of three were associated with negative Lasso coefficients (in blue, and brown colors respectively, corresponding to HPV- cohorts).

Figure 5 demonstrates normalized contrast enhanced CT images (first and third columns) after isotropic resampling, quantization, and re-segmentation of the data along with respective HPV probability maps (second and forth columns) for the segmented GTV cross sections of 10 example patients reconstructed by the inverse 3D-discrete WD analysis (see steps 7 to 11 of the Subfigure 2B).

As shown in Table 1, the NCV-based performance of the Lasso-GLM classifier for the prediction of HPV status was: 0.800, 0.800, and 0.910 for AUC, PPV, and NPV, respectively. The CV-based evaluation (with 500 RPS) of the optimal Lasso-GLM (with average L ₁ regularization parameter) performance were: AUC/PPV/NPV= 0.798/0.745/0.823.

A frequency-based spectral analysis was introduced to decompose the intensity and textural information of the tumor into different frequency sub-bands to reveal the frequency-based information most associated with characteristics of HPV. Figure 4B and 4C demonstrate spatial frequency sub-bands most associated with the characterization of HPV. Inclusion of such spatial frequency information in the development of radiomics-based models for prediction of HPV is likely to enhance model performance. Our group, among others ^{36 , 95 – 98} have studied the recruitment of spatial frequency-based information to extract image characteristics relevant to the target/outcome of interest. We have shown that frequency-based data can contain more relevant information compared to real-domain data. For instance, Martinez et al. ⁶⁷ utilized the Fourier frequency-spectral analysis method to evaluate different machine learning systems for automatic categorization of ovarian tumors on ultrasound images. Lofstedt et al. ⁹⁵ investigated GLCM textural analysis for generating textural features that are invariant from image pre-processing steps ( e.g. image quantization).

Other groups ^{95 – 98} have also recruited Fourier spectral analyses to find the optimal scale selection method for refinement of the GLCM textural analysis as well as for other tasks in image-based classification.

Despite the promise of the applications of radiomics analysis in the therapeutic and diagnosis fields for stratification and predictive modeling, many unanswered questions remain, as elaborated previously by the IBSI. ^{6 , 76} IBSI’s protocols for data processing are considered an important step toward standardization of algorithms used for extraction of radiomics features from medical images. These procedures will ensure that consistency and reproducibility are achieved among radiomics implementations at different institutions. The next step involves association between radiomics features with known patient outcomes such that reliable connections can be developed between significant radiomic features and outcomes for use in clinical modeling. These factors are considered key elements in the field of radiomics and much further investigation is needed to move this field forward. ^{75 , 76} This research provides an important step toward understanding of frequency-based parameters of the radiomics feature extractor cores based on the available information content and frequency spectra of image structures. Frequency-based features have the potential to improve association between radiomics features and outcomes by minimizing the chance of under/over-sampling of the image dataset during the feature extraction phase.

The conventional feature engineering processes along with multicollinearity analysis significantly reduced the number of correlated features from 22 to 13, suggesting a benefit of statistical-based feature engineering analysis prior to modeling to remove correlated features. Indeed, the feature engineering analysis results confirmed that suppression of the multicollinearity effects among different spatial frequencies prior to the Lasso analysis can improve model performance and enable more robust L ₁ regularization in the frequency-based feature space. We used L ₁ regularization to control the overfitting problem by penalizing the model coefficients in Lasso optimization. However, the Lasso technique is limited. For instance, in a group of variables where pairwise correlations are high (or where variables have high level of multicollinearity) Lasso will tend to select only arbitrarily and randomly one feature as the representative of the entire group of features. To address this limitation, we utilized the Belsley collinearity diagnostics test ⁹⁰ to remove the frequency-based features with high level of multicollinearity. This decreases the feature dimensionality, which improves model accuracy and generalizability. However, the Belsley test can still be limited in the small sample size setting which is the case in this study (N = 128).

In addition to the multicollinearity effect, the Lasso-GLM classifier may be prone to overfitting due to the small sample size (N = 128). The Lasso method is also sensitive to the chosen value of the tuning parameter (penalty term) which constrains the magnitude of the estimated coefficients. ^{91 – 93} The tuning parameter controls the amount of regularization, and selecting the optimal value is an important step since each parameter value corresponds to a different fitted model. ^{91 – 93} To address this issue, we utilized K-fold Cross validation (KFCV) in the inner loop of the NCV with collective deviance. ^{91 – 93} To minimize the effects of potential feature non-uniformities from different patient sub cohorts/folds, we used a higher number of folds (K = 10) in the KFCV and NCV analyses compared to the default number (K=5) recommended by the literature. ^{91 – 93}

As shown in Figure 3C and Table 1, the high variance of the prediction performance in the ROC analysis of the optimal Lasso-GLM classifier illustrates the contributions of different frequency bandwidths with insignificant information at an average L ₁-regularization toward characterization of HPV. The ROC family curves (corresponding to 10 folds) show a tight spread (see Figure 3B) suggests higher generalizability of the Lasso-GLM from contribution of distinct frequency bands at different L ₁-regularization levels. Such a tight spread implies that the dominant frequency features used to develop the model are sparse and less sensitive to tumor heterogeneity across different patients.

Model development, validation, and testing without nested CV uses the same data to train the classifier parameters and evaluate model performance. This could lead to information leakage from the training cohort into the testing cohort resulting in overfitting of the model. The overfitting magnitude is primarily dependent on the size of the dataset and the complexity of the model. To minimize this problem, we implemented a nested CV technique ⁹⁴ with a conventional KFCV in its inner loop, in which the testing dataset is held-out during the training and validation phases at each of multiple folds.

Prior to the wavelet analysis, the time-domain spectral information of the GTVs was interpolated to the reference dimensions. Such interpolation can increase the chance of additional up-sampling of the tumor information/resolution, adding uncertainty to the dominant frequency sub-bands identified as the discriminant features for the characterization of HPV. Finding an optimal interpolation technique with minimal effect on the tumor resolution/information is a challenging task and requires further investigation. ⁹⁹

Image artifacts have potential to impact on the wavelet analysis results and consequently impact the detection of the dominant frequency sub-bands. We inspected CT images of all 128 patients for identification of streaking artifacts impacting the visibility of the GTV’s. Given that the number of CT slices with dental filling artifacts was less than about 1% of the total number of 2D-CT cuts utilized (N=3596), we elected not to split the dataset into separate cohorts (with and without CT artifacts) for separate wavelet analyses and relevant modeling.

The wavelet decomposition technique is a mathematical transformation which decomposes a signal into multiple lower resolution levels by controlling the scaling and shifting factors of a single wavelet function (‘mother’ wavelet). Unlike the discrete Fourier transform (DFT), the wavelet analysis refers not just to a single transform, but rather a set of transforms, each with a different set of wavelet basis functions. Robustness and efficiency of different frequency-based methods are highly dependent on the properties of their transformation kernels. Here, a Symlet wavelet function (with two momentums) which was nearly symmetrical, orthogonal, and biorthogonal was used as the ‘mother’ wavelet of the analysis. The symlets are nearly symmetrical wavelets proposed by Daubechies as modifications to the Daubechies wavelet family. ^{84 – 86}

Nonsymmetrical filters can cause a shift variance on the frequency spectrum which can introduce a biasing effect on the analysis results. It has been shown ^{100 , 101} that the Symlet functions which belongs to the orthogonal wavelet family are compactly supported wavelets with minimal asymmetry and a maximum number of vanishing moments for a given support length. ^{100 , 101} One of the main advantages of using bi-orthogonal wavelet (such as Symlets) compared to orthogonal wavelet is their ability to filter phase information linearly. This means that the phase response of the filter is a linear function of the frequency which ensures preservation of energy of the signal during the transformation.

Indeed, in addition to the Symlets, there are many other ‘mother’ wavelet functions (such as Haar, Daubechies, Coiflets, Morlet, Mexican Hat, Meyer, and other real and complex wavelets) which can be used in the wavelet analysis. However, each of these basis functions has different properties and can produce different results when they are applied on the same dataset. Thus, the “mother’ wavelet function and its properties (such as supporting length, number of vanishing momentums, etc.) can be selected and fine-tuned according to the image specifications and the goals of the study. Thus, for a more accurate analysis of different image resolutions, these parameters should be carefully selected and optimized accordingly. For instance, the wavelet kernel with the shortest possible support length of one (such as the Haar wavelet) can capture more detailed information enabling a more precise linear transformation to generate the high-pass coefficients. Thus, it might be more beneficial to use an asymmetrical frequency-based ‘mother’ wavelet for the analysis of images in which axial sampling pitch ( i.e. slice thickness, intra slice distance) is much larger than in-plane sampling pitch (inter-slice distance), as it is often the case with CT and MR modalities. Of note, in our analyses, we used a nearly symmetrical ‘mother’ wavelet because the wavelet analysis was performed on a volumetric dataset which was isotropically resampled to 1 mm ³/voxel.

The results (see Figure 5) of this pilot study confirm that compared to the central-zones of tumors, peritumoral-regions have more information for characterization of HPV. The 3D-wavelet analysis of GTVs in these patients not only provided adequate information for building an HPV classifier, but also revealed important spatial frequency-based tumor information which can be used as the most representative features for HPV classification. The HPV association probability maps can be susceptible to the radiomics feature Log Energy Entropy (LEE) which was used to evaluate the importance of the wavelet coefficients as well as the properties of the ‘mother’ wavelet function. The number of vanishing moments, shape (symmetric or asymmetric), scaling function, minimum supporting length, orthogonality or biorthogonality conditions, etc. can impact on the reconstruction error of the inverse wavelet analysis, directly affecting the identified tumor zones associated with the HPV. Further investigation in this area of research is warranted.

The frequency decomposition analysis results revealed that any spatial frequency information of the tumor with up-sampling level higher than three (2 ^L with L = 3, equivalent to eight times of the image resolution, 8 mm) does not add any significant HPV-related information to the extracted features set. Results also imply that tumor zones with smaller image structures or components corresponding to higher details or higher spatial frequencies (components ~<2 mm corresponding to one decomposition level) contain important radiomics information for characterization of HPV+ patients. On the other hand, tumor zones with larger image structures or components corresponding to lower spatial frequencies (components ~>8 mm, corresponding to three decomposition levels) contain more information for classification of HPV- cases.

Our group and others ^{53 , 102 – 106} have previously shown that in OPSCC patients, as the mean breadth value (proportional to tumor size) of the tumor increases the chance the tumor for being classified as HPV negative increases. Interestingly, in this study, the LLL3 frequency sub-band (spatial frequency information approximation with three levels of up-sampling being more related to the overall structure of tumor) being identified as a discriminant feature associated with the HPV negativity agrees with this finding.

Further investigation, including adaptive filtering analysis (designing ‘mother’ wavelet filters according to the image specifications), using different radiomics features and related impact on reconstruction errors in images with different levels of frequency information are warranted.

The findings of this pilot study and similar studies can be beneficial toward establishing the role of filter-based radiomic analysis in studying the underlying imaging characteristics of clinical relevance. The main aim of this study was to apply a series of feature engineering processes in the frequency domain to suppress the frequency bandwidth information of the tumor less relevant to HPV characterization while enhancing the tumor frequency bandwidths information descriptive towards prediction of HPV.

Another limitation of this work pertains to the use of p-16 immunohistochemistry (IHC) testing for determining the HPV status of the OPSCC patients. Although it has become a standard-of-care method for HPV testing of OPSCC, the p-16 IHC testing method does have limitations which can impact the accuracy of HPV test results. ^{107 – 110} Recent systematic reviews of the diagnostic accuracy of p16-IHC in OPSCCs analyzed 24 and 31 different studies and concluded that p16-IHC has a high sensitivity (90.5%), but modest specificity (83.3%) to classify HPV status when used as a single test. ^{109 , 111 , 112} Such variations in the performance of the p16-IHC testing can directly affect the source of truth, which was used for the feature engineering and model development parts of the study. This would also add a degree of uncertainty to the model accuracy. Despite this limitation, the radiomics findings here agree with those of other studies and with radiological findings of tumor characteristics of HPV+/- patients with OPSCC’s. ¹¹³

We developed a frequency-based wavelet analysis to characterize and identify the spatial frequency sub-bands and regions of primary tumor responsible for the characterization of HPV in patients with OPSCC. The overall findings of this study, albeit subject to confirmation in a larger patient population, demonstrated encouraging results in support of frequency analysis of CT images toward tumor zone-based interpretation, characterization, radiomics-based modeling of HPV and differentiation between HPV+/- cohorts in this patient cohort.

This project contains the following underlying data which can be directly downloaded from GITHUB using the following link:

https://github.com/Ebadian-HFHS/F1000-DATASETS 1-

HPVClass_128_Cases.txt: The tumor biopsies immunohistochemistry P16 protein (IHC) testing results for all128 OPSCC patients.