This was a retrospective, multicenter, and diagnostic study. From January 2019 to December 2019, 153 patients, which were collected consecutively from Yantai Yuhuangding Hospital, were randomly assigned to a training set or an internal test set at a ratio of 7:3 ( Figure 1A ). Additionally, 46 consecutive patients at Qilu Hospital of Shandong University were recruited to serve as an external test set from June 2019 to December 2019 ( Figure 1B ). The inclusion and exclusion criteria can be found in Appendices .

The protocol for this retrospective study was approved by the Institutional Review Board of The Affiliated Yantai Yuhuangding Hospital of Qingdao University. The requirement for obtaining written informed consent from the patients was waived by the review board. A study work flow diagram is shown in Figure 1C .

CT assessments of ETE were performed by two experienced radiologists, one with 15 years of experience and another with 5 years of experience. Both radiologists were blinded to histopathological results. Any disagreement was resolved by discussion or consultation with a third radiologist who had 20 years of experience. ETE was reported on CT images when at least one of the following CT criteria was fulfilled: 1) tumor in contact with 180° or more of the tracheal, esophageal, muscle or vascular circumference; 2) loss of normal structure; 3) a clinical symptom (e.g., ipsilateral vocal cord palsy) was present that could be explained by the CT images; 4) tumor showed poorly defined margin with heterogeneous signal intensity in adjacent soft tissue; 5) focal bulging out or disruption of the thyroid capsule by tumor; 6) more than 25% of perimeter of the tumor was abutting the thyroid capsule (10). Kendall’s coefficient of concordance W (Kendall’s W) was used to evaluate inter-radiologist agreement. W scores range between 0 and 1, which was graded as very good (0.80 to 1.00), good (0.60 to 0.80), fair (0.40 to 0.60), moderate (0.20 to 0.40) or poor (<0.20) (14).

The Image Acquisition protocol is shown in Appendices . After that, a radiologist with 8 years of experience manually delineated the PTC regions of interest (ROI) (15) on non-contrast and venous phase CT images by using an ITK-SNAP (version 3.8.0; www.itksnap.org). Based on their traits, ROIs were manually delineated along the tumor contour on each transverse section and the ROIs in each slice constituted a volume of interest (VOI). A sample of segmentation process is shown in Figures 1D, E . Three months later, another two radiologists segmented the images of 30 patients who were randomly selected to assess the intra- and inter-observer reproducibility of radiomics features by intra- and inter-class correlation coefficients (ICCs). An ICC > 0.8 indicated an excellent agreement of radiomics features.

Before feature extraction, VOIs were prepossessed included gray value standardization, gray level discretization, and image resampling. A total of 2818 radiomics features were originally extracted from each VOI by using PyRadiomics (16) on Python (version 3.7). A total of 1409 sequential features based on time series were calculated as Delta _{(radiomics feature)}, which was defined as the difference between features seen on the venous phase and non-contrast CT images. Delta _{(radiomics feature)} was calculated as follows:

A total of 4227 features (containing radiomics features and sequential features) were identified for each patient containing 19 types of first order statistics, 13 types of Shape features, 28 types of Gray-Level Co-occurrence Matrix (GLCM) features, 16 types of Gray-Level Run-Length Matrix (GLRLM) features, 16 types of Gray Level SizeZone Matrix (GLSZM) features, and 56 types of sequential features.

The values of the extracted features were standardized with z scores by using mean and standard deviation values. Next a two-step procedure for dimensionality reduction and feature selection was devised. First, an analysis of variance (ANOVA) was performed to screen out discriminative features in the training set, with only features having a P < 0.05 being selected. Next, the least absolute shrinkage and selection operator (LASSO) method was used to reduce the dimensions of features, identify the most significant features (17), and make a second selection.

K-Nearest Neighbor (KNN), Logistics Regression (LR), Decision Tree (DT), Linear-support vector machine (Linear-SVM), Gaussian-SVM, and Polynomial-SVM methods were used to construct radiomics models, respectively. Area under the curve (AUC) was regarded as a performance indicator and used to evaluate the performance of the radiomics-based models built by using each classifier. Finally, a radiomics signature (Rad-score) was calculated by using a linear combination of selected features weighted by the corresponding LASSO coefficients.

In the training set, one-way ANOVA and multivariate logistic regression analysis were used to screen clinical risk factors including sex, age, primary site, tumor diameter, and radiologist interpretation to identify the clinical risk factors for ETE. To quantify the roles of the Rad-score and clinical risk factors in predicting ETE, both of those risk factors were included in a multivariate logistic regression analysis that was conducted using a backward-stepwise approach, where collinearity was considered and risk factors with a variance inflation factor (VIF) > 10 and a P > 0.05 were excluded. When the minimum Akaike information criterion was reached, the Akaike information was taken as the criterion, the stepwise procedure was stopped, and the final multivariate logistic regression constituted the nomogram.

The nomogram’s predictive performance was evaluated in the training set, internal test set, and external test set, respectively. The nomogram’s performance was evaluated in terms of the receiver operating characteristic (ROC) curve. In particular, the comparison between the nomogram and radiologists’ interpretation was illustrated in the same ROC curve. A decision curve analysis (DCA) was plotted for the entire set and applied in a clinical usefulness evaluation of the nomogram by calculating the net benefits at different threshold values in the training and test sets. In addition, a calibration (agreement between observations and prediction of ETE) curve was used to evaluate the agreement between actual status and ETE probabilities as predicted by the nomogram, accompanied by the Hosmer-Lemeshow test.

All statistical analyses were performed using Python (version 3.6, https://www.python.org), R software (version 4.0.3, https://www.r-project.org), and SPSS (version 26.0, IBM Corp.). Scikit-Learn, a Python library, was used for selecting radiomics features and constructing a radiomics model. The modules of “feature-selection,” “linear-model,” “svm,” “neighbors,” “tree,” and “metrics” were used for the entire procedure. The “rms” package in R software was used to select clinical risk factors, build the nomogram, and plot calibration curves. The “rmda” package was used to perform DCA. SPSS software was used to compare categorical variables (i.e., sex, T stage, primary site, radiologists’ interpretation, and lymph node metastasis) by using the χ² test or Fisher’s exact test. Continuous variables were compared by using the student’s t test or Mann-Whitney U test, when appropriate. Statistical significance was two-sided, and a P < 0.05 was considered to be statistically significant.

The numbers of enrolled patients in the training set, internal test set, and external test set were 107, 46, and 46, respectively. There were no significant differences in sex, age, diameter, Thyroid Stimulating Hormone (TSH), Free Thyroxine (FT4), Free Triiodothyronine (FT3), and Primary site in each data set, but there was a significant difference in T stage. Details of clinical characteristics in the 3 sets are shown in Table 1 . Radiologist A and radiologist B reported all CT imaging results with good inter-radiologist agreement (Kendall’s W, 0.823).

A total of 4227 radiomics features were extracted from each patient’s images. The intra-observer and inter-observer ICCs were 0.867 and 0.915, respectively, indicating an excellent agreement on radiomics features. Four radiomics features (Small Area Low Gray Level Emphasis, Gray Level Variance, Difference Entropy, and Busyness) in the training set were finally selected ( Figure 2 ; Table 2 ).

Predictive models built by using KNN, LR, DT, Linear-SVM, Gaussian-SVM, and Polynomial-SVM, showed favorable discriminatory abilities for radiomics parameters with AUCs of 0.841 (95% Confidence Interval [CI], 0.769-0.900), 0.774 (95% CI, 0.683-0.853), 0.841 (95% CI, 0.772-0.897), 0.774 (95% CI, 0.680-0.854), 0.805 (95% CI, 0.710-0.879), and 0.815 (0.708-0.888), respectively, in the training set (Figure 3 A). In the internal test set, the six radiomics models also showed good performance with AUCs of 0.669 (95% Confidence Interval [CI], 0.394-0.568), 0.701 (95% CI, 0.480-0.767), 0.680 (95% CI, 0.516-0.795), 0.694 (95% CI, 0.508-0.788), 0.642 (95% CI, 0.454-0.600), and 0.681 (0.466-0.700), respectively ( Figure 3B ). Interestingly, in terms of AUC, the LR, which was later used to build the nomogram, demonstrated the best performance of the 6 models.

The 4 features mentioned above were selected and their coefficients were plugged into a formula for calculating Rad-score as described below:

The ANOVA and multivariate logistic regression analysis only identified radiologists’ interpretation as being an independent ETE predictor ( Table 3 ). The values of Rad-score and radiologists’ interpretation for predicting the presence of ETE in PTC patients were quantified in our nomogram (score for each factor is shown in Figure 4 ). The radiomics nomogram showed a good prognostic capability, with AUCs of 0.860 (95% CI, 0.790-0.931), 0.750 (95% CI, 0.579-0.921), and 0.797 (95% CI, 0.665-0.926) in the training, internal test, and external test sets, respectively ( Figures 3C–E ; Table 4 ). Furthermore, the radiomics nomogram had a more powerful predictive efficiency than the interpretation of radiologists with 5 and 10 years of experience, whose AUC in the training, internal test, and external test sets were 0.762 (95% CI, 0.667-0.858), 0.718 (95% CI, 0.560-0.876), and 0.733 (95% CI, 0.615-0.841), respectively ( Figures 3C–E ). In the external test set, the nomogram’s accuracy (0.848 vs. 0.717) and specificity (0.909 vs, 0.697) were higher than those of radiologist’s interpretation, reflecting a better prognostic performance of the radiomics nomogram; however, the sensitivity of the nomogram (0.692, 95% CI 0.389-0.896) was similar with that of radiologists’ interpretation (0.769, 95% CI 0.460-0.938) ( Table 4 ). We tested this nomogram to distinguish minimal and grossly ETE in internal and external test set, and achieved relatively good performance with AUC of 0.800 (0.578-1.000) and 0.722(0.390-1.000), sensitivity of 0.800 (0.299-0.989) and 0.500(0.092-0.908), specificity of 0.667 (0.412-0.856) and 0.889(0.507-0.994), and accuracy of 0.696 (0.471-0.868) and 0.769(0.462-0.950) respectively. ( Figure 3F ). The DCA showed that radiomics nomogram could add more benefit than treat patients as all-ETE or non-ETE when the threshold probability ranged from 0.10 to 0.50 in internal test set and 0.12 to 0.75 in the external test set ( Figures 5A, B ). The calibration curve for the radiomics nomogram, when assessed by the Hosmer-Lemeshow test, yielded P of 0.48, 0.88, and 0.34 in the training set, internal test set, and external test set, respectively ( Figures 5C, D ), demonstrated good agreement between actual status and ETE probabilities as predicted by the nomogram in all three sets.