Predicting Human Epidermal Growth Factor Receptor 2 Expression in Breast Cancer Based on Radiomics of MRI Habitat and US

Introduction

Breast cancer is one of the most common cancers worldwide.¹ Overexpression of human epidermal growth factor receptor-2 (HER-2) is associated with high-grade tumors and lymph node metastasis.² Trastuzumab, the first of these targeted treatments, significantly improves survival rates and prognosis.³ Several targeted therapies have also emerged for specific molecules and pathways in triple-negative breast cancer, which have shown promising results in clinical trials.⁴ Consequently, assessing HER-2 status is crucial for selecting the appropriate treatment and predicting prognosis in breast cancer patients.

Currently, the gold standards for evaluating HER-2 expression are immunohistochemistry (IHC) and fluorescence in situ hybridization (FISH).⁵ However, both methods require invasive procedures to obtain pathological tissue samples.⁶ Additionally, HER-2 status may change in some patients during or after neoadjuvant chemotherapy, and such alterations more frequently manifest as a loss of HER2 expression (occurring in 18–42% of cases), which may be associated with a higher recurrence rate and a worse prognosis.⁷ Repeated assessments of HER-2 status during treatment are necessary to adjust medication promptly, but the invasiveness of biopsies limits their use.

Ultrasound (US) and magnetic resonance imaging (MRI) are commonly used in breast cancer evaluation.^8,9 While these imaging techniques can reveal features related to molecular subtypes of breast cancer, including HER-2 status, their predictive value based on visual markers remains limited.^10,11 Radiomics can enable more efficient data analysis through extracting quantitative features from images.¹² It has been used in various areas of breast cancer research, such as distinguishing between benign and malignant tumors, evaluating molecular subtyping, and predicting treatment efficacy and prognosis.^13–16 However, most radiomics-related studies currently focus on the entire tumor region, assuming a uniform distribution of heterogeneity within the tumor.^17,18

Recently, “habitat” imaging has emerged as a method for dividing tumors into subregions by identifying grayscale voxels, offering a more refined characterization of tumor heterogeneity. Studies have shown that specific tumor habitats are related to drug resistance and poor prognosis.^19,20 To our knowledge, few studies has combined multi-parametric MRI habitats, US, and clinical factors to model and predict HER-2 status in breast cancer. In this research, we developed a combined model based on multi-parametric MRI habitat radiomics, US imaging radiomics, and clinical features to predict HER-2 expression status in breast cancer.

Material and Methods

Patients

The patient inclusion process is shown in Figure 1. This retrospective study was approved by the First Affiliated Hospital of Guangxi Medical University Ethical Review Committee (No. 2024-E530-01) with a waiver of informed consent, as the research involved only anonymized medical records and posed no additional risk to participants. All data were de-identified and handled in strict compliance with the Declaration of Helsinki to protect patient confidentiality.

Figure 1 The patient recruitment in our study. Abbreviation: HER-2: Human epidermal growth factor receptor-2.

We included breast cancer patients from May 25, 2019 to April 15, 2025. The inclusion criteria were as follows: (1) patients confirmed breast cancer by pathology; (2) HER-2 status determined by IHC and/or FISH; and (3) those who underwent multi-parametric MRI and US examinations. The exclusion criteria were as follows: (1) partial surgical resection or biopsy of the tumor before MRI or US imaging; (2) poor or incomplete image quality; (3) presence of other malignant tumors; and (4) incomplete clinical data.

Pathological Examination

Tumor tissues for pathological examination were obtained through surgery or biopsy. HER-2 status was determined via IHC and/or FISH assays.⁵ HER-2 positivity was defined as IHC 3+, while HER-2 negativity was defined as IHC 0 or IHC 1+. For IHC 2+ results, FISH was performed, and HER-2 was considered positive if gene amplification was observed; otherwise, it was deemed negative.

Clinical Characteristics

Clinical characteristics were extracted from patient medical records. These included body mass index, age, menopause status, family history of breast cancer, breast appearance changes, palpation findings, and serum biomarkers.

Ultrasound Examination

All patients underwent preoperative US examinations of the breast. The US machines used included GE LOGIQ E9 and VOLUSON E9 (General Electric, Boston, USA), HITACHI ARIEETTA 70 (HITACHI Ltd., Tokyo, Japan), and Mindray Resona 7 and DC-7 (Mindray, Shenzhen, China), equipped with linear array probes operating at a frequency of 9–12 MHz. Patients were positioned supine, with their arms raised above their heads to expose the breast fully. Two US physicians with 5 years of experience analyzed the US features, and any discrepancies were resolved by a third US physician with 10 years of experience. They were all blinded to the pathological results. The recorded features included multifocality, lesion location, ADLER blood flow grading, calcification, boundary, morphology, margin, internal echo, rear acoustic shadow, Breast Imaging-Reporting and Data System classification, and the maximum tumor diameter.

MRI Examination

MRI scans were performed using Siemens Verio 3.0T, Siemens Prisma 3.0T (Siemens Healthineers, Erlangen, Germany), and GE 750 3.0T imagers (General Electric, Chicago, USA). To minimize motion, patients were positioned prone. The imaging sequences included axial T1-weighted imaging (T1WI), axial turbo-spin-echo T2-weighted imaging (T2WI) with fat suppression (FS), axial single-shot diffusion-weighted imaging (DWI), and axial dynamic contrast-enhanced (DCE) imaging. Prior to DCE, a pre-contrast T1WI-FS scan was conducted, followed by an intravenous injection of gadolinium contrast (Gadopentetate Dimeglumine, Magnevist, 0.1 mmol/kg at 2 mL/s) and a 20 mL saline flush. The DCE protocol for Siemens Verio 3.0T included one pre-contrast and five post-contrast phases, each lasting 68 seconds. The Siemens Prisma 3.0T protocol comprised one pre-contrast and four post-contrast phases, each 77 seconds long. The GE 750 3.0T protocol consisted of one pre-contrast and five post-contrast phases, with each phase lasting 66 seconds. Detailed MRI sequence information is provided in Table S1.

MRI imaging features included tumor length, short axis, ADC value, time-intensity curve type (persistent, plateau, washout), presence of lymph nodes with an axillary long axis ≥1 cm, and tumor volume. Tumor volume was obtained from the volume of interest (VOI) in ITK-snap; other imaging features were based on MRI reports.

Clinical Model Development

Univariate logistic regression analysis (p<0.05) was conducted to identify significant variables. Variables with a variance inflation factor (VIF) greater than 5 were removed to address multicollinearity. Subsequently, multivariate logistic regression analysis (p<0.01) was performed to identify risk factors associated with HER-2 positivity, forming the clinical predictive model.

MRI Habitat Generation

All MR images underwent N4 bias correction for intensity normalization (http://stnava.github.io/ANTs/) and were resampled to a resolution of 1 × 1×1 mm. Image segmentation and registration were performed using ITK-Snap (Version 4.0.2, http://www.itksnap.org/). Two radiologists with 5 years of experience, who were blinded to patient results, used a semi-automatic method to delineate the VOIs on each slice of the second phase of axial dynamic contrast-enhanced (DCE2) images, with the final VOI determined by consensus. The VOIs were then applied to other sequences. T1 and the first phase of axial dynamic contrast-enhanced (DCE1) images were non-rigidly registered to DCE2, while T2 and DWI images were rigidly registered to DCE1, followed by manual adjustment for optimal alignment.

Habitats were generated by clustering the VOIs into sub-regions based on the voxel values from all MRI sequences of each sample. To better characterize the habitats, two subtraction sequences were created: SUB_1 by subtracting T1 from DCE1 and SUB_2 by subtracting DCE1 from DCE2. Voxel values from seven sequences (T1, T2, DWI, DCE1, DCE2, SUB_1, and SUB_2) were used for patient-level VOI clustering using the K-means algorithm. SUB_1 and SUB_2 were not included in subsequent analyses. The optimal number of habitats was determined by testing 3 to 10 clusters and selecting the k-value with the highest mean Calinski–Harabasz (CH) score.²¹

US Image Preprocessing and Segmentation

US images underwent z-score normalization to eliminate variations in grayscale values. Tumor regions of interest (ROIs) were delineated by US Physician 1 using the open-source software ITK-Snap on the cross-section of the grayscale US image showing the largest tumor dimension. For inter-observer consistency, US Physician 2 delineated ROIs for 30 randomly selected patients. The intra-class correlation coefficient (ICC) was used to quantify the consistency of radiomic features. Features showing an ICC ≥ 0.75 were considered robust.

Feature Extraction and Selection

Radiomic features were extracted using the PyRadiomics package (v3.6.2) from US or MRI images respectively, including first-order features, shape features, and texture features. Only radiomic features from US images with an ICC ≥ 0.75 were retained.

All features were standardized using Z-scores. Student’s t-test was applied to retain radiomic features. Using the Spearman rank correlation coefficient, a greedy recursive elimination strategy was applied to filter features. This approach retained only one feature from any pair with a correlation coefficient > 0.9, ensuring the retention of descriptive features. The least absolute shrinkage and selection operator (LASSO) regression model, tuned for regularization parameter λ, selected important non-zero coefficient features to form the optimal feature set for model construction.

Radiomics Model and Combined Model Development

We developed several predictive models using the extremely randomized trees (ExtraTrees) method: five models based on the combined features from the three habitat regions within each MR sequence (Model_T1__H1+H2+H3, Model_T2__H1+H2+H3, Model_DWI__H1+H2+H3, Model_DCE1__H1+H2+H3, Model_DCE2__H1+H2+H3), one US-based model, three models on the combined features from the same habitat region across all MR sequences (Model_H1__{multi-parametric}, Model_H2__{multi-parametric}, Model_H3__{multi-parametric}), and a habitat fusion model (Model_Habitat__whole) incorporating features from all habitat regions across all MR sequences.

Using Averaging-Based Stacking, we developed two joint prediction models. Model__H+US integrates the best habitat and US model by averaging their out-of-sample predicted probabilities into a meta-feature to train an ExtraTrees meta-learner for final predictions. Model__H+US+Cli additionally incorporates the clinical model, similarly averaging probabilities as the meta-feature to train its ExtraTrees meta-learner. Each model’s parameters were fine-tuned for optimal performance (Table S2).²²

The overall study workflow is summarized in Figure 2.

Figure 2 Workflow of the study. Abbreviations: DCE1, the first phase of dynamic contrast-enhanced imaging; DCE2, the second phase of dynamic contrast-enhanced imaging; DWI, diffusion-weighted imaging; T1WI, T1-weighted imaging; T2WI, T2-weighted imaging; VOI, volumes of interest; ROI, regions of interest; LASSO, least absolute shrinkage and selection operator; Model_H1+US, a model that combines the Habitat 1 based model and the ultrasound-based model; Model_H1+US+Cli, a model that combines the Habitat 1 based model, ultrasound based model and clinical based model.

Statistical Analysis

Statistical analyses were conducted using SPSS (version 26.0) and MedCalc (version 16.2). For missing variable values, discrete variables were filled using the median and continuous variables using the mean. A variable was discarded if it had 20% or more missing values. Categorical variables were compared using the chi-square test or Fisher’s exact test. Continuous variables were compared using the Student’s t-test or Mann–Whitney U-test. Model performance was evaluated using receiver operating characteristic (ROC) curves, accuracy, F1 score, sensitivity, and specificity. Area under the curve (AUC) differences between models were compared using DeLong’s test. Incremental differences between combined and single models were measured using IDI, and clinical net benefits were determined using decision curve analysis (DCA). SHAP analysis was performed to assess the interpretability of the best habitat model and US model. Differences with p<0.05 were considered statistically significant.

Results

Baseline Data of Patients

Our study included 182 breast cancer patients with a mean age of 49.57 ± 10.49 years. Of these, 88 (48.4%) were HER-2-positive and 94 (51.6%) were HER-2-negative. Patients were randomly divided into a training set (n=145) and a testing set (n=37) at an 8:2 ratio. The baseline characteristics of the patients are shown in Table S3. Univariate and multivariate analysis identified calcifications as an independent predictor of HER-2+ status (p<0.01, Table 1).

Table 1 Univariate and Multivariate Analyses of Patient Baseline Characteristics to Predict HER-2 Positivity in Breast Cancer

The clinical model achieved AUCs of 0.645 (95% CI=0.568–0.723) and 0.635 (95% CI=0.478–0.792) in the training and testing sets, respectively (Figure S1).

Habitat Radiomics Feature Selection Based on MRI

The highest CH score was observed when the tumors were clustered into three habitats across the entire cohort (Figure S2). From each MRI sequence’s habitat region, 107 features were extracted, resulting in a combined feature set of 321 features per MRI sequence. Features from the same habitat region across all five sequences were combined to form three feature sets, each containing 535 features, while all habitat-extracted features were amalgamated into a single set containing 1605 features. Further details on feature selection for the habitat 1 models are provided in Figure S3A, C, E.

Habitat Model Establishment and Evaluation Based on MRI

We developed nine radiomic models using the ExtraTrees algorithm. Table S4 presents the specific performance metrics of the nine models, while Figure S4 displays the ROC curves and DeLong’s test results. The Model_H1__{multi-parametric} achieved AUCs of 0.880 (95% CI=0.825–0.935) and 0.801 (95% CI=0.655–0.948) in the training and testing sets, respectively. As a result, the Model_H1__{multi-parametric} was ultimately selected as the best habitat model.

Radiomics Feature Selection and Model Establishment Based on US

There were 1561 radiomic features extracting from the US ROIs, with 1499 robust features retained after ICC testing. Details on feature selection for the US models are provided in Figure S3B, D, F. The US model’s AUC was 0.729 (95% CI=0.647–0.811) and 0.631 (95% CI=0.432–0.831) in the training and testing set, respectively (Figure S5).

Evaluation of the Best Habitat Model, US Model, Clinical Model, and Combined Model

Employing Averaging-Based Stacking, we developed Model__H1+US and Model__H1+US+Cli (Table 2).

Table 2 Performance of Combined Models in Training and Testing Sets

Detailed results of DeLong’s test, DCA curves, and calibration curves for Model_H1__{multi-parametric}, US model and the clinical model are shown in Figure S6. Figure 3 shows the feature importance, distribution and individual analysis through SHAP for both Model_H1__{multi-parametric} and US model.

Figure 3 Continued.

Figure 3 Continued.

Figure 3 SHAP-based model visualization. (A and B) Bar charts illustrating the weights of the most important features in Model_H1_multi-parametric (A) and ultrasound model (B). (C and D) The bees-warm plots depicting the impact of each feature on the prediction probability for Model_H1_multi-parametric (C) and ultrasound model (D). Red and blue points represent positive and negative impacts on the prediction outcome, respectively. (E and F) SHAP force plots demonstrating correct predictions by Model_H1_multi-parametric for representative HER-2 negative (E) and HER-2 positive (F) cases. (G and H) SHAP force plots demonstrating correct predictions by the ultrasound model for representative HER-2 negative (G) and HER-2 positive (H) cases. Abbreviations: SHAP, Shapley Additive explanation; h1.1, from the habitat 1 region of the second phase of dynamic contrast-enhanced images; h1.4, from the habitat 1 region of the T2-weighted imaging.

The ROC curves, DeLong test, and IDI results among Model_H1__{multi-parametric} and the two combined models are shown in Figure 4. Calibration curves showed good overall agreement across the three models (Figure 5A and B). DCA showed significant net benefits for all three models (Figure 5C and D). The radiomics prediction scores for each patient in Model__H1+US+Cli are presented in Figure S7.

Figure 4 Receiver operating characteristic curves, DeLong Test and Integrated Discrimination Improvement (IDI) Test Comparisons of Three Models on Training Set (A–C) and Test Set (D–F). Abbreviations: AUC, area under the curve; Model_H1_multi-parametric, radiomics model based on habitat 1; Model_H1+US, a model that combines the Habitat 1 based model and the ultrasound-based model, Model_H1+US+Cli, a model that combines the Habitat 1 based model, ultrasound based model and clinical based model.

Figure 5 Calibration Curves and Decision Curve Analysis for Model_H1_multi-parametric, Model_H1+US and Model_H1+US+Cli on Training Set (A and B) and Test Set (C and D).

Discussion

Our study developed a combined model by integrating the Model_H1__{multi-parametric}, US model, and clinical model using Averaging-Based Stacking, which demonstrated excellent performance in predicting HER-2 expression status in breast cancer (AUCs of 0.945 and 0.835 in the training and test sets, respectively).

Unlike previous studies that extracted radiomic features from the entire tumor region based on MRI, this study extracted features from habitat regions based on MRI. Habitat imaging, which emphasizes sub-regional radiomics analysis, may capture specific information related to gene or protein expression within different habitat regions.^23,24 In the meantime, Model_H1__{multi-parametric} exhibited high accuracy, F1 scores, sensitivity, and specificity, demonstrating better performance than Model_H2__{multi-parametric}, Model_H3__{multi-parametric}, or those based on multimodal multi-habitat clustering features or unimodal multi-habitat clustering features. This improved performance may be associated with the specific characteristics of HER-2-positive breast cancer in Habitat Subregion 1. The Model_Habitat__Whole, which incorporates features from multiple habitats, may suffer from diluted or masked information due to noise, reducing its predictive performance. Similar conclusions were drawn by Wang et al.²⁵

The primary methods for habitat generation are the k-means algorithm and the Otsu algorithm.^21,26,27 The Otsu algorithm divides an image into two categories: minimizing intra-class variance and maximizing inter-class variance. When habitat generation is based on n types of imaging data, the number of habitats generated is 2ⁿ. If we had used the Otsu algorithm, it would have yielded 32 habitat regions, making the results difficult to interpret and apply. Therefore, we chose the k-means algorithm. Notably, most previous studies have performed clustering at the patient group level,^21,26 ensuring the same physical meaning for each habitat category across different tumors. However, group-level clustering means that voxel values from all sample VOIs are divided into three clusters, and at the individual sample level, some clusters may be missing, resulting in incomplete radiomics features for that region. This issue has rarely been addressed in previous studies. Although missing radiomics features can be handled using methods for missing data, it may still decrease predictive performance and reduce model interpretability.²⁸ To avoid this, we used individual-level clustering in this study, ensuring each sample had three habitat regions. We addressed the potential issue of differing physical meanings for the same habitat region across samples by using the radiomics features of one sample as an anchor and rearranging the radiomics features of the three habitat regions for all samples according to their similarity to the anchor, ensuring consistency across samples.

Our Model__H1+US, which integrated the Model_H1__{multi-parametric} and the US radiomics model, achieved AUCs of 0.917 and 0.830 in the training and test sets, respectively, outperforming Model_H1__{multi-parametric}. This result indicates that adding US information enhances the model’s predictive ability. US imaging provides supplementary information distinct from MRI, making the model more comprehensive and accurate in capturing tumor characteristics. Multimodal radiomics has been applied to distinguish benign from malignant tumors,²⁹ predict treatment responses,³⁰ analyze prognoses,³¹ and predict biomarkers,³² all achieving higher predictive performance than single-modal approaches. To date, our study is the first to use a multimodal radiomics model to predict HER-2+ in breast cancer.

Multivariate analysis revealed that the presence of calcifications is an independent predictor of HER-2+ in breast cancer (p<0.01). Microcalcifications have been associated with higher malignancy levels in breast cancer.³³ By incorporating the clinical model, Model__H1+US+Cli demonstrated even better performance compared to Model__H1+US. In studies involving similar tasks, our Model__H1+US+Cli outperformed the convolutional neural network (CNN)-based deep learning radiomics model developed by Quan et al (AUCs: 0.917 [training set] and 0.810 [test set], evaluated on breast ultrasound videos from 445 patients).³⁴ This reconfirms the value of habitat imaging, which is the cornerstone of our innovative approach. Whereas the CNN model in Quan et al’s study analyzed the entire tumor region, habitat imaging provides deconstruction capability for the spatial ecosystem of tumors by identifying distinct habitats and characterizing their spatial relationships.

Interpretable machine learning techniques, such as SHAP, have effectively addressed the “black box” issue in model predictions.³⁵ Our SHAP analysis identified original_firstorder_Skewness as the top contributor to the Habitat 1 model, while wavelet_LLH_gldm_SmallDependenceLowGrayLevelEmphasis was key for the ultrasound model. A negative original_firstorder_Skewness value (quantifying asymmetry in the voxel value histogram) promoted a HER2-positive prediction. This indicates that within Habitat Region 1, a negatively skewed DWI distribution (characterized by a long left tail and most voxels at higher signals) suggests HER2-positivity. This distribution reflects high heterogeneity: predominant mid-high DWI signals correspond to a dense, proliferative core tumor, while the low-signal tail likely indicates micro-necrotic foci.³⁶ These foci may associate with HER2-positive cancer’s high proliferation, known to cause vascular compromise, elevated pressure, and hypoxia, leading to necrosis, as has been described in previous studies. Correspondingly, the wavelet_LLH_gldm_SmallDependenceLowGrayLevelEmphasis feature reflects the distribution characteristics of local texture details within low gray-level regions. The feature’s high value (positively contributing to the HER2-positive prediction) indicates scattered, isolated dark pixel areas on ultrasound. Interestingly, this imaging manifestation similarly points to potential necrotic regions. In addition, force plots for both positive and negative sample cases illustrated the contribution of each feature in individual patients, enabling the calculation of final Shapley values and supporting personalized prediction outcomes for each patient.

Our integrated model demonstrated high predictive accuracy in both training (AUC=0.951, Acc=0.892) and test sets (AUC=0.951, Acc=0.892), showing strong concordance with the IHC/FISH gold standard. Critically, unlike biopsies that sample only focal regions and fail to capture tumor-wide heterogeneity, our habitat imaging generates comprehensive heterogeneity maps, particularly highlighting subregions exhibiting higher association with HER2-positive expression. While not replacing gold-standard testing, the model serves as a valuable adjuvant tool aimed at enhancing efficiency, reducing costs, and alleviating clinical burden. For instance, upon successful prospective validation and clinical integration, the model could guide biopsy strategies: High-risk patients (probability >0.8) should receive prioritized biopsy with concurrent initial IHC/FISH testing to reduce repeat procedures, while low-risk patients (probability <0.3) could selectively forgo routine FISH testing though biopsy remains indicated. Moreover, during treatment monitoring, our approach enables real-time, repeatable assessment of HER2 status changes, facilitating timely therapeutic adjustments and reducing the need for serial biopsies – thereby alleviating patient trauma and financial strain. Furthermore, for equivocal pathology, extreme model probabilities (>0.9 or <0.1) provide strong radiological evidence; and in cases of insufficient/poor-quality samples, the model aids judgment on repeat biopsy necessity.

We acknowledge several limitations in our study. First, the retrospective, single-center design may limit the generalizability and external validity of our findings, which need to be validated through future prospective, multi-center, and larger sample size studies. Second, we used a semi-automatic method for tumor segmentation. In future research, we plan to explore fully automated segmentation methods. Third, the biological interpretation of radiomic habitat subregions remains theoretical without direct histopathological validation. We plan to validate these habitats through image-guided biopsies and correlate them with tumor microenvironment features using digital pathology.

Conclusion

In conclusion, this study successfully developed and validated an effective non-invasive tool for predicting HER-2 expression status in breast cancer. Specifically, the key innovation focuses on the utilization of advanced MRI habitat analysis to characterize tumor spatial heterogeneity. Building upon this foundation, we further developed an integrated multimodal model (Model__H1+US+Cli) that combines information derived from this novel MRI habitat approach with ultrasound and clinical data. Our findings demonstrate that this multimodal framework, leveraging the unique insights from MRI habitats, exhibited superior predictive performance compared to models based solely on MRI habitats, US, or clinical data. Therefore, exploring the combination of insights derived from advanced imaging phenotyping, such as MRI habitats, with other modalities and clinical data represents a promising direction for optimizing breast cancer radiomics models.