The Cancer Genome Atlas (TCGA) database (http://gdc.cancer.gov) was used to retrieve TCGA-breast invasive carcinoma (BRCA) gene expression and clinical data. After data pre-processing, 741 TCGA breast cancer samples were obtained for subsequent analysis. In addition, clinical information from the GSE7390 dataset was acquired from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/). After excluding specimens with an overall survival (OS) time of 0 days or missing survival data (due to data unavailability or inability to obtain), a total of 198 breast cancer samples were included in the analysis.

Kaplan-Meier survival curves were generated using the R package ‘survival’ (Version 3.2-13; http://bioconductor.org/packages/survivalr/) (11), with statistical comparisons performed by log-rank test. The Kaplan-Meier analysis was employed to evaluate the association between different HER2 sample groups (negative and positive) and OS. In addition, the association between HER2 and clinical data of breast cancer cases was evaluated by χ² or Fisher's exact tests.

The Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) (https://cibersort.stanford.edu/index.php) (12), Single-sample Gene Set Enrichment Analysis (ssGSEA) (http://www.bioconductor.org/packages/release/bioc/html/ssGSEA.html) (13) and xCELL (https://github.com/dviraran/xCell) (14) algorithms were adopted to assess the score of immune cells according to the mRNA expression matrix of breast cancer samples. The ‘ESTIMATE’ package (http://127.0.0.1:29606/library/estimate/html/estimateScore.html) (15) in R was employed to obtain stromal, immune and ESTIMATE fractions. The expression data of checkpoint genes, human leukocyte antigen family genes and chemokine genes were extracted from the breast cancer expression data, and the differences in expression levels between negative- and positive-HER2 sample groups were compared using the Wilcoxon rank-sum test.

Based on breast cancer sample mutation information, the mutation status of each gene in the samples was counted, genes were sorted in descending order by the number of mutations and then the TOP 20 genes with the highest mutations were selected for mutation display. Furthermore, the mutation frequency of the top 20 genes was analyzed using the ‘maftools’ package (Version 2.8.0) (https://bioconductor.org/packages/release/bioc/html/maftools.html) (16). In addition, the tumor mutation burden (TMB) of all cancer samples was determined, and the discrepancy in TMB values between the HER2-negative and -positive groups was analyzed.

The sensitivity of all cases to chemotherapy agents was assessed based on data obtained from the Genomics of Drug Sensitivity in Cancer database (https://www.cancerrxgene.org/), and the half maximal inhibitory concentration (IC₅₀) was quantified with the ‘pRRophetic’ package (https://github.com/paulgeeleher/pRRophetic) (17). Differences in the IC₅₀ values of 138 chemotherapeutics between HER2-negative and -positive groups were analyzed using the Wilcoxon rank-sum test.

The Tumor Immune Dysfunction and Exclusion (TIDE) database (http://tide.dfci.harvard.edu/) was used to analyze the response of patients to immune checkpoint treatment, which was represented as the TIDE score. In addition, the immunophenoscore (IPS) was used to determine the scores of the four different immune phenotypes (namely, inhibitory cells, effector cells, antigen presentation and checkpoints). Furthermore, the Gene Set Variation Analysis algorithm was employed to evaluate the tertiary lymphoid structure (TLS) scores of TLS feature genes [such as CC motif chemokine ligand (CCL)2, CCL3 and CCL4]. Finally, the differences in TIDE, IPS and TLS scores between HER2-negative and -positive groups were analyzed using the Wilcoxon rank-sum test.

GSEA (http://bioconductor.org/packages/release/bioc/html/GSVA.html) (18) was used to assess the significant enrichment pathways of hallmark gene sets (h.all.v7.4. symbols) between HER2-negative and -positive groups with P<0.05 and |(Normalized Enrichment Score|>1.

To identify HER2-specific genes, differential expression analysis was performed between HER2-positive and HER2-negative samples using the ‘limma’ package (Version 3.34.7) (https:// bioconductor.org/packages/release/bioc/html/limma.html) (19). This analysis generated gene-specific P-values, logFC (log₂ fold change), and other relevant metrics. Additionally, the limma package was employed to perform an empirical Bayes moderated t-test/F-test. The resulting P-values were sorted and adjusted using the Benjamini and Hochberg method to obtain false discovery rate (FDR)/q-values, thereby controlling the overall FDR. DEGs meeting |log2 fold-change|≥0.585 and adj. P<0.05 were used for subsequent analyses. In addition, enrichment analysis of DEGs was performed through the ‘clusterProfiler’ package (Version 4.0.5) (http://bioconductor.org/packages/release/bioc/html/clusterProfiler.html) (20) in R, with a threshold of count >2 and adj. P<0.05.

Univariate Cox regression analysis of DEGs was performed using the Kaplan-Meier ‘survival’ package (http://bioconductor.org/packages/survivalr/) (11) with a critical value of P<0.05 to identify genes associated with prognosis. The least absolute shrinkage and selection operator (LASSO) is a dimensionality-reduction method that has shown advantages over regression analysis in evaluating high-dimensional data. LASSO is an improvement of linear regression that achieves feature selection and complexity control by adding a penalty term. For parameter optimization, LASSO provides two criteria: i) ‘min’ selects the model achieving minimum cross-validation error within allowable variance; and ii) ‘1se’ (one standard error rule) chooses the most parsimonious model whose performance is statistically indistinguishable from the ‘min’ model, effectively balancing accuracy and simplicity. This algorithm can determine the optimal penalty coefficient based on the minimum likelihood deviation and 10-fold cross validation (21). Therefore, LASSO (21) was used to identify key genes. The risk score was constructed using stepwise Cox regression analysis with the survminer package in R (Version 0.4.9) (https://cran.rstudio.com/web/packages/survminer/index.html). The formula is as follows: RiskScore=h (t,X)=h₀ (t) × exp (β₁X₁ + β₂X₂ + ... + β_nX_n), where β indicates the regression coefficient; h₀(t), the benchmark risk function; and h (t,X) the risk role related to X (covariate) at time t. All specimens were classified into high- and low-risk groups according to the median RiskScores. Survival analysis was performed using the Kaplan-Meier curve method with the log-rank test for statistical comparison between groups. To evaluate the diagnostic performance of feature genes, receiver operating characteristic (ROC) curves and the area under the curve (AUC) were plotted across all datasets, including both training and validation sets.

The association between RiskScores and clinicopathological data (such as age and stage) was investigated, and univariate and multivariate Cox regression analyses were performed to screen for useful prognostic features, with a threshold of P<0.05. Next, the ‘rms’ package (Version 6.2-0) (https://cran.r-project.org/web/packages/rms/index.html) (22) was further utilized to construct a nomogram.

The ‘ssGSEA’ and ‘CIBERSORT’ algorithms were employed to evaluate the proportional distribution of immune cells between different risk groups. Spearman correlation analysis was performed using the ‘ggcor’ package (Version 0.9.8.1) (https://pan.baidu.com/s/1S6w93IjfO6sU8IHGvvCa1w) (23) to assess the association between RiskScore and immune cell infiltration.

The associations between the expression levels of key genes were determined with the ‘ggplot2’ package (Version 3.3.5) (https://github.com/tidyverse/ggplot2) (24), while the ‘RCircos’ package (Version 1.2.2) (https://github.com/hzhanghenry/RCircos) (25) was employed to draw the mutation gene plot of key genes. The Kaplan-Meier curve was used to assess the survival of patients. The log-rank test was employed to determine the statistical significance of differences between survival curves with a P-value threshold of <0.05.

Based on a single-cell RNA sequencing database focused on the TME accessed through the Tumor Immune Single-cell Hub (TISCH) database (http://tisch.comp-genomics.org), scRNA-seq data from GSE161529 (26) were analyzed. Detailed cell type annotation was performed at the single-cell level for subsequent analysis of specific gene expression across distinct cell populations (27).

Six HER2-negative and six HER2-positive breast cancer samples were collected from May 1, 2023, to December 31, 2023, at Weifang Hospital of Traditional Chinese Medicine (Weifang, China). The enrolled patients ranged in age from 32 to 69 years, and all samples were obtained through surgical resection. The Ethics Committee of Weifang Hospital of Traditional Chinese Medicine approved the present study (approval no. WF2023-428). All the participants provided written informed consent.

For RT-qPCR, total RNA was extracted from tissue samples using an RNeasy Kit (Qiagen, Inc. cat. no. 74104) according to the manufacturer's protocol. According to the manufacturer's instructions, the RNA samples were reverse transcribed into cDNA using an QuantiTect Reverse Transcription Kit (cat. no. 205311; Qiagen GmbH), and qPCR was performed using a CFX96 Bio-Rad system (Bio-Rad Laboratories, Inc.) with SYBR Green I fluorophore (cat. no. 1708882; Bio-Rad Laboratories, Inc.). The thermocycling conditions were as follows: Initial denaturation at 95°C for 3 min, followed by 40 cycles of 95°C for 10 sec and 60°C for 30 sec, with a final melting curve analysis from 65°C to 95°C in 0.5°C increments. The primer sequences, which were designed and synthesized by Sangon Biotech Co., Ltd., are listed in Table I. The internal reference gene for normalization was GAPDH. Each sample was analyzed in three technical replicates, and the relative levels of key genes were analyzed by the 2^−ΔΔCq method (28).

The HER2-positive breast cancer cell lines BT474 (cat. no. TCHu143) and SKBR-3 (cat. no. TCHu225) (29), and the HER2-negative breast cancer cell line MDA-MB-468 (cat. no. TCHu136) (30) were purchased from The Cell Bank of Type Culture Collection of The Chinese Academy of Sciences. All cells were cultured in RPMI 1640 medium (cat. no. 11875093; Gibco; Thermo Fisher Scientific, Inc.) containing 10% fetal bovine serum (FBS; cat. no. 10099141; Gibco; Thermo Fisher Scientific, Inc.) in an incubator with 5% CO₂ at 37°C. Cells at logarithmic proliferation stage were transfected. The ELOVL2 overexpression plasmid (OE-ELOVL2) and its negative control (NC-ELOVL2) were constructed using the pCMV6-AC-GFP vector backbone (cat. no. YC-13849RJ) by Shanghai Yaji Biotechnology Co., Ltd., and then transfected into cells using Lipofectamine^™ 2000 (cat. no. 11668019; Invitrogen; Thermo Fisher Scientific, Inc.) at a concentration of 2 µg per 1×10⁶ cells. After 6 h of transfection, the complete medium was replaced and the cells were further cultured for 48 h. The transfected cells were then collected for subsequent experiments.

BT474 and SKBR-3 cells were inoculated into 96-well plates with 1×10⁴ cells per well in a ~100 µl suspension, and cultured in a 5% CO₂ incubator at 37ºC in triplicate wells. After the cells reached 80–90% confluence (as observed by phase-contrast microscopy), CCK-8 reagent (cat. no. CK04; Dojindo Laboratories, Inc.) was added, and incubated at 37°C for 2 h. Next, the absorption value at 450 nm was measured using a microplate reader (ELx800; BioTek Instruments; Agilent Technologies, Inc.), and the cytotoxicity curves of the cells at different times were plotted. The experiment was repeated three times independently.

Cells were lysed using RIPA buffer (cat. no. P0013B; Beyotime Institute of Biotechnology) supplemented with 1 mM PMSF to extract the total proteins from each group. SDS-PAGE was performed with 30 µg protein loaded per lane on 10% polyacrylamide gels after determining the protein concentration in each group using a BCA kit (cat. no. 23225; Thermo Fisher Scientific, Inc.). After electrophoresis, the proteins were transferred onto a PVDF membrane, which was then blocked with 5% bovine serum albumin (MilliporeSigma) for 1 h at 37°C. Next, primary antibodies against PI3K (p110 alpha) (1:1,000; cat. no. 4255; Cell Signaling Technology, Inc.), phosphorylated (p)-PI3K recognizing p85 (Tyr458)/p55 (Tyr199) (E3U1H) (1:1,000; cat. no. 17366; Cell Signaling Technology, Inc.), AKT (pan) (11E7) (1:1,000; cat. no. 4685; Cell Signaling Technology, Inc.), p-AKT (Ser473) (D9E) XP^® (1:1,000; cat. no. 4060; Cell Signaling Technology, Inc.) and GAPDH (used as a loading control for normalization) (1:10,000; cat. no. 2118; Cell Signaling Technology, Inc.) were added and incubated at 4°C overnight. After that, the membranes were washed thoroughly to remove unbound primary antibodies, followed by incubation with the corresponding secondary antibody [Goat Anti-Rabbit IgG(H+L) (peroxidase/HRP conjugated); 1:10,000, cat. no. E-AB-1003; Wuhan Elabscience Biotechnology Co., Ltd.; and Goat Anti-Mouse IgG (H+L) (peroxidase/HRP conjugated); 1:10,000; cat. no. E-AB-1001; Wuhan Elabscience Biotechnology Co., Ltd.] at room temperature for 1 h. Next, a laser imager (Typhoon FLA 9500; GE Healthcare) was used for scanning with ECL Plus Western Blotting Substrate (Thermo Fisher Scientific, Inc.), and the gray values were analyzed after normalization to GAPDH as a loading control.

R software (SPSS 22.0) and GraphPad software (version 7.0; GraphPad; Dotmatics) were used for data processing, statistical analyses and plotting. Survival curves were constructed using the Kaplan-Meier method, and the log-rank test was employed to compare survival differences between groups, provided the proportional hazards assumption held. For time-to-event data with intermediate events, the two-stage survival analysis method from the Two Stage Hazard Rate Comparison package (Version 0.1-6) was applied when the proportional hazards assumption was violated. One-way analysis of variance was performed to assess overall differences, followed by pairwise comparisons between groups using Tukey's honestly significant difference test. In Table II, continuous variables with normal distribution are presented as the mean ± standard deviation, with intergroup comparisons assessed using unpaired Student's t-test; while skewed data are reported as median (interquartile range) and were compared via Mann-Whitney U test or Kruskal-Wallis H test followed by Dunn's multiple comparison post hoc test. Categorical variables are expressed as counts and percentages, with comparisons analyzed using χ² or Fisher's exact tests. Unpaired All statistical tests were two-sided and P<0.05 was considered to indicate a statistically significant difference.

As shown in Fig. 1A-D, no statistically significant differences were observed in prognosis between HER2-negative and -positive groups. Analysis of the association between clinical traits and HER2 suggested that HER2 was significantly associated with ER status, PR status, N stage, American Joint Committee on Cancer (AJCC) stage, mutation count and TMB (Table II).

Fig. 1E shows that the Stromal and ESTIMATE scores in HER2-positive group were higher than those in the HER2-negative group. Moreover, the ‘CIBERSORT’ algorithm revealed that the fraction of M2 macrophages and resting dendritic cells showed significant differences between the different HER2 groups (Fig. 1F). Additionally, the ‘ssGSEA’ and ‘xCELL’ algorithms revealed that there were 9 types (Activated CD4 T cell, T follicular helper cell, Gamma delta T cell, Type 17 T helper cell, CD56 bright natural killer cell, Myeloid derived suppressor cell, Plasmacytoid dendritic cell, Macrophage and Eosinophil) and 20 types of immune cells, respectively, which showed marked differences between different HER2 groups (Fig. 1G and H, respectively). The present study also identified 19 immune regulatory genes (Fig. 1I) and 11 chemokine genes (Fig. 1J) that exhibited remarkable differences between the HER2-negative and -positive groups. Fig. 1K demonstrates the differential expression of HLA family genes across distinct HER2 status groups.

Based on the TCGA-derived mutational profiles, gene-specific mutation frequencies across breast cancer cohorts were systematically analyzed. Fig. 2A and B illustrates the hierarchical distribution of the top 20 most frequently mutated genes, with Fig. 2A and B demonstrating distinct mutational landscapes between HER2-negative and HER2-positive subtypes, respectively. Mutation analysis showed that, compared with those of the HER2-negative group, the TMB values of the top 20 genes in the HER2-positive group were higher (Fig. 2C). Fig. 2D shows the top three chemotherapeutic drugs with significant differences between the two groups. In addition, no significant difference was observed in the TIDE score among the HER2-positive and -negative groups, whereas the interferon γ score in the HER2-positive group was higher than that in the HER2-negative group (Fig. 2E). The data showed remarkable discrepancies in the effector cell, immunosuppressive cells and immune checkpoint IPS scores among the different HER2 groups (Fig. 2F). In addition, the TLS scores showed remarkable differences among the HER2-positive and -negative groups (Fig. 2G).

GSEA identified six pathways of hallmark gene sets between different HER2 groups, namely, mTORC1 signaling, cholesterol homeostasis, glycolysis, heme metabolism, androgen response and KRAS signaling (Fig. 3A). A total of 251 DEGs between different HER2 groups were screened (Fig. 3B), and these DEGs were enriched in the GO terms of urogenital system development and RAGE receptor binding, and were involved in the following KEGG pathways: Estrogen signaling pathway, chemical carcinogenesis-receptor activation and PI3K-AKT signaling pathway (Fig. 3C).

A total of 39 prognostic related genes were screened from the 251 DEGs (Fig. 4A). Subsequently, LASSO regression analysis was performed using univariate Cox regression with a P-value threshold of <0.05, which narrowed down the candidates to 25 genes (Fig. 4B). Moreover, eight key genes were identified after stepwise Cox regression analysis (Fig. 4C): ETFA, RAPGEFL1, KRT7, CD24, PRR15L, ALOX15B, ELOVL2 and CXCL9. A risk model was constructed based on the aforementioned key genes. As shown in Fig. 4D and E, compared with the low-risk group, patients in the high-risk group had a notably poorer prognosis (P<0.01), and AUC of the ROC curve was 0.735, 0.784 and 0.752 for the 1-, 3- and 5-year survival, respectively, in TCGA dataset (Fig. 4D), vs. 0.787, 0.796 and 0.771, respectively, in the GSE7390 dataset (Fig. 4E).

Fig. 5A shows marked differences in Risk score between different HER2 groups (positive vs. negative), M stage (M0 vs. M1), T stage (T1 vs. T4, T2 vs. T4, T3 vs. T4) and AJCC stage (I vs. II, I vs. IV, II vs. IV, III vs. IV). Fig. 5B demonstrates that the high-risk group exhibits a significant tendency toward HER2-negative cases, whereas the low-risk group shows a higher proportion of HER2-positive cases. The HER2 status indirectly influences survival outcomes by affecting AJCC staging. HER2-negative cases are more likely to be classified as early-stage (I/II), associated with higher survival rates, while HER2-positive tumors tend to progress to advanced stages (III/IV), linked to elevated mortality risks. AJCC stage is strongly associated with survival outcomes: stages I, II, and III primarily exhibit NO (survival), whereas stage IV predominantly presents YES (death). Univariate Cox regression analysis was performed on clinical factors and risk score of the samples. Factors with P<0.05 were selected for subsequent multivariate Cox regression analysis to identify significant independent prognostic factors. Age, risk score and M stage were considered as independent prognostic factors (Fig. 5C and D). These factors were then used to build a nomogram, demonstrating that each variable independently contributed to survival probability prediction. Increased age corresponded with lower survival rates, advanced M stage significantly reduced survival probability, and higher risk scores were associated with poorer prognoses. Total points provided an intuitive prediction of outcomes, such as the 5-year survival rate. Model validation revealed high discriminatory ability (C-index=0.748) (Fig. 5E). The survival rate model predicted by the nomogram shows overall consistency with the actual survival rates for 1-, 3- and 5-year predictions. Specifically, the 1-year survival rate predictions are closely clustered around the ideal dashed line, indicating better consistency. The 3-year survival rate data points are more dispersed but still remain around the dashed line, demonstrating good consistency. However, the 5-year survival rate data points exhibit increased dispersion, with some deviating from the dashed line, leading to a decline in consistency. Overall, the consistency decreases as the predicted time span increases (Fig. 5F). Fig. 5G shows the marked association of these factors with patient outcomes. The ROC for the nomogram suggested that the AUCs at 1, 3 and 5 years were 0.849, 0.779 and 0.790, respectively (Fig. 5H).

A significant difference was observed in the proportions of 7 and 24 immune cell types between the high- and low-risk groups according to the results of ‘CIBERSORT’ and ‘ssGSEA’ analyses, respectively (Fig. 6A and B, respectively). The associations of 22 and 28 immune cell types with eight key genes are shown in Fig. 6C and D, respectively.

Spearman's correlation test was utilized to analyze the association between the expression levels of the aforementioned eight key genes and the prognosis of patients with breast cancer was analyzed revealed a positive trend in general (Fig. 7A). Based on the circular graph, the chromosomal localization distribution and copy number variation (CNV) status of 8 model genes were visually mapped, clearly presenting the positions of each gene on the chromosome and their CNV alteration characteristics (Fig. 7B). The gain and loss trends are displayed in Fig. 7C. Fig. 7D shows heatmaps of the aforementioned eight key genes. Survival analysis showed that the prognosis of patients with high expression of CD24, ETFA and PRR15L was worse than that of patients with low expression, whereas RAPGEFL1, KRT7, ALOX15B, ELOVL2 and CXCL9 showed the opposite trend (all P<0.05; Fig. 7E).

The GSE161529 single-cell dataset was used to explore the expression of the aforementioned eight crucial genes in the immune microenvironment, and the results revealed the involvement of 11 cell types (Fig. 8A). The distribution and percentage of these cell types are displayed in Fig. 8B. The largest proportions corresponded to malignant and epithelial cells. Fig. 8C shows the distribution of the RiskScore in the 11 cell types. In addition, immune microenvironment analysis revealed that the aforementioned eight key genes were expressed in various immune cells (Fig. 8D).

RT-qPCR was performed to validate the expression of eight crucial genes (ETFA, RAPGEFL1, KRT7, CD24, PRR15L, ALOX15B, ELOVL2 and CXCL9). As shown in Fig. 9, the assays confirmed that ETFA, RAPGEFL1, KRT7, CD24 and PRR15L were significantly upregulated in the HER2-positive groups, whereas ELOVL2, CXCL9 and ALOX15B were significantly downregulated, compared with the HER2-negative group. ELOVL2 is a novel tumor suppressor with low expression in breast cancer (31). In prostate cancer, high expression of ELOVL2 suggests improved prognosis, and small hairpin RNA targeting ELOVL2 promotes cell proliferation, colony formation, migration and invasion, as well as the growth of subcutaneous xenografts by activating the PI3K/AKT signaling pathway (32). Therefore, ELOVL2 was selected in the present study for follow-up experiments.

The RT-qPCR results showed that the expression of ELOVL2 was significantly lower in BT474 and SKBR-3 cells than in MDA-MB-468 cells (P<0.001; Fig. 10A). In addition, RT-qPCR revealed that OE-ELOVL2 significantly increased the expression of ELOVL2 in BT474 and SKBR-3 cells, suggesting successful transfection (P<0.001; Fig. 10B). Furthermore, the CCK-8 proliferation assay showed that OE-ELOVL2 significantly inhibited the proliferation of BT474 and SKBR-3 cells (P<0.01; Fig. 10C). Western blot analysis showed no significant difference in the expression levels of PI3K or AKT in the OE-ELOVL2 group compared to the negative control group, whereas overexpression of ELOVL2 could decrease the levels of p-PI3K and p-AKT proteins. In addition, compared with the corresponding control group, the p-AKT/AKT and p-PI3K/PI3K ratios in the OE-ELOVL2 groups of both BT474 and SKBR-3 cell lines were significantly decreased (P<0.001; Fig. 10D). In conclusion, ELOVL2 overexpression inhibits HER2-positive breast cancer cell proliferation by inhibiting the PI3K-AKT pathway.