This study was approved by our institutional review board and written informed consent was obtained from all participants (Approval Nos. 2020AS0113 and 2021AS0318). This prospective study was registered at clinicaltrials.gov (NCT06095414). Sample size estimation was described in the Supplementary Information. Between June 2020 and April 2022, we enrolled 215 consecutive women who were scheduled to undergo ultrasound-guided tissue biopsy for suspicious breast masses assessed as category 4C or 5 according to the Breast Imaging-Reporting and Data System (30). If a patient had multiple suspicious masses, the most suspicious mass was selected as the representative lesion before biopsy. Ultrasound-guided 14-gauge core-needle biopsy (Bard, Tempe, Arizona, USA) was performed, including the boundary of the tumor. The eligibility criteria were as follows: (a) histologically confirmed invasive breast cancer, (b) participant consent, and (c) MRI performed before treatment. We excluded 94 women for the following reasons: benign histopathology results (n = 38), ductal carcinoma in situ (n = 19), invasive tumor size < 10 mm (n = 8), infiltrative nonmass lesion on MRI (n = 7), refusal to participate (n = 10), no MRI before treatment (n = 4), inadequate volume histological imaging (n = 3), and MRI with a different protocol at an external hospital (n = 5). Ultimately, 121 women were included in this study (Figure 1).

We performed FxClear volume histology, which is an acrylamide-free electrophoretic tissue-clearing protocol with a fast reaction time and high immunoreactivity (31). Volume imaging preserves the structural integrity and allows visualization of the biological architecture in thick tissues and organoids (32, 33). The immunohistochemical reactivity of TME components of each cancer was quantified. Primary antibodies were used to assess the following TME components: collagen type 1 (1:4500, AB34710, Abcam), laminin (1:4500, L9393, Sigma), nidogen-1 (1:9000, NBP1-97701, Novus), FOXP3 (1:600, AB20034, Abcam), and CD8 (1:25, AB75129, Abcam). FOXP3 antibody was used to evaluate regulatory T cells, and CD8 was used to evaluate cytotoxic T cells. The expression of ECM components, including collagen, laminin, and nidogen, was measured as the amount of immunoreactivity. The expression of regulatory and cytotoxic T cells was assessed by counting the number of cells (Figure 2). Details of the histological evaluation of the TME are provided in the Supplementary Information. We also reviewed the histological reports to evaluate prognostic factors and subtypes. The histological factors were dichotomized by respective median value unless specified: hormone receptor status including estrogen or progesterone (positive vs. negative), HER2 (positive vs. negative), and Ki67 (high [> 20%] vs. low [≤ 20%]). The molecular subtype was classified into four types according to the St Gallen classification criteria: luminal A, luminal B, HER2-enriched, or triple-negative (34).

MRI was performed with the patient in the prone position using a 3T magnet (MAGNETOM Vida; Siemens Healthineers, Erlangen, Germany) with a dedicated 18-channel breast coil. 3D T2, 3D dynamic contrast-enhanced T1, and DWI were obtained in each patient. Table 1 demonstrates acquisition parameters for MRI. All images were obtained with bilateral axial views and fat saturation. For dynamic contrast-enhanced T1, five postcontrast phases were acquired at 60, 120, 180, 240, and 300 seconds after intravenous contrast administration. DWI was performed using the intravoxel incoherent motion technique.

For radiomics, we analyzed five MRI inputs: T2, DWI (b = 800 s/mm²), precontrast T1, initial-phase postcontrast T1 (60 s), and delayed-phase postcontrast T1 (300 s). In this manuscript, ‘phase’ denotes the two time points (initial and delayed) of the same dynamic contrast-enhanced T1 sequence. Although five DCE time points were acquired (60, 120, 180, 240, 300 s), we pre-specified 60 s and 300 s as the representative phases for radiomics because they correspond to the Breast Imaging-Reporting and Data System (BI-RADS) MRI lexicon windows (initial enhancement vs. delayed persistent/plateau/washout) (35). The initial-phase emphasizes perfusion/permeability, whereas the delayed-phase reflects delayed contrast retention. Intermediate phases (120–240 s) were used only to derive BI-RADS kinetic labels for descriptive reporting and were not included as radiomic inputs to avoid multicollinearity and unnecessary dimensional inflation.

For 3D segmentation of tumors, regions of interest (ROIs) along the entire enhancing tumor margin of cross-sectional area were drawn at axial views of initial-phase contrast-enhanced T1 from top to bottom of each tumor. On precontrast and delayed-phase contrast-enhanced T1, ROIs were drawn on the lesion corresponding to initial-phase contrast-enhanced T1. Because the slice thickness differed for DWI (4 mm) and T2 (3 mm) relative to T1 (1 mm), ROIs for DWI and T2 were drawn separately. Segmentation was performed on axial images by a radiologist (E.S.K. with 5 years of experience in breast MRI) using a semi-automated method with MRIcro software (version 1.40, https://www.nitrc.org/projects/mricro/) under the supervision of a senior radiologist (B.K.S. with 23 years of experience in breast MRI), while blinded to the clinicohistological data except for information about the diagnosis of breast cancer (Figure 2). The peritumoral regions were generated by extending the intratumoral regions with 4 mm thickness in 3D. A distance of 4 mm from the tumor was chosen to evaluate the peritumoral area based on previous breast peritumoral radiomics studies (22, 23). In addition, we performed a margin-sensitivity analysis using an 8-mm annulus on the best-performing input from the main analyses to assess robustness.

To assess intra-reader segmentation reproducibility, we performed a stratified random sampling of 50 cases according to tumor-size distribution; each case was re-segmented by the same reader using the identical semi-automated protocol and compared with the original segmentation. Agreement between the two segmentations by the same reader was quantified using the Dice similarity coefficient (DSC) and the Jaccard similarity coefficient (JSC) for volumetric overlap (36).

For feature robustness, isotropic resampling was performed prior to feature extraction. For each of the five MRI inputs, 1618 features were extracted from the intratumoral region and combined intratumoral and peritumoral regions. Thus, a total of 16180 radiomic features were extracted for each cancer. The 1618 features were categorized into four groups: (a) histogram-based first-order statistical features (n = 17), (b) shape and volume features (n = 7), (c) textural features using the gray-level co-occurrence matrix (GLCM) and gray-level run-length matrix (GLRLM) (n = 162), and (d) wavelet-transformed features (n = 1432) (26). We extracted MRI radiomic features in compliance with IBSI (27, 37) using MATLAB R2023b (MathWorks) and the Pyradiomics 3.1.0 library (https://www.radiomics.io/pyradiomics.html) in Python 3.8. The radiomics quality score was 25 out of 36 (69%) (Supplementary Table 1) (25). By context, prior oncology radiomics studies have reported a mean radiomics quality score of 9/36 (26%). Details of feature extraction are provided in the Supplementary Information.

After z-score normalization of all features, the dataset was randomly divided into training and testing cohorts in a 7:3 ratio. For each MRI input independently, feature selection was performed using the least absolute shrinkage and selection operator (LASSO) method with fivefold cross-validation to determine the optimal lambda (0.05). The LASSO tool is a statistical method used to select the most relevant features from complex data sets. The LASSO process was iterated 25 times, and features that were selected 20 times during each iteration were considered significant radiomic features. After the process iterations, we identified two sets of five top features from intratumoral features alone and from both intratumoral and peritumoral features and then used these radiomic signatures as predictors in a linear regression analysis (38, 39). Beyond the primary 7:3 split, we performed fivefold cross-validation to assess internally stability and robustness of the predictive models. We report the mean area under the receiver-operating characteristic curve (AUC) and 95% confidence intervals (CI) across folds as the cross-validated performance. Ultimately, we constructed two predictive models: intratumoral-only model based on intratumoral radiomic signatures and a combined intratumoral–peritumoral model combing both intratumoral and peritumoral signatures. To aid interpretability, we visualized the relative contributions of the top 10 selected features as bar plots of standardized coefficients with bootstrap 95% CIs, and assessed redundancy or complementarity using a Pearson correlation heatmap.

To assess potential imbalance between the training and testing cohorts, we compared baseline clinicopathological characteristics: age, tumor size, hormone receptor status, HER2 status, Ki67, and molecular subtype. For continuous variables, normality was assessed with the Shapiro–Wilk test, and between-group differences were evaluated using Welch’s t test or the Mann–Whitney U test, as appropriate. For categorical variables, chi-square or Fisher’s exact tests were applied. We also calculated effect sizes (Cohen’s d for continuous variables; absolute differences in proportions for categorical variables), interpreting 0.2, 0.5, and 0.8 as small, medium, and large effects, respectively (40).

In addition, we examined the associations between MRI radiomic features and histological TME components. Model performance was assessed using AUC, accuracy, sensitivity, specificity, positive predictive value, and negative predictive value. To assess the model’s discriminative performance, the AUC was calculated with its 95% CI determined using the bootstrap method, focusing on the 2.5% and 97.5% percentiles of the distribution. The DeLong test was conducted to compare the performance of models using only intratumoral radiomic signatures with models that combine intratumoral and peritumoral radiomic signatures. Additionally, test compared predictive performance across five MRI sequences. A P value < 0.01 (0.05/5), 0.0125 (0.05/4), or 0.025 (0.05/2) from the DeLong test was considered significant, applying Bonferroni correction for multiple comparisons. All reported performance metrics were derived from confusion matrices of our cross-validated predictions.

We evaluated clinical utility using decision-curve analysis. For each model, net benefit was computed across prespecified threshold probabilities (Pt = 0.10–0.40) and compared with treat-none and treat-all strategies. We compared intratumoral-only and combined intratumoral- peritumoral models for two representative endpoints: collagen (ECM) and regulatory T cells (immune cells), using the best performing sequence.

Furthermore, we performed subtype-stratified analyses using the MRI sequence that showed the best overall predictive performance in the main analyses. Within each subtype, intratumoral-only and combined models were evaluated, and pairwise AUC comparisons were conducted using DeLong’s test with Bonferroni correction. To compare performance between subtypes, we used bootstrap resampling (1,000 iterations) to estimate P values and CIs, accounting for unequal and relatively small sample sizes.

All statistical metrics were computed using scikit-learn in Python 3.8 (41). The overall design of the study is presented in Figure 3. Analysis code is available at [ https://github.com/KUAH-rad/Radiomics-analysis].

We evaluated samples from 121 women (mean age, 54 ± 11 years; age range, 30–83 years) with invasive breast cancer. Of the 121 cancers, 114 were invasive ductal carcinoma (Figure 2), five invasive lobular carcinoma, one medullary carcinoma, and one metaplastic carcinoma. Table 2 summarizes baseline characteristics, and Supplementary Figure 1 compares the training and testing cohorts. No significant differences were observed for age, tumor size, estrogen receptor, progesterone receptor, hormone receptor, HER2, Ki67, or molecular subtype (all P > 0.05). Effect-size analyses indicated negligible-to-small differences across all variables, with the largest for age (Cohen’s d = 0.188) and Ki-67 (d = 0.176). Additionally, standardized mean differences for the top 10 radiomic features in the best-performing models—predicting collagen (ECM representative) and regulatory T cells (immune representative)—showed no material imbalances between the training and testing cohorts. (Supplementary Figure 2).

Tables 3, 4 demonstrate the performance of models using intratumoral radiomics signatures and combing intratumoral and peritumoral radiomic signatures from each MRI sequence. Because model performance ultimately depends on reliable ROI definition, we evaluated intra-reader segmentation reproducibility in 50 randomly selected cases; volumetric overlap was excellent (mean DSC = 0.97 ± 0.01, mean JSC = 0.95 ± 0.02). The radiomic signatures of intratumoral regions included 60 first-order and 65 texture features (62 GLCM features and three GLRLM features), and those of combined regions included 59 first-order, 53 texture features (52 GLCM features and one GLRLM feature), and 13 shape and volume features (Supplementary Tables 2, Supplementary Tables 3). Among these signatures, the most frequently selected feature for intratumoral regions was the first-order standard deviation and in the combined regions, it was the first-order energy. These two features were found to increase in cancers associated with higher ECM and regulatory T cells, and lower cytotoxic T cells. These conditions, when the TME components were dichotomized based on the median, correlated with a poor prognosis (P < 0.05) (Supplementary Table 4).

Table 5 shows the comparison of AUCs between two models using intratumoral radiomic signatures alone and combined intratumoral and peritumoral radiomic signatures to predict TME. Models using combined radiomic signatures showed significant improvements in predicting laminin on delayed postcontrast T1 (AUC [95% CI] for combined model vs. intratumoral model, 0.83 [0.78, 0.88] vs. 0.80 [0.76, 0.85]), nidogen on DWI (0.73 [0.68, 0.78] vs. 0.68 [0.64, 0.72]), regulatory T cells on DWI (0.69 [0.64, 0.74] vs. 0.66 [0.63, 0.69]), regulatory T cells on T2 (0.82 [0.77, 0.87] vs. 0.77 [0.72, 0.81]), and cytotoxic T cells on precontrast T1 (0.82 [0.76, 0.87] vs. 0.79 [0.74, 0.83]) (P < 0.01).

The bar plot of standardized coefficients highlighted a subset of top features with the greatest contributions (Supplementary Figure 3). The Pearson correlation heatmap demonstrated clustered groups of highly correlated features and low inter-feature correlations (Supplementary Figure 4).

When comparing performance between five MRI sequences in models using combined intratumoral and peritumoral radiomic signatures, delayed postcontrast T1 yielded the highest AUCs (AUC [95% CI] for collagen 0.82 [0.79, 0.86], laminin 0.83 [0.78, 0.88], nidogen 0.83 [0.76, 0.89], regulatory T cells 0.83 [0.78, 0.88], and cytotoxic T cells 0.83 [0.79, 0.89]) (Table 4). The AUCs of delayed postcontrast T1 were significantly superior to those of DWI in predicting all the TME components, precontrast T1 in predicting collagen, and T2 in predicting nidogen (P < 0.0125) (Table 6). However, there was no difference in AUCs between delayed postcontrast T1 and initial postcontrast T1 for all TME components (P ≥ 0.0125). The combination of initial postcontrast T1 and delayed postcontrast T1 did not show a significant improvement over using either initial or delayed postcontrast T1 alone (P ≥ 0.025) (Supplementary Table 5).

Using fivefold cross-validation, the AUCs of the overall performance in predicting the abundance of ECM and immune cells using combined features were 0.82 [0.78, 0.87] and 0.82 [0.78, 0.88] on initial postcontrast T1 and 0.82 [0.77, 0.87] (Figure 4) and 0.83 [0.78. 0.88] on delayed postcontrast T1, respectively (Figure 5). As shown in these figures, the “Actual AUCs” listed in the figure panels represent the values obtained from each fold of the fivefold cross-validation. Table 7 shows radiomic signatures for predicting overall ECM and immune cell abundance on initial or delayed postcontrast T1. Nineteen (95%) out of 20 radiomic signatures were first-order and GLCM texture features.

Decision-curve analysis demonstrated that the combined intratumoral-peritumoral models provided consistently greater net benefit than intratumoral-only models across the clinically relevant threshold range (P_t = 0.10-0.40) (Supplementary Figure 5). For collagen, the combined model outperformed the intratumoral-only model across all threshold, and both exceeded the treat-all and treat-none strategies. Similarly, for regulatory T cells, the combined model consistently yielded higher net benefit than the intratumoral-only model within the same threshold range. In the combined models on delayed-phase postcontrast T1, we performed correlation analyses between the top radiomic features and quantitative histology metrics, collagen (ECM) and regulatory T cells (immune cells). GLCM contrast was significantly correlated with collagen proportion, and first-order standard deviation was correlated with regulatory T-cell density (both P < 0.001) (Supplementary Figure 6).

In subtype-stratified analyses using the delayed postcontrast T1-weighted sequence, luminal A cancers showed AUCs of 0.59-0.61 for intratumoral-only models and 0.64-0.67 for combined intratumoral–peritumoral models (Supplementary Tables 6, Supplementary Tables 7). Luminal B cancers showed AUCs of 0.59-0.60 for intratumoral-only models and 0.62-0.65 for combined models. Triple-negative breast cancer showed AUCs of 0.56-0.59 for intratumoral-only models and 0.58–0.64 for combined models. Across subtypes, combined models tended to outperform intratumoral-only models for TME prediction, but differences did not reach statistical significance after correction for multiple comparisons (P > 0.016) (Supplementary Table 8). In the combined models, ECM prediction for luminal A cancers was higher than for triple-negative cancers (collagen P = 0.031; nidogen P = 0.28); however, these pairwise differences did not meet the Bonferroni-corrected significance threshold (P = 0.025 for two comparisons, 0.05/2) (Supplementary Table 9). HER2-enriched tumors (n = 4) were excluded from statistical testing because such a limited cohort renders receiver-operating characteristic curve estimation unstable and prevents reliable inference.

Additionally, we performed a peritumoral margin-sensitivity analysis (4 mm vs. 8 mm) on delayed postcontrast T1. the 8 mm annulus yielded lower AUCs, supporting 4 mm as the primary peritumoral margin (Supplementary Table 10).