This retrospective case-control study was conducted at two institutions. Its aim was to develop and validate a model for the individual prediction of CLNM in patients with PTMC. Data were collected from 327 patients at Cangzhou Hospital of Integrated TCM-WM·Hebei for the training cohort between January 2018 and December 2020. Data were obtained from Hebei Medical University Health System (n = 153) during the same period for the validation cohort. The study was approved by each institution’s institutional review board. Written informed consent was obtained from all patients. This study was conducted in accordance with the Declaration of Helsinki (2013 revision) and followed the Strengthening the Reporting of Observational Studies in Epidemiology reporting guideline.

The inclusion criteria were as follows: (1) age > 18 years; (2) histologically confirmed PTMC that was clinical N0 assessed by preoperative US; and (3) having undergone standard thyroidectomy and PCLND.

The exclusion criteria were as follows: (1) concurrent malignant disease of other organs; (2) recurrent or metastatic thyroid cancer; (3) history of previous surgery or radiotherapy of the neck; and (4) incomplete US or clinicopathological data (Figure S1).

All clinicopathological data, including age, sex, and BRAF status, were collected from the medical records. For the BRAF mutation analysis, an AmoyDx® BRAF Mutation Detection Kit (V2) (ADx-BR02; Amoy Diagnostics Co., Ltd., Xiamen, China) was utilized. The detection of the mutation was performed using a next-generation sequencing method, followed by a real-time fluorescence polymerase chain reaction–amplification refractory mutation system.

US examinations were performed using a 5–14-MHz transducer (Siemens, ACUSON Sequoia, Siemens Medical Solutions USA, Inc., Malvern, PA, USA) by radiologists with at least 8 years of experience performing thyroid US evaluations. The US characteristics of each lesion, including diameter, texture, echo, boundary rule, presence of calcification, and capsule invasion, were evaluated by an independent radiologist (W.L.).

US images were retrieved from the picture archiving and communication system (Carestream, Toronto, Canada) for further feature extraction. The region of interest (ROI) of each lesion was manually segmented on the largest diameter image using ITK-SNAP software. To measure the interobserver agreement, all the manual segmentations were performed by two experienced radiologists who were blinded to patients’ characteristics. Moreover, one of the radiologists delineated the ROIs again after two weeks to measure the intraobserver agreement. The ROIs delineated by this radiologist in the second round were used for subsequent feature extraction. The radiomics feature extraction was performed using the open-source platform Pyradiomics (version 3.1.0). This platform allows the extraction of 851 radiomics features, which can be classified into shape features, first-order features, gray-level co-occurrence matrix features, gray-level size zone matrix features, gray-level run length matrix features, and gray-level dependence matrix features. The interclass correlation coefficient (ICC) was used to evaluate the inter- and intraobserver agreements of the feature extraction. Features with good consistency (ICC > 0.75) were subjected to further analysis.

Before the feature selection, the values of the extracted features were standardized with z scores. A three-step procedure was performed to select the robust radiomics features in the training cohort [14]. First, a univariable logistic regression analysis was performed to identify significant CLNM predictors with P < 0.05. Second, the Pearson correlation coefficient for each of the two features was calculated, and we excluded the one with a higher P value for those feature pairs with a strong correlation (Pearson r > 0.90). Overall, 107 features were screened out for the last selection. Third, the least absolute shrinkage and selection operator (LASSO) logistic regression model was performed to determine the optimal combination of radiomics features and calculate a radiomics (Rad) score by 10-fold cross-validations via the 1–standard error criteria.

Categorical variables are expressed as frequencies and percentages and were compared using the chi-squared test. We calculated the hazard ratios (HRs) and 95% confidence intervals (CIs) of CLNM using the logistic regression model with uni- and multivariate analyses. Pearson correlation coefficients were calculated to evaluate correlations among the parameters.

To provide a quantitative tool to predict the individual probability of CLNM, we generated the radiomics nomogram based on the multivariate analysis of the training cohort. The discrimination of the nomogram was assessed using receiver operating characteristic curves by calculating the area under the curve (AUC). The model’s calibration was assessed using calibration curves by comparing the predicted and actual probability. A decision curve analysis was utilized to assess the clinical usefulness of the nomogram.

All statistical analyses were performed using SPSS software (version 22.0; IBM Corporation, Armonk, NY, USA) and R software version 4.1.3 (R Foundation for Statistical Computing, Vienna, Austria). Statistical significance was set at a two-tailed value of P < 0.05.