Between January 2014 and December 2014, 793 patients underwent surgery for breast cancer at our institution. The exclusion criteria were as follows: (1) 175 patients with ductal carcinoma in situ, (2) 51 patients with masses larger than 4 cm at preoperative US that could not be fully included in a single standard US image, (3) 21 patients who underwent re-operation or surgery for recurrence, (4) 20 patients referred to our institution after excisional biopsy at an outside clinic, (5) 17 patients with non-mass lesions with uncertain boundaries due to vague regions of interest, and (6) 13 patients who underwent neoadjuvant chemotherapy without initial histological confirmation of the axillary lymph nodes. After exclusion, 496 patients (mean age, 52.5±10.9 years) with breast cancer were included in our study. Among them, 306 patients who underwent surgery between January 2014 and August 2014 were enrolled as the training cohort, and 190 patients who underwent surgery between September 2014 and December 2014 were enrolled as the validation cohort (Fig. 1).

US examinations were performed by 10 radiologists using two different ultrasound machines (iU22, Phillips Medical Systems, Bothell, WA, USA; LOGIQ E9, GE Healthcare, Milwaukee, WI, USA) with linear array transducers. If a patient underwent multiple US examinations prior to surgery, we selected the US examination taken at the time at which a suspicious mass was detected at our institution. The median interval between the initial US examination and surgery was 15 days (range, 2 to 335 days). A radiologist (E.K.K.) retrospectively reviewed the US images and collected data regarding mass size, tumor location, multiplicity in a single breast, and skin-to-tumor distance.

We reviewed post-biopsy pathologic reports to investigate cancer type (ductal, lobular, or other) and estrogen receptor, progesterone receptor, human epidermal growth factor receptor 2 (HER2), and Ki67 status. Mixed ductal and lobular cancer, mucinous cancer, invasive micropapillary carcinoma, tubular carcinoma, and other mixed types were classified as "other" with regard to cancer type. Estrogen receptor and progesterone receptor positivity were defined as immunoreactivity of 1% or higher for tumor cell nuclei, and Ki67 positivity was defined as immunoreactivity of 14% or higher. In cases of equivocal HER2 overexpression, an amplification ratio of 2 or higher on fluorescence in situ hybridization testing was considered to indicate HER2 positivity. We also collected data from electronic medical records regarding whether each patient underwent neoadjuvant chemotherapy, and we obtained information regarding lymphovascular invasion and histologic grade from postoperative pathologic reports; these data were used for baseline comparison of clinicopathologic features between the training and validation cohorts. In the development of the clinicopathologic model, these variables were excluded because we aimed to develop the model in a preoperative setting. The standard reference for lymph node status was based on the results of sentinel lymph node biopsy or axillary lymph node dissection. In patients who received neoadjuvant chemotherapy, we instead referred to the results of fine-needle aspiration before starting treatment.

A radiologist with 1 year of experience in breast imaging (S.E.L.) selected one axial image among the US images of each breast mass and cropped the image to remove the space used for informative text. After the image was resampled to a pixel size of 0.2 mm, a region of interest along the mass margin was imaged semi-automatically using MIPAV software version 8.0.2 (NIH, Bethesda, MD, USA; open-source, https://mipav.cit.nih.gov) and converted into mask files for feature extraction by the same radiologist. Another resident radiologist (Y.S.), a third-year resident, independently performed segmentation of 50 randomly-chosen masses to evaluate interobserver reproducibility.

A radiologist with 4 years of experience in data science (S.K.) extracted features from the mask files using Pyradiomics software (version 2.0.0, open-source, https://pyradiomics.readthedocs.io/en/latest). A total of 444 radiomic features were extracted from the original and derived (wavelet-transformed) images. For each radiomic feature of the 50 randomly-selected masses, the intraclass correlation coefficient was calculated between two radiologists, and 39 features with coefficients of less than 0.75 were excluded. For the remaining 405 radiomic features, z-score normalization was applied to standardize the values. Features with Spearman correlation coefficients greater than 0.95 were represented by a single feature that showed the widest range among clustered features through the hierarchical clustering process. In total, 125 features were selected to optimize reproducibility and redundancy. They consisted of 40 features from the original images (4 shape features, 10 first-order features, 15 gray-level co-occurrence matrix features, 4 gray-level run-length matrix features, 5 gray-level size-zone matrix features, and 3 gray-level dependence matrix features) and 85 features from the wavelet-filtered images. Image processing was performed using ITK Python packages (version 4.13.2, open-source, https://itk.org/ITK/resources/software.html). Fig. 2 shows the process from US image acquisition to radiomics model development.

Finally, we selected radiomic features using penalized logistic regression under the least absolute shrinkage and selection operator (LASSO) model with 5-fold cross-validation in the training cohort. A rad-score was computed via a linear combination of the selected features weighted by each coefficient. The area under the curve (AUC) was calculated in the training cohort using the selected features with a 95% confidence interval (CI). A preoperative clinicopathologic model was established using multivariable logistic regression with the variables that had P-values less than 0.05. We calculated the predictive performance levels of the clinicopathologic and combined clinicopathologic-radiomics models to evaluate the incremental value of the radiomic model via the Delong test for two receiver operating characteristic curves. Performance was independently evaluated in the validation cohort. Differences in clinicopathologic characteristics between the training and validation cohorts were assessed using the Mann-Whitney U test and the chi-square test.

Statistical analyses were performed using R software version 3.6.1 (R Foundation for Statistical Computing, Vienna, Austria; http://www.R-project.org). P-values of less than 0.05 were considered to indicate statistical significance.

The clinicopathologic characteristics of the patients in the training and validation cohorts are summarized in Table 1. Lymph node positivity was present in 30.1% (92 of 306) of the patients in the training cohort and in 32.1% (61 of 190) of those in the validation cohort, which did not constitute a statistically significant difference (P=0.689). Similarly, no significant difference was observed between the training and validation cohorts for any other factor.

Based on multivariable logistic regression, the independent preoperative clinicopathologic factors identified as predictors of axillary lymph node metastasis were mass size on US, tumor location (outer, medial, or subareolar), tumor type (ductal, lobular, or other), and multiplicity (Table 2). Age, skin-to-tumor distance, distance from the nipple, and receptor status showed no significant association with lymph node metastasis. The predictive performance of the clinicopathologic model was moderate, with AUC values of 0.760 (95% CI, 0.703 to 0.817) in the training cohort and 0.708 (95% CI, 0.630 to 0.786) in the validation cohort.

Of the 125 features that were originally chosen, 23 were selected in the training cohort using the LASSO logistic regression model (Table 3, Fig. 3). Among the 23 radiomics features, 'first order_kurtosis' was a dominant feature in our radiomics model, as it was associated with the highest value of the coefficient. A clinicopathologic model was developed with four factors: tumor size, location, subtype, and multiplicity. The predictive performance of the radiomics model was comparable to that of the clinicopathologic model, with AUCs of 0.812 (95% CI, 0.760 to 0.864) in the training cohort and 0.831 (95% CI, 0.773 to 0.889) in the validation cohort. The radiomics model showed significantly better predictive performance than the clinicopathologic model in the validation cohort (P=0.013) (Table 4).

To evaluate the incremental value of the radiomics model, we developed a combined model using the radiomics score and the four aforementioned clinicopathologic factors. The AUC of the combined model was 0.858 (95% CI, 0.814 to 0.902) in the training cohort, which was significantly better than the performance of the clinicopathologic model alone (AUC, 0.760; P=0.007).

When we applied these models to the validation cohort, the AUC of the combined model was 0.810 (95% CI, 0.745 to 0.876). The combined model performed significantly better than the clinicopathologic model in the prediction of axillary lymph node metastasis (AUC, 0.708; P=0.048) (Table 4, Fig. 4).