The superior value of radiomics to sonographic assessment for ultrasound-based evaluation of extrathyroidal extension in papillary thyroid carcinoma: a retrospective study

This study included 460 consecutive patients (575 nodules) who were diagnosed with PTC pathologically from January 2018 to September 2019 in our hospital. The inclusion criteria were as follows: 1) preoperative US examination within two weeks of surgery; 2) initial surgical resection (total thyroidectomy, near total thyroidectomy or hemithy-roidectomy) with cervical dissection, and the ETE state was confirmed pathologically and/or by intraoperative evidence. The exclusion criteria were as follows: 1) incomplete clinical characteristics; 2) multifocal lesions in the same lobe of thyroid; 3) radiofrequency ablation before surgery. Finally, 310 patients (337 nodules) with mean age of 45.95 ± 11.69 years (from 12 to 77 years) and male-to-female ratio of 0.40 (89:221) were included in this study. Among them, 27 patients had one lesion on each lobe of the thyroid. According to previous study, minimal ETE was defined as an extension into the sternothyroid muscle or parathyroid soft tissues.⁹ And it was confirmed pathologically in this study. Meanwhile, gross ETE was defined as cancer invasion into the subcutaneous soft tissue, trachea, larynx, esophagus or recurrent laryngeal nerve.⁹ It was determined during surgery and further confirmed pathologically. Thus, nodules were divided into those with minimal ETE (n = 99), those with gross ETE (n = 4) and those without ETE (n = 234) (ETE negative group). Because of the very small sample of gross ETE, nodules with minimal ETE and gross ETE were categorized as ETE positive group.

Clinical characteristics comprised of age, sex, biochemical results and US findings. Standardized biochemical examination included those concerning preoperative assessment and those concerning thyroid function. Liver function analysis, renal function analysis, differential blood count and routine urianlysis were classified into former ones. Blood calcium ion, thyroid function analysis (free triiodothyronine 3 and 4, total triiodothyronine 3 and 4 and thyroid stimulating hormone) (FT3, FT4, TT3, TT4 and TSH), thyroglobulin (TG), antithyroglobulin antibodies (ANTITGAB) and antithyroid peroxidase antibody (ANTITPOAB) were classified into latter ones. The composition, echogenicity, shape, margin and echogenic foci consisted of the US findings of each nodule. They were measured, scored and classified according to the Thyroid Imaging Reporting and Data System (TI-RADS) criteria of American College of Radiology.¹⁹ A total of 34 clinical characteristics were enrolled. The flowchart of this study was shown in Figure 1.

Figure 1.

The flowchart for this retrospective study.

AUC = area under the curve; ETE = extrathyroidal extension; PTC = papillary thyroid carcinoma; US = ultrasonic; xgboost = the extreme gradient boosting

The assessment of ETE was performed by two US experts (YY, HL) who had 24 and 16 years of experience in US thyroid examination. Consensus were made when disagreements appeared and each of them reviewed half of the nodules without knowing the pathological state of ETE. Based on previous studies, sonographic ETE was suspected on US images when a capsular abutment or protrusion presented.¹³ Capsular abutment referred to nodules' contact of the thyroid capsule, which could be graded by their perimeter ratio (< 25%, 25–50% or > 50%). While protrusion referred to nodules protruding thyroid capsule and invading the muscles or soft tissue around the thyroid gland, including the trachea, fat space between the trachea, esophagus, esophagus sulcus, cervical sheath vessels, et al.^10,13 The predictive performance of sonographic method was calculated using all nodules and compared by the AUC of the receiver operating characteristic (ROC), sensitivity, specificity, accuracy, positive predictive value (PPV) and negative predictive value (NPV).

Different US systems were employed in nodules' examination, including Esaote MyLab Class C (Esaote, Italy), Siemens ACUSON OXANA 2 (Siemens Medical Solutions, USA), GE Volume E8 (GE Medical Systems, USA), Hitachi HI VISION Preirus (Hitachi-Aloka Medical, Japan) and Mindray Resona 7T (Mindray Medical International, China). High-frequency linear probes (5 MHz to 14 MHz) were used depending on US systems. In order to have the best image for nodule, the depth and width of US systems were adjusted individually. Images of transverse and longitudinal section of nodules were saved as JPG files before the largest cross section of nodules were chosen for feature extraction. The US experts who were blinded to the pathological state of ETE reviewed the images, measured the US findings of nodules and delineated the region of interest (ROI) manually by Paint Win10. Consensus were made when disagreements appeared. Each expert review half of the nodules, respectively.

From those ROIs of US images, a total of 1769 US radiomic features were extracted using PyRadiomics (version 3.7, http://pyradiomics.readthedocs.io/en/latest/index.html). The features included 9 shape features, 360 first-order statistical features and 1400 textural features. Gray-level size-zone matrix (GLSZM), gray-level co-occurrence matrix (GLCM), gray-level dependence matrix (GLDM) and gray-level runlength matrix (GLRLM) consisted of the gray matrices in textural features. What's more, high-dimensional features were acquired by filters including Laplacian of Gaussian (LOG) with different sigma values (1.0 mm–10.0 mm with step 1.0 mm), square root, square, exponential, wavelet with 2D transform (low-pass/high-pass, LH; low-pass/low-pass, LL; high-pass/high-pass, HH; high-pass/low-pass, HL;), gradient and logarithm.

All nodules were randomly assigned to training group (237 nodules) and validation group (100 nodules) by a ratio of 7:3. In order to show the randomness of grouping, the clinical characteristics between training group and validation group were compared. Then all features were enrolled in the following progression, including 1769 US radiomic features and 34 clinical characteristics. To eliminate redundant and irrelevant features as well as reduce variables recruited into model, a two-step features selection method was applied after standardization. Firstly, the minimum redundancy maximum relevance (mRMR) was used to select the most relevant features. Secondly, they were further selected by least absolute shrinkage and selection operator (LASSO) algorithm and candidate features were obtained. Ten-fold cross validation was used to avoid overfitting.²⁰

In total, nine predicting models were constructed using these candidate features, including k nearest neighbors (KNN), binary logistics regression (LR), support vector machine (SVM), naivebayes (NB), randomforest (RF), decision tree (DT), adaptive boosting (adaboost), the extreme gradient boosting (xgboost) and gradient boosting machine (GBM). Their predictive performances were compared by AUC, sensitivity, specificity, accuracy, PPV and NPV. And the model with the highest AUC and accuracy was chosen as US radiomics method. Finally, the predictive performances of US radiomics method and sonographic method were compared. The statistical process was shown in Figure 1.

Continuous variables were presented as mean ± standard deviation (SD) and were measured by Mann-Whitney U test or Student's t-test depending on the results of normality analysis. Categorical variables were presented as number (frequency) and Chi-square analysis or Fisher exact test were performed afterwards. SPSS software (version 19.0, IBM) was employed for the calculation mentioned above. The random allocation, standardization, mRMR, LASSO algorithm, KNN, binary LR, SVM, NB, RF, DT, adaboost, xgboost, GBM, DeLong's test and other statistical analysis were carried out with R software (version 4.0.3, MathSoft, http://www.r-project.org). DeLong's test was employed for statistical significance in AUC comparison.²¹ A value of p < 0.05 was considered statistically significant.

Clinical characteristics of training group and validation group were summarized in Table 1, including age, sex, biochemical results and US findings. The training group comprised of 237 nodules with a positive ETE rate of 30.38%. And the validation group comprised of 100 nodules with similar positive ETE rate of 31.00%. No significant differences were found in any characteristics between the two groups. Of all the nodules, images acquired by Esaote, Siemens, GE, Hitachi and Mindray systems accounted for 70.33% (237/337), 15.13% (51/337), 8.90% (30/337), 4.15% (14/337) and 1.48% (5/337).

The predictive performance for sonographic method was summarized in Table 2. And its ROC was shown in Figure 2A. A relatively high accuracy (0.70), specificity (0.7279) and NPV (0.9145) was found for this method. However, its sensitivity (0.5349), PPV (0.2233) and AUC (0.569) was at a low level.

Figure 2.

The ROC curves of sonographic method and US radiomics method for ETE predicting. (A) The ROC curve of sonographic method in all nodules. (B) The ROC curve of US radiomics method in validation group.

ETE = extrathyroidal extension; ROC = the receiver operating characteristic; US = ultrasonic

After standardization, the top-99 most relevant features to ETE were retained using mRMR. And 81 candidate features were selected using LASSO algorithm afterwards. Among the nine predicting models, xgboost model had the highest accuracy (0.77) and AUC (0.813) with p-value <0.05 compared to other models (Table 2). Meanwhile, with its relatively high specificity (0.8116), and NPV (0.8485), it was chosen as US radiomics method. Then the predictive performance of US radiomics method and sonographic method was compared. And the results showed that the AUC of US radiomics method was significantly higher than that of sonographic method with p-value <0.05. Meanwhile, the accuracy, sensitivity, specificity and PPV of US radiomics method was higher than that of sonographic method (Figure 2B, Table 2). Thus, we believe US radiomics method surpassed sonographic method in ETE predicting (Figure 1).

The image of feature importance for xgboost model was shown in Figure 3. And the top-6 features were chosen at a cut-off of 0.05 (Gain value). The partial dependence profiles of the them were revealed in Figure 4, in which the variation trends between features and ETE states were presented.

Figure 3.

The image of feature importance for xgboost model. (A-D) the top 1–20, top 21–40, top 41–60 and top 61–80 feature importance of xgboost model.

xgboost=the extreme gradient boosting; GLSZM=Gray-level size-zone matrix; GLCM=gray-level co-occurrence matrix; GLDM=gray-level dependence matrix; GLRLM=gray-level runlength matrix; LH=low-pass/high-pass; LL=low-pass/low-pass; HH=high-pass/high-pass; HL=high-pass/low-pass; ANTITPOAB=anti-thyroid peroxidase antibody; AST=asparate aminotransferase; RBC=red blood cell; FT3=free triiodothyronine 3; Ca=calcium ion; TT4=total triiodothyronine 4; TG=thyroglobulin; TSH=thyroid stimulating hormone; NEUT=neutrophil; PLT=platelets.

Figure 4.

Partial dependence profile of the top-6 important features in xgboost model. (A-F) Partial dependence profile of logsigma6.0mm3D-GLCM-MCC, logsigma5.0mm3D-GLCM-Imc1, square-firstorder-RootMeanSquared, waveletLL-GLSZM-SmallAreaLowGrayLevelEmphasis, logarithm-GLSZM-SmallAreaHighGrayLevelEmphasis and logsigma6.0mm3D-firstorder-90Percentile in xgboost model in training group. xgboost, the extreme gradient boosting; GLSZM, Gray-level size-zone matrix; GLCM, gray-level co-occurrence matrix; LL, low-pass/low-pass.