Development and validation of a radiomic signature to predict HPV (p16) status from standard CT imaging: a multicenter study.

OBJECTIVES: Human papillomavirus (HPV) positive oropharyngeal cancer (oropharyngeal squamous cell carcinoma, OPSCC) is biologically and clinically different from HPV negative OPSCC. Here, we evaluate the use of a radiomic approach to identify the HPV status of OPSCC.
METHODS: Four independent cohorts, totaling 778 OPSCC patients with HPV determined by p16 were collected. We randomly assigned 80% of all data for model training (N = 628) and 20% for validation (N = 150). On the pre-treatment CT images, 902 radiomic features were calculated from the gross tumor volume. Multivariable modeling was performed using least absolute shrinkage and selection operator. To assess the impact of CT artifacts in predicting HPV (p16), a model was developed on all training data (Mall) and on the artifact-free subset of training data (Mno art). Models were validated on all validation data (Vall), and the subgroups with (Vart) and without (Vno art) artifacts. Kaplan-Meier survival analysis was performed to compare HPV status based on p16 and radiomic model predictions.
RESULTS: The area under the receiver operator curve for Mall and Mno art ranged between 0.70 and 0.80 and was not significantly different for all validation data sets. There was a consistent and significant split between survival curves with HPV status determined by p16 [p = 0.007; hazard ratio (HR): 0.46], Mall (p = 0.036; HR: 0.55) and Mno art (p = 0.027; HR: 0.49).
CONCLUSION: This study provides proof of concept that molecular information can be derived from standard medical images and shows potential for radiomics as imaging biomarker of HPV status. Advances in knowledge: Radiomics has the potential to identify clinically relevant molecular phenotypes.

RCT Entities:   Population Interventions Outcomes

Introduction

Over the last years, the incidence of oropharyngeal squamous cell carcinoma  (OPSCC) has shown a dramatic increase relative to other head and neck cancers, with a substantial proportion of OPSCC being linked to human papillomavirus (HPV) infections.[1]

HPV positive OPSCC is biologically and clinically different from HPV negative OPSCC, which is often related to alcohol and tobacco consumption. HPV positive OPSCC has been shown to have superior response to radiochemotherapy. Approximately, 80% of HPV positive OPSCC patients achieve locoregional control and 5 years overall survival, in comparison to less than 50% of patients with HPV negative OPSCC and non-oropharyngeal head and neck cancers.[2, 3] This favorable outcome makes HPV positive OPSCC in particular interesting for de-escalation protocols.[4]

Widely accepted methods for detection of HPV infection are in situ hybridization for viral DNA, HPV DNA or RNA PCR, and immunohistochemical investigation of the level of p16 expression, which strongly correlates with HPV infection.[5]

Radiomics is a rapidly emerging field, introduced in 2012, which concerns with the high-throughput mining of large amounts of quantitative features, derived from (standard-of-care) medical imaging, for knowledge extraction.[6-9] Radiomics is in particular promising within decision support systems for precision medicine[10-12] and its potential to predict HPV status in head and neck cancer has been recognized.[13] Indeed, previous studies have reported radiologic differences between HPV positive and negative cases in terms of qualitative radiologist’s readout[14] or perfusion CT.[15] Furthermore, exploratory radiomic studies have shown that heterogeneity of image-based density is potentially associated with HPV in OPSCC.[16, 17]

Most of the studies investigating imaging phenotypes of tumors are based on single center data, which introduce bias to a model and limits its applicability.[8] In particular, factors such as CT scanner, tube voltage, tube current, reconstruction kernel and contrast agent influence the results of quantitative analysis. In this multicenter study, we further investigate if a quantitative CT-based radiomic approach can objectively identify the HPV (p16) status of OPSCC, by developing and validating a radiomic signature on a large and international collection of patient data from four different institutions. This study does not intend to develop methodology to replace existing HPV tests, yet aims to provide a proof of concept that radiomics is able to derive molecular information from standard medical images.

Methods and materials

Patients and CT imaging

Four independent cohorts, with a total of 778 OPSCC patients with HPV status determined by p16 immunohistochemistry and treated with curative intent by radiation therapy with/without concurrent chemotherapy, were collected from the Princess Margaret Cancer Center (N = 427), the VU University Medical Center (N = 158), the University Hospital Zürich (N = 100) and MAASTRO clinic (N = 93). All patients underwent pre-treatment contrast enhanced CT imaging of the head and neck. The gross primary tumor volume (GTV) was manually delineated for each patient for treatment planning purposes. The delineation was not standardized and was performed according to the clinical protocols, separate for each institute. Images were visually assessed for the presence of CT artifacts (e.g. streak artifacts due to dental fillings) within the GTV. A more detailed description of acquired CT images for each of the included cohorts can be found in the Supplementary Material 1. Institutional review board approval was obtained for each of the participating centers. Patients provided informed written consent, unless the need for written consent for this retrospective study was waived by the participating center.

Image analysis

Prior to analysis, all images were resampled to isotropic voxels of 2 mm, using linear interpolation.[18] A total of 902 radiomic features were calculated, divided into five groups: tumor intensity, shape, texture, Wavelet and Laplacian of Gaussian. All features were extracted using in-house developed software, using Matlab 2014a (MathWorks, Natick, MA). Feature descriptions and mathematical definitions can be found elsewhere.[8, 19] To calculate wavelet features, we used the low pass approximation and the high pass decomposition (i.e. applying either a low or high pass filter in each direction, respectively), since these decompositions are directionally invariant. For Laplacian of Gaussian features, the texture size (fine to coarse) was highlighted by modifying the Gaussian radius parameter from 2 to 7 mm with 1 mm increments. Textural features were computed discretizing image intensities into bins, using both a bin width of 10 and 25 Hounsfield unit.[20]

Statistical analysis

We randomly assigned 80% of all data for model training (N = 628) and 20% for validation (N = 150), with balanced HPV status, institution, and number of patients with visible CT artifacts.

Highly correlated features were first removed from further analysis by evaluating all pairwise correlations in the training data set. For each highly correlated feature pair (Pearson correlation coefficient ρ > 0.9), the variable with the largest mean absolute correlation with all remaining features was removed.

Multivariable logistic regression was performed using the least absolute shrinkage and selection operator model selection technique,[21] with 100 times repeated 10-fold cross-validation to select the optimal tuning parameter (λ). To further reduce the chance of overfitting on the training data, we selected the simplest candidate model (i.e. the model with fewest non-zero coefficients), i.e. within one standard error of the best performing model. The area (AUC) under the receiver operator curve (ROC) was used to assess model performance in predicting HPV (p16) status.

Finally, we compared Kaplan–Meier survival curves between patients with positive and negative HPV status, based on conventional p16 immunohistochemistry and based on radiomic model HPV predictions, for all validation patients. Model class predictions were made with a probability cutoff of 0.5. Overall survival was defined as the time from start of treatment to death as a result of any cause. A log-rank test was applied to test for significant differences between survival curves.

To assess the impact of CT artifacts, a model was also developed on the subset of patients in the training cohort for which there were no visible CT artifacts within the GTV. All model validation was subsequently performed on the entire validation data (Vall), and the subgroups of validation patients with (Vart) and without (Vno art) CT artifacts. AUC values for paired ROC curves were compared using DeLong’s test.[22] Model calibration was measured by the intercept and slope of the logistic calibration curve.[23]

To further compare the two models, confusion matrices for HPV (p16) predictions by Mall and Mno art on all validation data Vall were determined.

Statistical analysis was performed in R (R Foundation for Statistical Computing;v. 3.3.3).

Results

Radiomic models for HPV prediction

Patient characteristics, including HPV status, presence of CT artifacts, and follow-up time are summarized in Table 1.

Table 1.

HPV status, presence of CT artifacts and median follow-up time for the PMH, the VUmc, the USZ, MAASTRO clinic (MAASTRO), the 80% training data (training) and the 20% validation data (validation)

Variable	PMH (n = 427)	VUmc (n = 158)	USZ (n = 100)	MAASTRO (n = 93)	Training (n = 628)	Validation (n = 150)
HPV (p16) status
Positive	303 (71%)	34 (22%)	56 (56%)	33 (35%)	344 (55%)	82 (55%)
Negative	124 (29%)	124 (78%)	44 (44%)	60 (65%)	284 (45%)	68 (45%)
CT artifacts
Yes	219 (51%)	69 (44%)	57 (57%)	26 (28%)	300 (48%)	71 (47%)
No	208 (49%)	89 (56%)	43 (43%)	67 (72%)	328 (52%)	79 (53%)
Overall survival
Median follow-up (months)	71.6	74	44.5	51.8	69.4	65.1
Cohort
PMH	–	–	–	–	343 (55%)	84 (56%)
VUmc	–	–	–	–	128 (20%)	30 (20%)
USZ	–	–	–	–	82 (13%)	18 (12%)
MAASTRO	–	–	–	–	75 (12%)	18 (12%)

PMH, Princess Margaret Cancer Center; UHZ, University Hospital Zürich; VUmc, VU University Medical Center.

For the training and validation data sets, the amount of patients from each individual cohort are given as well. Median follow-up for overall survival was determined by “reverse” Kaplan–Meier analysis (i.e. inversed censoring).

A subset of 165 uncorrelated features was preselected for the full training cohort, whereas a subset of 173 uncorrelated features was identified for the training data without CT artifacts. The models developed on all training data (Mall; 37 degrees of freedom) and on the subset of training data without CT artifacts (Mno art; 50 degrees of freedom). For model performance on the training data, we observed an AUC of 0.824 [95% CI (0.791–0.856)]. and 0.868 [95% CI (0.830–0.906)] for Mall and Mno art, respectively. Both models were subsequently validated on Vall, Vno art and Vart. The resulting AUC values, logistic calibration intercepts and slopes are summarized in Table 2. The corresponding ROC plots are shown in Figure 1.

Table 2.

AUC values, logistic calibration intercepts and slopes for the model developed on all training data (Mall) and the model developed on the subset of training patients without CT artifacts (Mno art), validated in all validation data (Vall), the subset of validation data without CT artifacts (Vno art) and the subset of validation data with CT artifacts (Vart)

Model	Validation dataset	AUC	Intercept	Slope
Mall	Vall	0.7636 [95% CI (0.6874–0.8399)]	0.034	1.041
Vno art	0.7658 [95% CI 0.6592–0.8724)]	−0.238	1.191
Vart	0.7521 [95% CI (0.6378–0.8665)]	0.37	0.852
Mno art	Vall	0.7391 [95% CI (0.6582–0.8199)]	0.408	0.561
Vno art	0.8005 [95% CI (0.6967–0.9044)]	0.057	1.103
Vart	0.7017 [95% CI (0.5775–0.8259)]	0.767	0.341

AUC, area under the curve.

Figure 1.

ROC plots for the model developed on all training data Mall and the model developed on the subset of training data without CT artifacts Mno art, validated on all validation data Vall (a), the subset of validation data without CT artifacts Vno art (b) and the subset of validation data with CT artifacts Vart (c). AUC, area under the curve; ROC, receiver operator curve.

Kaplan–Meier survival curves, including numbers at risk, for all validation data Vall are shown in Figure 2. For HPV determined by p16, there was a significant split between survival curves for HPV (p16) positive and negative cases (p = 0.007), with a hazard ratio of 0.46 [95% CI (0.26–0.82)]. For HPV (p16) predictions by Mall (p = 0.036) and Mno art (p = 0.027), we observed a similar significant split between survival curves, with hazard ratios of 0.55 [95% CI (0.31–0.97)] and 0.49 [95% CI (0.26–0.93)], respectively.

Figure 2.

Kaplan–Meier curves and number of patients at risk for HPV predictions by Mall vs p16 (a) and Mno art vs p16 (b). Survival times are in months. HPV, human papilloma virus.

Comparison of HPV models based on training data with and without CT artifacts

AUC values for HPV (p16) predictions made by Mall and Mno art were not significantly different for all validation data sets. Confusion matrices for HPV (p16) predictions by Mall and Mno art on all validation data Vall are shown in Tables 3 and 4.

Discussion

In this multicenter study, we developed and validated a CT based radiomic signature to predict the HPV status of OPSCC patients. In the context of radiogenomics,[24, 25] our study provides a proof of concept that molecular information can be inferred from standard medical images by means of radiomics.

Previous exploratory radiomic studies that indicated a correlation between HPV infection and heterogeneity of imaging-based tumor density in OPSCC[16, 17] either were performed on small populations without validation, or only used single institution data for both model development and validation. This is a major issue in radiomic studies, as can be learned from recent literature describing the process and challenges of radiomics.[8,13,26-28] In this multicenter study, we used a large collection of imaging data from four different institutions for model development and validation.

Previously published studies report that HPV positive tumors are more homogenous in CT density.[16, 17] Although a full interpretation of the complex radiomic signature for HPV is difficult, we have analyzed the meaning of several features selected both in the model trained on all training data and the model trained on the subset of training data with no CT artifacts. We have summarized our observations in Table 5. Our models selected different features than the ones previously published. However, as many radiomic features are correlated with each other, it is more relevant to compare features interpretation than name. Our multicenter results confirm that HPV tumors are more homogenous in CT density. This study provides also an additional insight into HPV imaging phenotype. We have observed that the HPV positive tumors seem to be characterized by lower contrast uptake, lower minimum density, and higher changes in the intensity of adjacent voxels.

Table 5.

Interpretation of selected features in the HPV radiomic signatures

	Radiomic feature	Interpretation
Decreased in the HPV tumors in comparison to HPV negative	Gray level size zone matrixSmall zone emphasis	Higher homogeneity
Gray level co-occurrence matrix Inverse variance	Higher changes in the intensity of adjacent voxels
Laplacian of Gaussian (4 mm)10th percentile	Lower minimum HU value
Increased in the HPV tumors in comparison to HPV negative	Gray level size zone matrixLow gray level large size emphasis	Lower contrast uptake
Laplacian of Gaussian (3 mm) kurtosis	More outliers

HPV, human papillomavirus; HU, Hounsfield unit.

Histopathology analysis shows differences between HPV+ and HPV- microscopic traits. For example, HPV+ tumors are characterized by lobular growth, infiltrating lymphocytes and well differentiated cells.[29] A direct link between these histopathology traits and radiomic signature of HPV positive tumors is not possible at this stage, and would require further investigation on a surgical cohort with full histology data.

Including data from different institutions introduces variety in image acquisition and reconstruction, which has been shown to affect radiomic features.[30, 31] Shafiq-ul-Hassan et al[18] investigated voxel-size dependency of radiomic features and found that the robustness of radiomic analyses can be improved by resampling to a nominal voxel size or by normalizing the voxel size. All images in this study were, therefore, resampled to isotropic voxels of 2 mm, which was approximately the average slice spacing, using linear interpolation. Furthermore, as shown previously, textural features and their interpretation are affected by the bin width used to discretize image intensities[20] . Therefore, features calculated for different bin widths may provide additional predictive information. To account for this, textural features were computed using both a bin width of 10 and 25 Hounsfield unit. In our HPV radiomic signatures, features calculated using both bin sizes were selected and no preference for a bin size was observed.

Besides variability in CT imaging, demographic differences also have to be considered. Developing a model on a single, independent cohort is, therefore, unlikely to sufficiently capture the variability that exists across datasets, resulting in a model with poor generalizability. We, therefore, performed our model development on more heterogeneous data, by randomly assigning 80% of all included data for model training and 20% for testing, with balanced HPV (p16) status, institution and number of patients with visible CT artifacts.

A common concern in the analysis of CT images of head and neck cancer are metallic dental fillings or other high atomic number material implants, which result in imaging artifacts.[32] An existing radiomic signature for overall survival[7] has previously been shown to have prognostic power regardless of CT artifacts.[33] Another recent study exploring the link between HPV status and CT radiomics, preprocessed images by completely removing artifacts affected slices from analysis.[17] However, such a process neglects potentially relevant three-dimensional information. To investigate the impact of CT artifacts on HPV prediction, we developed a model on all data (Mall) and a model on the subset of data without artifacts (Mno art). What can be observed from our results is that there is no significant difference in discriminative power of both models. However, overall calibration of Mall was better than that of Mno art. It has to be noted that the extent of CT artifacts and the impact on radiomic features will vary between patients. For an individual patient, model accuracy will, therefore most likely depend on the amount of the tumor region that is obscured by artifacts. This would have to be further investigated, preferably including techniques for metal artifact reduction in CT.

Since HPV-related OPSCC have been shown to have superior response to radio-chemotherapy,[2, 3] we compared Kaplan–Meier survival curves between patients with positive and negative HPV status, based on p16 and model class predictions by Mall and Mno art for all validation patients. Indeed, we observed a significant split (p < 0.05) between HPV positive and HPV negative patients based on p16. For HPV (p16) predictions based on both models, we obtained survival curves similar to that of p16, with significantly different survival for HPV positive and HPV negative patients, indicating that model predictions are indeed in line with p16.

It has previously been shown that part of the OPSCC patients that test positive for p16 immunohistochemistry are in fact HPV DNA negative.[34, 35] Since HPV testing for patients included in our study was performed by p16, the likelihood of false-positives has to be acknowledged. Furthermore, model class predictions (i.e. predicting either HPV positive or HPV negative), were made with a probability cut-off of 0.5, meaning that the costs for false-positives and false-negatives were considered equal. In clinical practice, this will not be the case and false-positives should be avoided (i.e. have a high cost), in order not to unjustly deescalate any patient’s treatment. To achieve a clinically acceptable level of accuracy, further development and validation would be needed, including HPV DNA. It is, therefore, important to note that radiomic HPV prediction models are not meant to replace HPV testing and we acknowledge that clinical decision making should always be made on the universally accepted most accurate testing (i.e. p16 positive tumors should be subjected to HPV DNA testing).

Considering HPV testing is routinely performed for OPSCC patients in western countries, the clinical usefulness of a radiomic biomarker could be considered to be limited. However, our results show there is potential for radiomics to serve as a cost-effective, complementary method for HPV screening, which may also be useful in non-oropharyngeal SCCs.[36] Another potential application for a reliable radiomic biomarker could be to perform (retrospective) HPV analyses when no tissue samples are available, or in countries where it is not routinely done.

In this study, we only considered the primary tumor. However, HPV-associated OPSCC, commonly present with a relatively smaller primary tumor, and relatively more advanced nodal disease. Severity of the disease may then be overestimated by the resulting higher tumour, node and metastasis and overall stage, as these are related to other head and neck SCC.[37] HPV has also been shown to affect the morphology of affected lymph nodes.[38] Including radiomics of involved lymph nodes could potentially provide additional value in predicting HPV status. Furthermore, additional improvement in inferring tumor HPV status may be achieved when combining radiomics with clinical features.[14]

Table 3.

Confusion matrix for HPV (p16) predictions by Mall on all validation data Vall

		Reference
HPV–	HPV+
Predictions	HPV–	44	19
HPV+	24	63

HPV, human papillomavirus.

Table 4.

Confusion matrix for HPV (p16) predictions by Mno art on all validation data Vall

		Reference
HPV–	HPV+
Predictions	HPV–	57	34
HPV+	11	48

HPV, human papillomavirus.