This study was approved by the Ethics Committee of the Second Affiliated Hospital of South China University of Technology (Guangzhou First People’s Hospital) Hospital, with informed consent being waived due to the retrospective nature of this study. A total of 535 patients who underwent breast mpMRI for preoperative assessment at our hospital between January 2017 and April 2022 were included. The inclusion criteria were as follows: (1) histopathological confirmation of BC by surgical resection or needle biopsy; (2) patients who underwent a routine mpMRI including T1WI, T2WI, DWI (with b values of 0 s/mm², 600 s/mm² and 800 s/mm²), DWI-derived ADC map and DCE-MRI (with 6 continuous enhancing phases) within one week prior to pathological examinations; (3) no additional therapy prior to MRI. The exclusion criteria were: (1) recurrent BC (n=11); (2) incomplete pathological results, such as those lacking IHC results and Ki-67 scores, or unclear histological types (n=15); (3) cases with Volumes of interest (VOI) that were difficult to delineate due to images artifacts (n=39). (4) patients with breast implants (n=4). In cases of multicentric or multifocal tumors, only the largest malignant lesion was selected. For bilateral disease, the largest lesions of both breasts were selected according to pathological results. Finally, 460 patients with 466 lesions were enrolled in this study. The lesions were categorized into HR+ (n=336) and HR- groups (n=130), with the HR- group further divided into HEBC (n=76) and TNBC (n=54) subgroups. Based on sample size calculations (15, 16), a required sample size of 210 (42 cases of TNBC and 168 cases of non-TNBC) was sufficient to detect differences between various molecular subtypes of BC with a power of 95%. Appendix 1 showed the detailed information on the sample size calculation process. All lesions were divided randomly into a training/validation cohort (n=337) and a test cohort (n=129) at a ratio of ~3:1, in which a random selected training cohort (n=337) was established to determine the optimal single MR sequence for subsequent experiments, as shown in Figure 1 .

All preoperative routine mpMRI examinations were performed on a 1.5-T MRI system (uMR 560, United Imaging) using a dedicated 4-channel SENSE breast coil. The standard examinations included T1W, T2W, DWIs (with three b values of b=0 s/mm², b=600 s/mm² (DWI₆₀₀) and b=800 s/mm² (DWI₈₀₀)) and six continuous DCE-MRI (DCE_1-6) scans with ~62 seconds per phase. ADC maps were derived from DWI using the b=0 s/mm² and b=800 s/mm² data through the embedded immediate post-processing software. The DCE-MRI protocol involved injecting gadolinium-diethylenetriamine pentaacetic acid (GD-DTPA, 0.1 mmol/kg) and acquiring images in the five consecutive phases after pre-contrast T1WI. The six phases were named as: pre-contrast (DCE₁), super-early-contrast (DCE₂), early-contrast (DCE₃), and delayed-contrast (DCE_4-6). The detailed MRI scanning parameters are provided in Supplementary Materials Table S1 .

The volume of interest (VOI) was defined on all images that were stored in DICOM format. In order to standardize the extracted image biomarkers from mpMRI, we followed the major procedure outlined by the Image Biomarker Standardization Initiative (IBSI) (17). Before VOI delineation, we used the General registration (elastix) method, available as the “SlicerElastix” plugin in the open-source image analysis platform 3D Slicer (https://www.slicer.org), to register all sequences’ images. This alignment enabled us to better handle morphological variations and structural differences in breast tissue, particularly when aligning the other sequence images with the DCE₂ image. Additionally, we resampled all MRI sequences to a standard resolution of 1.096 x 1.096 x 1.2, ensuring isotropic voxels and reducing variations caused by differences in scanning equipment, protocols, and patient positioning. Furthermore, we normalized the intensity levels of all images to a range of 0-255 to reduce the influence of contrast and brightness variations, which might otherwise affect the quantification of radiomics features (18).

Slice-wise delineation of the VOI was carried out using the ITK-SNAP software (http://www.itksnap.org) on T2W, DWI₆₀₀, DWI₈₀₀, ADC maps, and DCE_1-6 images. The process started with manual delineation of the visible tumor margins on the DCE₂ images with the most distinguishable lesion boundary. The contoured VOIs on DCE₂ were then replicated to the remaining DCE sequences, resulting in 6 VOIs based on the DCE-MRI data. Similar steps were repeated for the DWI₆₀₀, DWI₈₀₀, and ADC maps. Subsequently, VOI delineation was performed on the T2W images according to the position and shape of the VOI completed above. Slight adjustments were allowed on all images in order to obtain tailored VOIs. The VOI delineation was performed by two radiologists (WZ and JL, with 6 and 4 years of experience in radiological diagnosis, respectively) who were blind to all prior patient information. The interobserver correlation coefficient (ICC) value of the two radiologists was assessed.

The radiomics features were extracted from ten VOIs of each lesion using the open-source software toolkit Pyradiomics (19). A total of 109 features were extracted from three categories of features: 1) intensity features (n=19); 2) morphology features (n=15); texture features (n=75). Only the extracted radiomics features with ICC > 0.75 were then fed into 150 classification base models, which were built using 10 classifiers and 15 feature selection methods. Detailed definitions of the above-mentioned features can be found in Pyradiomics documentation and IBSI (17). The full list of radiomics features and the methods employed in this study are summarized in Tables S2 , S3 , respectively.

Based on the newly developed mpMRI-based RadioFusionOmics model by our lab, we constructed a feature fusion radiomics (R_FF) model that integrated radiomics information from different MRI sequences to produce more discriminative fused features. A random selected training cohort (n = 337) was used to analyze all radiomics features from each MRI sequence, analogous to a radiologist’s initial reviewing of a patient’s complete set of MR images. According to the highest cross-validation AUC obtained in the training/validation process, the optimal single sequences that can identify hormone receptor positive (HR+) vs. HR- BC, TNBC vs. HEBC, as well as TNBC vs. non-TNBC were determined and regarded as the single sequence-based radiomics (Rss) model.

Subsequently, the radiomics features from the top four high-performing single sequences were combined to perform multiple sequence feature fusion, similar to a radiologist’s final reviewing focusing on specific sequences after a preliminary review. The best combination of sequences (combination of two, three or four sequences, a total of 11 types of combinations) was then identified to develop the R_FF models. Utilizing feature-level fusion, the R_FF model conducted a feature-wise fusion strategy by finding a transformation to map the feature matrix with a set of MRI sequences (e.g., dimension = 10) to a lower dimensional space (e.g., dimension = 1). By integrating the class structure information (i.e., information on the molecular receptor status of memberships of the training samples) in the calculation of the transformation, the R_FF was able to eliminate the between-class correlations and strengthen the within-class correlations during the feature fusion, which can effectively enhance the discriminative power of fused features. Various base models (n=11*150 = 1650) were trained using the fused features and their performances were evaluated and ranked via a stratified ten-fold cross-validation. The optimal base models for Rss and R_FF were verified in the training/validation cohort and test cohort. Technical details related to the R_FF are shown in Appendix 2 . The flow chart of this study was displayed in Figure 2 .

All surgical or biopsy specimens were examined by two pathologists (YZ and WD, with 6 and 16 years of experience in the pathological diagnosis of BC, respectively). The following pathological biological markers of BCs were assessed and recorded: tumor maximal diameter, affected side in the breast, number of tumors, histology type, and IHC status of estrogen receptor (ER), progesterone receptor (PR), HER-2, and Ki-67 index. Tumors with ER or PR positive expression (> 10% of tumor nuclei staining) were classified as HR positive (HR+) (20). Positive HER-2 expression was defined as a 3+ IHC score or 2+ accompanied by fluorescence in situ hybridization positive (FISH+) result (21). The Ki-67 scores were classified into two groups: < 14% as low Ki-67 level and ≥14% as high Ki-67 level. The molecular subtypes of BCs were classified as follows: luminal A (ER and/or PR positive, HER-2 negative, and Ki-67 < 14%), luminal B (ER and/or PR positive, HER-2 negative, and Ki-67 ≥ 14% or ER and/or PR positive, HER-2 positive, regardless Ki-67 expression), HER-2 enriched (ER and PR negative, HER-2 positive), which was recorded as HEBC, and triple negative cancer (ER, PR and HER-2 negative), named as TNBC. The luminal A and luminal B comprised the HR+ group. The Ki-67 expression was scored as the percentage of positive invasive tumor cells with any nuclear staining, with the mean percentage of positive cells recorded (4). Four cases of different molecular subtypes of breast cancer were presented in the supplementary materials Figures S3 - S6 .

The Chi-square Test and Fisher’s Exact Test were used for categorical variables, the One-way ANOVA analysis was used for normally distributed continuous variables, and the Kruskal-Wallis H test was used for non-normally distributed continuous variables to compare demographic and pathological characteristics between different molecular subtypes. The normality of data distribution was evaluated by the Shapiro-Wilk test. The results for normally distributed continuous variables were reported as mean ± SD, while non-normally distributed continuous variables were reported as median (interquartile range, IQR). Categorical variables were presented as numbers and proportions. The performance of each Rss and R_FF base models were evaluated via the area under the receiver operative characteristic curve (AUC), sensitivity (SEN), specificity (SPE) and accuracy (ACC) among different subtypes of BC. The performance of the Rss and R_FF was compared using the paired samples Wilcoxon signed rank test. Two-sided p < 0.05 was considered statistically significant. All statistical analyses were conducted using the SPSS 25.0 software (IBM SPSS Corporation, USA) and python 3.6.2 (Python Software Foundation (USA, https://www.python.org/downloads/).

The clinical pathological characteristics of the 460 patients with 466 lesions (6 patients had bilateral lesions) enrolled in the study are presented in Table 1 . Among the 466 lesions, 336 lesions (72.1%) were classified as HR+ BCs, with 142 lesions being luminal A and 194 lesions being luminal B. Additionally, 76 lesions (16.3%) were classified as HEBCs, and 54 lesions (11.6%) were classified as TNBCs. The median tumor size of TNBCs (26.0 mm) and HEBC (27.0 mm) was found to be significantly larger than that of HR+ (21.0 mm) (p = 0.000). TNBCs showed a higher prevalence of mass enhancement in DCE MRI (81.5%) and invasive carcinoma (96.2%) compared to HR+ and HEBCs (p < 0.001). TNBCs also had a higher Ki-67 index (> 14%) in comparison with HR+ and HEBCs. Moreover, the age of patients and number of tumors among HR+, HEBC and TNBC groups were significantly different (p < 0.05). Baseline characteristics were not significantly different between both training/validation and test cohorts ( Table S4 ).

All the discriminative base models established based on single mpMRI sequence were compared to determine the optimal sequences among HR+ vs. HR-, TNBC vs. HEBC and TNBC vs. non-TNBC. Supplementary Figure S1 demonstrated the discrimination comparison results on ten sequences of the three classification tasks.

By analyzing the dominant radiomics features of each sequence, the optimal sequence for discriminating HR+ vs. HR- was DWI_600, the optimal Rss model, namely Rss (DWI₆₀₀), achieved the highest AUC of 0.787 in the random training cohort ( Figure 3 ), and similar performance in the training/validation cohort (AUC=0.767) and test cohort (AUC=0.768), respectively ( Table 2 ).

The optimal sequence for identifying TNBC and HEBC was DWI-derived ADC map, the best Rss model, recorded as Rss (ADC), yield the highest AUC of 0.788 in the random training cohort ( Figure 3 ), and the best AUC of 0.769 and 0.718 in the training/validation cohort and test cohort, respectively ( Table 2 ).

Regarding TNBC vs. non-TNBC discrimination, the ADC map was also the best sequence, the optimal Rss model (Rss [ADC]) demonstrated the highest AUC of 0.809 in the random training cohort ( Figure 3 ), and the best AUC of 0.784 and 0.735 in the training/validation cohort and test cohort, respectively ( Table 2 ).

We selected the top four superior sequences for molecular receptor status classification to build the R_FF model. As shown in Figure 3 and Figure S1 , the top four superior sequences for HR+ vs. HR- were DWI₆₀₀, DWI₈₀₀, DWI-derived ADC map and DCE_5, with all AUCs > 0.77 in the random training cohort ( Figure S1A ). Similarly, DWI-derived ADC map, DCE₂, DCE₃ and DCE₄ were the top four dominant sequences for TNBC vs. HEBC, yielding all AUCs > 0.72 ( Figure S1B ). While the four most predominant sequences for TNBC vs. non-TNBC were DWI-derived ADC map, DWI₆₀₀, T2WI and DCE₂, achieving all AUCs greater than 0.73 ( Figure S1C ).

Subsequently, the performances of each combination of the top two, three or four high-performance mpMRI sequences in random training cohort (a total of 11 types of combinations during each classification task) were compared and displayed in Figure S2 . Our results illustrated that the model R_FF (DWI₆₀₀+DWI₈₀₀+DCE₅), R_FF (ADC+DCE₂+DCE₄) and R_FF (ADC+DWI₆₀₀+T2WI+DCE₂) were superior over the other sequences combinations in the random training cohort, yielding the maximal AUC of 0.809, 0.805 and 0.847, respectively. Similar performances were obtained in the training/validation cohort and test cohort, outperforming the Rss model with an AUC of 0.778 and 0.726, 0.787 and 0.773, 0.818 and 0.773, respectively (both p<0.05 except HR+ vs. HR-), as shown in Table 2 . Among R_FF (DWI₆₀₀+DWI₈₀₀+DCE₅), R_FF (ADC+DCE₂+DCE₄) and R_FF (ADC+DWI₆₀₀+T2WI+DCE₂), the base model (classifier + feature selection method) were respectively “Logistic Regression + Multi-Cluster Feature Selection” (MCFS), “Logistic Regression + Discriminative Feature Selection” (UDFS) and “Logistic Regression + trace_ratio”. The MpMRI-based feature fusion method employed in the task of TNBC vs. non-TNBC achieved the optimal discriminative capability, yielding AUC, ACC, SEN and SPE of 0.818, 0.718, 0.705, 0.721 in the training/validation cohort and 0.773, 0.767, 0.636, 0.780 in the test cohort, respectively.

The top-ranked features associated with the three classification tasks were also sieved by the proposed R_FF model and their discriminative capabilities were analyzed. Based on the feature selection procedure of each base model, we counted and ranked the occurrence of each selected feature (only for base models with AUC > 0.6). The fifteen most frequently selected features of the three classification tasks were displayed in Tables S5 - S7 . Most dominant features were texture features in HR+ vs. HR- (8/15) and TNBC vs. HEBC (8/15), while intensity-based features were the superior discriminative features of TNBC vs. non-TNBC (11/15). The top 5 most frequently selected radiomics features associated with the discrimination of HR+ and HR- included three morphology-based features and two gray level co-occurrence matrix (GLCM) features, while intensity-based features accounted for 80% (4/5) and 100% (5/5), respectively among the top 5 radiomics features of TNBC vs. HEBC and TNBC vs. non-TNBC ( Table 3 ). All the features showed statistically significant differences between HR+ and HR-, TNBC and HEBC, TNBC and non-TNBC with p-values < 0.001. The mean feature values of each group were used as the threshold to identify different molecular receptor statuses. In the task of discriminating TNBC from non-TNBC, the top 5 features outperformed other two tasks, with ~75% of the non-TNBC having larger feature values, while ~65% of the TNBC group had smaller values in all top 5 features ( Table 3 ).