Machine learning magnetic resonance imaging radiomics predicts axillary lymph node metastasis in invasive breast cancer

In current clinical practice, the standard evaluation for axillary lymph node (ALN) status in breast cancer is based on the invasive procedure and many patients will suffer from operative associated complications. Hence, a novel signature incorporated tumor and lymph node magnetic resonance imaging (MRI) radiomics, clinical and pathological characteristics, and molecular subtypes based on the machine learning approach was established to accurately identify ALN metastasis in early-stage invasive breast cancer patients. Although the misjudgment of ALN status by clinicians according to preoperative MRI are common during clinical practice and even the senior radiologists make mistakes sometimes, this multiomic radiomic signature showed the superiority over clinicians and could precisely discriminate ALN metastasis among different molecular subtype patients. Furthermore, the association between MRI radiomic features and tumor-microenvironment features including immune cells, long non-coding RNAs, and types of methylated sites were found, which revealed the potential biological underpinning of MRI radiomics.


Introduction
The importance of axillary lymph node (ALN) status in early-stage invasive breast cancer was highlighted in 10-year follow-up of the ACOSOG Z0011 and IBCSG 23-01 randomized clinical trials. [1][2] Since sentinel lymph node biopsy (SLNB) was included in the breast cancer standard diagnosis procedure from 2005, it have greatly improved the accuracy of ALN status judgment and changed the traditional concept that women could be away from axillary lymph nodes dissection (ALND) if they had fewer than three positive sentinel lymph nodes. [3][4] However, SLNB is still an invasive procedure that may bring complications such as lymphedema and upper limb numbness and its accuracy is also limited by the radiopharmaceuticals and modalities of radiocolloid injection. [5][6] As a consequence, it's urgent and indispensable to develop a reliable non-invasive evaluation approach of ALN status in patients who had early-stage invasive breast cancer.
Nowadays, the development of radiomics in the aspects of tumor diagnosis, treatment decision and prognosis prediction is encouraging and breast cancer is one of the pioneers in the exploration. [7][8][9] Magnetic resonance imaging (MRI) radiomics-based approaches for predicting molecular subtype, recurrence risk and survival, pathologic complete responses of breast cancer patients showed powerful predictive ability and the prospect of extended application. [10][11][12][13] Moreover, the association between MRI radiomics with breast tumor microenvironment have been corroborated in previous studies. [14][15] As MRI examination has been recommended for preoperative estimation of breast cancer patients, in this retrospective multicenter study, a new preoperative non-invasive approach based on MRI radiomics was established to predict ALN metastasis in early-stage invasive breast cancer patients, which reduces unnecessary SLNB or ALND and would decrease the incidence of associated postoperative complications to improve the life quality of patients.

Results
Patients' characteristics A total of 1,161 early-stage invasive breast cancer patients were retrospective recruited, of which 1,088 patients were eligible for this study. The study design was shown in Supplementary Fig. 1 and the study work ow was shown in Fig. 1. The Table 1 showed the clinicopathological characteristics of patients in the training cohort (n=803), the prospective-retrospective validation cohort (n=106), and the external validation cohort (n=179). 389 (35.75%) of 1,088 patients were initially diagnosed as positive ALN status by radiologists though MRI, but 106 (27.25%) of 389 patients didn't have ALN metastasis (ALNM) and were con rmed as negative ALN status with pathologically examination. Rather, 190 (28.27%) of 672 patients with clinical negative ALN status were found to have pathological positive ALN status.
Distinguishing ALN metastasis by tumor and lymph node radiomic signature The key radiomic features of ALN and tumor regions were selected from a total of 5,178 quantitative features by the random forest algorithm. Overall distribution of key radiomic features from contrastenhanced T1-weighted imaging (T1+C), T2-weighted imaging (T2WI), and diffusion-weighted imaging quantitatively measured apparent diffusion coe cients (DWI-ADC) sequences among patients with and without ALNM in the training cohort was demonstrated in Fig. 2A. Incorporating these key features of ALN region to predict ALNM yielded AUC values of 0.85, 0.61, and 0.81 in the training cohort, the prospective-retrospective validation cohort, and the external validation cohort, respectively ( Supplementary Fig. 2). Simultaneously, incorporating three-sequence key features of tumor region for ALNM prediction achieved an AUC of 0.78, 0.59, and 0.63 in the training cohort, the prospectiveretrospective validation cohort, and the external validation cohort, respectively ( Supplementary Fig. 3).
The AUC values for ALNM prediction of combining multi-sequence features were higher than incorporating single-sequence features in the training and the validation cohorts (Supplementary Table  1).
When combined both ALN and tumor regions features, the ALN-tumor radiomic signature was constructed and illustrated good performance in detecting ALNM in the training cohort (AUC, 0.88), the prospective-retrospective validation cohort (AUC, 0.87), and the external validation cohort (AUC, 0.87), which outperformed the ALN or tumor radiomic signature alone (Fig. 2B). The detailed evaluation indicators for model performance including sensitivity and speci city were summary in Table 2.
Distinguishing ALN metastasis by multiomic radiomic signature In the univariate analysis, which was presented in Supplementary Table 2, ve differentially expressed clinical characteristics were found to be sassociated with ALN status in the training cohort, including age (P = 0.012), clinical T stage (P < 0.001), clinical N stage (P < 0.001), Ki67 expression (P = 0.010), and molecular subtype (P = 0.033). To develop a more precise and clinically applicable method that could predict an individual's ALN status, the multiomic radiomic signature incorporated all key radiomic features of ALN and tumor region with clinical characteristics, pathological characteristics and molecular subtype that signi cantly associated with ALNM was built and showed better performance of ALNM prediction, which achieved the higher AUC (0.90) in the training cohort, the prospective-retrospective validation cohort (AUC, 0.93), and the external validation cohort (AUC, 0.91), respectively (Fig. 2C). The sensitivity and speci city of the multiomic radiomic signature were list in Table 2. As for combination of these clinicopathological characteristics and molecular subtype, it achieved an AUC of 0.74, 0.68, and 0.72 for predicting ALNM in the training cohort, the prospective-retrospective validation cohort, and the external validation cohort, respectively ( Supplementary Fig. 4).
The multiomic radiomic signature also presented the ability of discriminating ALNM patients with 1, 2,  Table 3).
Besides, to further assess the added value of the multiomic radiomic signature to the ALN status, we conducted subgroup analysis within patients with different molecular subtypes. Encouragingly, the multiomic radiomic signature could identify ALNM patients in the subgroups of Luminal A (AUC, 0.91, 0.91, respectively), Luminal B (AUC, 0.89, 0.92, respectively), human epidermal growth factor receptor 2 (Her-2) positive (AUC, 0.90, 0.76, respectively), and triple-negative breast cancer (TNBC) (AUC, 1.00, 1.00, respectively) patients in the training and external validation cohorts (Supplementary Table 3). In addition, when strati ed by other factors such as age, Ki67 expression level, and clinical T stage, the AUC value of the multiomic radiomic signature remained at 0.86-1.00 in these subgroups (Supplementary Table 3).
Decision curve analysis (DCA) was conducted to determine the clinical usefulness of the multiomic radiomic signature by quantifying the net bene ts at different threshold probabilities. The decision curve showed that if the threshold probability is > 5%, using the multiomic radiomic signature to predict ALNM adds more bene t than either the ALN-tumor radiomic signature or the radiologists' diagnosis (Fig. 2D).
According to the radiomic score of the multiomic radiomic signature, an optimal cutoff value (0.334) was generated to classify patients into high-and low-score groups in the training cohort. High-score patients with high risk of ALNM had signi cantly shorter DFS compared with the low-score group (HR 0.43, 95% CI 0.21-0.86, P=0.014; Supplementary Fig. 5).
During clinical practice, patients' clinical ALN status judged by the preoperative MRI are commonly inconsistent with the pathological one and even the senior radiologists make mistakes sometimes. The multiomic radiomic signature could precisely recognize ALNM among patients with different clinical tumor stage in the training and the validation cohorts (Supplementary Table 3). As shown in Fig. 3, patient 1 had pathological positive ALN status but was considered as a non-ALNM patient by radiologists though MRI before surgery. On the contrary, patient 2 was initially diagnosed as ALNM by radiologists but found to be a non-ALNM patient by pathologically examination. The ALN status of two patients could be accurately assessed through the multiomic radiomic signature by the cutoff value of 0.334.

Radiomics associatied with tumor microenvironment
According to the cutoff values of radiomic score from T1+C and T2WI sequences signatures of tumor region in the training cohort, 91 breast cancer patients from The Cancer Genome Altas (TCGA) and The Cancer Imaging Archive (TCIA) with T1+C and T2WI sequences MRI were classi ed into two group. In total, 1381 T1+C and T2WI sequences-based differentially expressed genes (DEGs) were obtained among low-score and high-score patients. Next, the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were performed to further examine the biological functions of the identi ed radiomic-based genes. The GO enrichment analysis indicated that the radiomic-based genes were enriched in various physiological metabolic processes, such as affection of transmembrane transporter activity, NADP binding and ATPase complex (Supplementary Table 4 and Supplementary Fig. 6). The KEGG pathway enrichment analysis found these genes were involved in the oxidative phosphorylation signaling pathway (Supplementary Table 4).
We further explored the association between MRI radiomic features and tumor microenvironment

Discussion
Unlike previous studies just focusing on the tumor region for predicting ALNM, [16][17][18] in this study, ALNtumor radiomic signature that combined multi-sequence key radiomic features of ALN and tumor regions showed high predictive ability and could be applied to predict the ALN status precisely and noninvasively. The multiomic radiomic signature integrated radiomic features above and ALNM-associated clinicopathologic characteristics, molecular subtype demonstrated better performance for predicting ALNM. In addition, the multiomic radiomic signature contributed to improve the accuracy of judgement for ALN status in patients who had inconsistent clinical ALN status with pathological ALN status, and it could accurately recognize ALNM in patients with different number of positive lymph nodes, different clinical T stage, different age, different Ki67 expression level, and even different molecule subtypes. Besides, the key radiomic features were also found to be associated with tuomr microenvironment. Therefore, the multiomic radiomic signature was an effective tool for ALNM prediction that provide useful message for diagnosis and treatment decision preoperatively and meanwhile could preventing unnecessary ALND and even SLNB in early-stage invasive breast cancer patients.
Breast cancer patients have distinct ALN status when they are newly diagnosed, and the treatment plans are so as different depend on their ALN status. It has been demonstrated that patients with positive ALN status had poorer outcome than patients with negative ALN status. 19 It matters clinician's decisionmaking of therapy that patients with positive ALN status are consider as high-risk patients and need to undergo adjuvant chemotherapy according to NCCN guideline. 20 Parts of TNBC patients with lager tumor burden and clinical positive ALN status are consider receiving neoadjuvant therapy. It should be noted that the misjudgment of the ALN status is a cause of patients' overtreatment and waste in health care. As for the patients with negative ALN status and without other risk factors, they could be treated by endocrine therapy alone, which means lower cost and more comfortable treatment experience. Thus, it's very important to accurately and clearly distinguish the ALN status for patients with early-stage breast cancer.
Besides MRI examination, ultrasound is also favored by clinicians in clinical practice due to the wide application and easily accessible data for study. Previous research demonstrated that an approach based on ultrasound radiomics also showed acceptable performance for ALN status prediction in early-stage breast cancer patients. 7,21 However, compared with MRI examination, operators' experience and heterogeneity from the angle of incidence of the ultrasound beam often affect the standardization of ultrasound image data, which result in the lack of reproducibility of model. As MRI examination possesses the advantage of showing the three-dimensional spatial position of each organ and tissue with high resolution, it could provide more comprehensive information for diagnosis and TNM staging, and therefore has been recommended for lesions di cult for the ultrasound and mammography and for preoperative estimation of breast cancer patients. 22 Additionally, the ALN-tumor signature based on MR radiomic in this study illustrated the superiority over general radiologists in detecting ALNM. Hence, it may be bene cial to combine MR radiomics with ultrasound radiomics, which might achieve better performance for ALNM prediction, especially improve the diagnosis and treatment level of hospitals in rural area and make up the disadvantages of areas with insu cient medical resources in the future.
Previous studies demonstrated that the tumor microenvironment harbored variety of immune cells, blood vessels and extracellular matrix, and the changes in the distribution of immune cells and blood vessel formation could facilitate the development and metastasis of tumors. 23 Therefore, many tumor microenvironment related factors like cancer-associated broblasts, myeloid-derived suppressor cells, and tumor associated macrophages have been used to predict tumor metastasis. [24][25][26] It has been reported that radiomic features could provide large amount of information of tumor microenvironment. 27 In this study, the ALN-tumor radiomic signature based on the integration of radiomic features from the tumor and the ALN regions could effectively predict ALNM, which also con rmed the correlation between MRI radiomics and the breast tumor microenvironment. Combination of radiomics with tumor microenvironment is promising to improve the predictive ability of ALNM.
Several limitations still need to be addressed in our study. First, the clinical information and MRI sequences of some patients were missing due to retrospective nature, which led to a decrease in sample size. Second, we failed to constructed a ALNM prediction model with the combination of tumormicroenvironment features, MR radiomics and ultrasound radiomics, but further exploration could be conducted in the future to realize more precisely prediction for ALNM.
In conclusion, this study presented a multiomic radiomic signature that incorporated MRI multi-sequence key radiomic features of ALN and tumor regions with clinicopathological characteristics and molecular subtype, which could be conveniently used for identifying ALNM patients among different molecular subtypes in early-stage invasive breast cancer and be promising to eventually result in a noninvasive approach to guide future clinical practice. Furthermore, the association between MRI radiomic features and tumor-microenvironment features including immune cells, long non-coding RNAs, and types of methylated sites were found, which revealed the potential biological underpinning of MRI radiomics. A total of 1,161 early-stage invasive breast cancer patients were recruited from four institutions in China, of which 1,088 patients passed quality control. The inclusion criteria included the followings: (a) female patients aged at least 18 years, who had histologically con rmed staged I-III invasive breast cancer; 29 (b) patients who had been treated with surgery and SLNB or ALND, and had been pathologically con rmed ALN status; (c) preoperative MRI scan of breast tumor and/or ALN were conducted, including contrastenhanced T1-weighted imaging (T1+C), T2-weighted imaging (T2WI), or diffusion-weighted imaging quantitatively measured apparent diffusion coe cients (DWI-ADC). The exclusion criteria included the followings: (a) patients underwent biopsy at an external institution and pathological results were not available; (b) patients suffering from other tumor diseases before or at the same time. The primary outcome was ALN metastasis status.  30 Patients with estrogen receptor (ER) and progesterone receptor (PR) positive, human epidermal growth factor receptor 2 (HER2) negative and Ki-67's level of < 14% breast tumors were de ned as Luminal A subtype patients. ER positive and HER2 over-expressed or ampli ed patients, or ER positive and HER2 negative patients with at least one of situations: Ki-67's level of > 14%, PR negative or low, were considered as Luminal B subtype patients. While ER and PR absent, HER2-positive subtype patients were distinguished by HER2 over-expressed or ampli ed, and HER2 negative breast tumors were classi ed as triple negative subtype. These biomarkers were measured with immunohistochemical methods or in situ hybridization.

Radiomic Features extraction
The protocol of multi-sequence MRI acquisition across four institutions and parameters of MR scanners for patients were shown in Supplementary Table 5. For all cohorts, multi-sequence MRI images from all cohorts were retrieved from Picture Archiving and Communication System, and radiomic features corresponding to the quantitative data obtained after computational translation of images were extracted from T1+C, T2WI, and DWI-ADC sequences imaging. Regions of interest of the breast tumor area (ROI-1), and ALN area (ROI-2) were semi-automatically delineated on each slice obtained via T1+C, T2WI and DWI-ADC (delineated with b value of 800 or 1000 s/mm2 and then copied to the corresponding ADC maps) sequences imaging by 3D Slicer software method (https://www.slicer.org/, version 4.10.2). 31 N4ITK Bias Field Correction were applied to obtained standard normal distribution of image intensities in ROI-1 and ROI-2, then the voxel-based features of ROIs were extracted using the SlicerRadiomics extension in 3D Slicer software, the in-house texture extraction platform developed based on the python package "PyRadiomics". A total of 5,718 quantitative radiomic features were extracted from ROI-1s and ROI-2s separately, including six groups of radiomic features: shape, rst-order, the gray-level cooccurrence matrix (GLCM), the gray-level size zone matrix (GLSZM), the gray-level dependence matrix (GLDM), and the neigbouring gray tone difference matrix (NGTDM). All patients were separately reassessed by two radiologists (Lu N and Li XH) blinded to the patients' clinical outcomes, MRI reassessed under the guidance of two senior radiologists (Xie CM and Wu Z) who major in MRI interpretation more than 20 years.
All ROIs of each MRI scan for each patient was normalization separately using Z-score to obtain a standard normal distribution of image intensities before radiomic features selection. The Random forest algorithm 32 was used to select the top 30 radiomic features from T1+C, T2WI, and DWI-ADC sequence of ROI-1, and ROI-2, respectively, in the training cohort.

Tumor and lymph node radiomic signature construction and validation
The signature were built based on the training cohort with MRI radiomic features that extracted from the breast primary tumor and the ALN regions of T1+C, T2WI, and DWI-ADC sequences. First, t-test was used to detect the associations between each feature and patients' ALN status. Next, to reduce the risk of bias and potential over tting, MRI radiomic features that achieved significance at P < 0.05 were entered into the further selection and top 30 features identi ed by random forest algorithm were selected from the each T1+C, T2WI, and DWI-ADC sequences for the construction of the signature via the support vector machine (SVM) algorithm. 33 The predictive accuracy of the ALN-tumor radiomic signature for predicting ALNM was initially assessed in the training cohort and then validated in the prospective-retrospective validation and external validation cohorts using receiver operating characteristic (ROC) curve analysis.

Multiomic radiomic signature construction and validation
The independent-sample t-test was used to assess the association between clinical characteristics, pathological characteristics, molecular subtype and ALN status in the training cohort. To provide the clinician a quantitative tool to predict individual probability of ALNM, all these characteristics with P < 0.05 and key MRI radiomic features were used for the construction of the multiomic radiomic signature via the SVM algorithm. ROC curve analysis was also performed to evaluate the performance of the multiomic radiomic signature in the training cohort and was validated in the prospective-retrospective and the external validation cohorts.

MRI radiomics features associated with tumor microenvironment
The t-test was utilized through to identify differentially expressed genes 34 associated with the radiomic score in 91 patients with T1+C and T2WI sequences MRI from TCGA and TCIA. The GO and KEGG analysis were performed using the clusterPro ler R package 35 . The GO terms and KEGG pathways were considered statistically signi cant with P values and false discovery rates less than 0.05.
The CIBERSORT algorithm 36 and the LM22 gene signature were used for highly sensitive and speci c discrimination of 22 human immune cell phenotypes. CIBERSORT is a deconvolution algorithm that uses a set of reference gene expression values (a signature with 547 genes) that is considered a minimal representation for each cell type. Based on those values, CIBERSORT infers cell type proportions in data from bulk tumor samples with mixed cell types using support vector regression. Gene expression pro les were prepared using standard annotation les, the data were uploaded to the CIBERSORT web portal (http://cibersort.stanford.edu/), and the algorithm was run using the LM22 signature at 1,000 permutations.
12,578 long non-coding RNAs data of transcriptome RNA sequencing based on Illumina platform and 23,381 types of methylated sites data from both Illumina Human Methylation 27 and Illumina Human Methylation 450 platform were used for analysis. The top 30 ALNM-associated long non-coding RNAs and types of methylated sites were selected using the random forest algorithm. Correlation analysis and weighted linear regression models were used to estimate the strength of the correlations with Pearson ρ.
Statistical analysis.
The χ2 test was performed to examine the differences in categorical variables, and the independent t-test was applied to compare the differences in continuous variables between two groups. The predictive accuracy of the signatures was assessed by using ROC analysis. The area under ROC curve (AUC) was used to evaluate the performance of sensitivity and speci city in each signature. Additionally, DCA was performed to assess the clinical utility of the prediction model by quantifying the net bene ts when different threshold probabilities were considered. 37 Patients were categorized into high-and low-score groups with the optimal cutoff values de ned by the R package cutpointr. Heatmaps were generated to show the distribution and expression levels of radiomic features, immune cells, long non-coding RNAs and types of methylated sites. Survival was calculated using the Kaplan-Meier method and the log-rank test, and hazard ratios (HRs) and 95% con dence intervals (Cls) were calculated using a Cox regression analysis. For all the analyses, two-sided P-values less than 0.05 were considered statistically signi cant. Statistical analyses were performed using R software (version 4.0.2). This study is registered with ClinicalTrials.gov, number NCT04003558. Drafting of the manuscript: All authors.
Critical revision of the manuscript for important intellectual content: All authors.
Statistical analysis: All authors.  Abbreviations: ALN, axillary lymph node; PPV, positive predictive values; NPV, negative predictive values; CI, con dence interval; AUC, area under the receiver operating characteristics curve.