This study was conducted with patients retrospectively enrolled from a public database and our hospital. The public database was originally submitted to The Cancer Genome Atlas (TCGA) and the Cancer Imaging Archive (TCIA) by the contributing institutions under an Institutional Review Board-approved protocol. The second part of the study was approved by the ethics committee of our hospital, and the requirement for individual consent for this retrospective analysis was waived between January 2021 and December 2021. Clinical, breast MRI, and transcriptome data for all patients from TCGA and TCIA, and our hospital were integrated and analyzed. A general overview of the analysis protocol is shown in Figure 1 .

The inclusion criteria were as follows: (i) patients had histologically confirmed unilateral primary breast cancer; (ii) patients who underwent breast surgery without neoadjuvant chemotherapy, and sentinel lymph node biopsy or ALN dissection with curative intent; (iii) patients with suspected positive ALN by clinical and/or imaging examinations and plan to receive neoadjuvant therapy, lymph node status was determined by ultrasound-guided core needle biopsy before neoadjuvant therapy (20, 21); (iv) availability of clinical data (age), T staging, complete pathological data for ALN and molecular subtype; (v) breast DCE-MRI was conducted before anti-tumor treatment and core needle biopsy for evaluating tumor extension, intraductal spread, and the presence of multicentric or multifocal lesions; and (vi) availability of transcriptome data. The exclusion criteria were as follows: (i) insufficient MRI quality to obtain measurements and (ii) patients with multifocal lesions. (iii) patients with negative ALN by clinical and/or imaging examinations and plan to receive neoadjuvant therapy. Details are provided in Figure 2 .

Fifteen patients with breast cancer from our hospital were used for the external validation group. Postoperative tumor samples were extracted from the primary tumor sites, and snap-frozen breast cancer samples were collected by the Tumor Samples and Genetic Information Resource Library. Next, the sample from the hospital was processed using RNA extraction and sequencing. Total RNA was extracted from tissue samples with TRIzol reagent (Invitrogen, Carlsbad, CA, USA). RNA purity was monitored on 1% agarose gel. Then, cDNA was synthesized from total RNA using the PrimeScript RT reagent Kit with gDNA Eraser (TaKaRa, RR047A). Then, all patients’ transcriptomic data were analyzed for differentially expressed genes (DEGs). The DEGs between the groups with and without ALNM were screened using DESeq (http://bioconductor.org/packages/release/bioc/html/DESeq.html) to detect DEGs with gene symbol annotation with thresholds of P < 0.05 and |fold change| > 2.

The public MRI data of breast cancer were downloaded from TCIA. The breast cancer MRI was performed at our hospital using a commercially available 1.5-T system (Philips Ingenia, the Netherlands) with an eight-channel breast array coil. Routine precontrast MRI included T1-weighted imaging, fat-suppressed T2-weighted imaging, and diffusion-weighted imaging. The gadolinium contrast agent (Hengrui, Jiangsu, China) was injected at a rate of 2 mL/s for a dose of 0.1 mmol/kg, followed by a 20 mL saline flush at a rate of 2 mL/s. One precontrast phase image and seven postcontrast phase images were acquired using the axial e-THRIVE polyphase sequence with the following parameters: TR/TE = 5.1/2.5 ms; matrix = 252 × 348; flip angle = 10°; pixel size = 1 × 1 mm; and slice thickness = 2 mm.

The regions of interest (ROIs) were delineated manually in each slice of the T1+C data (the peak-enhanced phase of the multiphase contrast-enhanced MRI selected according to the time-intensity curve) by excluding the air, necrosis, and calcification areas via the Dr. Wise Multimodal Research Platform (https://keyan.deepwise.com,V1.6.3) (Beijing Deepwise & League of PHD Technology Co., Ltd., Beijing, China). Each DCE-MRI case was reviewed by two experienced radiologists (XXW and LL with 9 and 8 years of experience, respectively). After manual segmentation, consistency intraclass correlation coefficient (ICC) (22) analysis was performed to measure the features’ observer repeatability. ICC scores greater than 0.8 are generally considered to indicate good repeatability. We standardized the image processing procedure according to the image biomarker standardization initiative reference manual (23). Next, 1651 radiomics features were extracted from the normalized image using Pyradiomics (http://www.radiomics.io/pyradiomics.html). Then, they were normalized with Z scores to obtain a standard normal distribution of image intensities. The radiomics features were composed of the following eight types of features: first-order statistics, shape-based (3D), shape-based (2D), gray level cooccurrence matrix (GLCM), gray level run length matrix (GLRLM), gray level size zone matrix (GLSZM), neighboring gray-tone difference matrix (NGTDM), gray level dependence matrix (GLDM). All of these features have generally been used in previous radiomics studies (24, 25).

We used the Boruta method (R3.6.1 with Boruta version 5.2.0) from a coarse to fine feature reduction strategy to reduce both genomics and radiomics features (24). It was divided into the two steps described below. (i) Univariate analysis was performed using Student’s t test or the Mann–Whitney U test to compare genomics and radiomics features between the groups with and without LNM. All features were ranked in ascending order according to the P value, and the top 5% of features were retained for further analysis. (ii) Spearman correlation analysis was used to eliminate redundant features. All genomics and radiomics features with correlation coefficients > 0.85 were detected, and the features with lower P values were retained. Then, the selected features were used for further model construction and validation.

Patients from the TCIA-TCGA were randomly divided into training and testing group at a ratio of 3:1, and 15 patients from our hospital were enrolled as the external validation group. Logistic regression analysis (Python3.7 with sklearn version 1.1.2) was used to construct radiomics, genomics and radiogenomics models for predicting ALN status. The radiogenomics model was conducted through dimensionality reduction by integrating the features of the radiomics model and the features of the genomics model. All models were trained based on TCGA datasets, and fivefold cross validation was used to determine the parameters of the logistic regression models. We evaluated the prediction performance of the radiomics, genomics and radiogenomics models in the training, testing and external validation groups by constructing receiver operating characteristic (ROC) curves and calculating areas under the curves (AUCs). Decision curve analysis was applied to estimate the clinical utility of the three models. In addition, a calibration curve was generated to evaluate the consistency between the predicted value and the true value.

Statistical analyses were performed by commercially available statistical software (SPSS software, version 25.0; Armonk, US). Descriptive statistics were summarized as the means ± standard deviations. Categorical variables were expressed as numbers. Continuous clinical variables were compared using Student’s t test or the Mann–Whitney U test. For categorical variables, chi-square tests were used to test differences between groups. The performances of the radiomics model, genomics model and radiogenomics model were compared using the Delong test. For all tests, P < 0.05 was considered as statistically significant.

The baseline characteristics in the training, testing, and external validation group are presented in Table 1 . A total of 126 patients with breast cancer from one hospital and public database were enrolled in the study. The training group (53.39 ± 11.44 years) included 84 patients, and the testing group (55.96 ± 10.56 years) included 27 patients from the TCGA and TCIA datasets. The external validation group (50.53 ± 7.24 years) included 15 patients who were from our hospital. There were no significant differences in age or LNM among the three groups (P = 0.91 and 0.92).

For differential gene expression analysis, 136 DEGs were found between patients with and without ALNM ( Figure S1 ). As shown in Figure S1A , expression levels of 94 genes were increased and 42 genes were decreased between patients with and without ALNM. Heatmap representation of DEGs revealed a similarity between patients with and without ALNM ( Figure S1B ). Then, 9 genomics features were selected through Boruta analysis ( Figure 3A ). The 9 genomics features were used to develop the genomics model for predicting ALNM in breast cancer. The AUCs of the training, testing, and external validation group were 0.80, 0.67 and 0.52, respectively ( Figure 4A ). The AUC and corresponding sensitivity, specificity, positive predictive value (PPV), and negative predictive value (NPV) values of the genomics model are displayed in detail in Table 2 . The calibration curves for the probability of the genomics model for ALNM in the training group, testing group and external validation group are shown in Figures S2A–C .

A total of 1651 features were extracted from DCE-MRI from every patient. After the reduction by Boruta ( Figure 3B ), 3 radiomics features were used to construct and validate the radiomics model. The AUCs of the prediction performance of ALNM in breast cancer of the radiomics in the training and testing group were 0.72 and 0.68, respectively ( Figure 4B ). In addition, the performance of the radiomics model was validated in an independent external validation group, and the AUC of the validation group was 0.71. The AUC and corresponding sensitivity, specificity, PPV and NPV values of the radiomics model are displayed in detail in Table 2 . The calibration curves for the probability of the radiomics model for ALNM in the training group, testing group and external validation group are shown in Figures S2D–F.

After stepwise logistic regression with both direction, 5 genomics features and 3 radiomics features were finally selected to develop a radiogenomics model. We performed correlation analysis between the 3 radiomics and 5 genomics features, as shown in Table S1 . The results showed that the correlations between the genomics and radiomics features were slightly weaker. The AUCs of predicting ALNM for the radiogenomics model in the training group, testing group and external validation group were 0.84, 0.75 and 0.82, respectively ( Figure 4C ). The AUC and corresponding sensitivity, specificity, accuracy, PPV, and NPV values of the three groups are detailed in Table 2 . Good agreement of the radiogenomics model between the observation and prediction was assessed by the calibration curve, which showed that the bias-corrected line lay close to the ideal curve in the training, testing, and validation group ( Figures S2G–I ). In addition, the decision curve showed that the radiogenomics model could add more benefit to the prediction of ALNM than the genomics model and radiomics model in the three groups ( Figure 5 ).

The prediction profiles of the radiogenomics model compared with the genomics model and radiomics model by the DeLong test in the training group, testing group and external validation group are shown in Figure 6 . The results showed that the radiogenomics model significantly improved the AUC compared with the radiomics model in the training group (P = 0.01) and showed superior performance compared with the genomics model in the external validation group (P = 0.02) ( Table 3 ).