Bioinformatics-based multi-omics and machine learning analysis identifies stemness-associated molecular subtypes and a prognostic index in breast cancer

Introduction

Breast cancer (BC) is the most widespread cancer in women worldwide and the cause of the largest number of female cancer-related deaths. By 2022, there will be approximately 2.3 million new cases of BC worldwide, representing 23.8% of new female cancer cases (1). Simultaneously, about 666,000 women will die from BC, accounting for 15.4% of female cancer deaths (1). The impact of BC on female health has become a global public health issue. Currently, the prognosis in BC patients relies mainly on pathological features, such as tumor stage and molecular subtypes, for prediction, but they often have large errors (2). For example, low-risk patients may receive unnecessary treatment, whereas high-risk patients who are not properly treated often face the risk of cancer recurrence or metastasis. Existing treatments for BC include surgery, chemotherapy, radiotherapy, and endocrine therapy. However, the therapeutic effect and prognosis of BC patients, especially metastatic patients, are unsatisfactory. Treatment resistance, as well as tumor recurrence and metastasis, remain major challenges in the treatment of BC (3). BC is a highly heterogeneous disease. While conventional molecular subtypes provide a foundation for prognostic evaluation and therapeutic decisions, they remain inadequate in accounting for the diverse treatment responses and clinical outcomes observed among patients within a given subtype (2). These limitations are particularly evident in predicting tumor recurrence, metastatic progression, and therapeutic resistance. In light of these challenges, efforts are underway to develop more refined classification frameworks. For example, Lu et al. developed a prognostic model based on sulfur metabolism-related genes, revealing the crucial role of metabolic reprogramming in BC heterogeneity and its correlation with immune microenvironment features (4). Similarly, Fang et al. established a molecular subtype system centered on coenzyme Q10 metabolism, which not only serves as a prognostic indicator but also discriminates responders to immunotherapy, thereby providing valuable insights for personalized treatment strategies (5). Together, these studies underscore the potential of biology-driven classification systems to transcend the limitations of conventional subtyping and highlight the importance of molecularly tailored approaches in the pursuit of precision medicine for BC patients. Therefore, the search for innovative, reliable, and effective novel biomarkers to distinguish different subgroups of BC patients is important for individualized BC diagnosis and treatment.

Many tumors, including BC, contain a small percentage of cells, called cancer stem cells (CSCs), with extremely high tumorigenicity and extremely low differentiation. Several investigations have shown that breast cancer stem cells (BCSCs) induce tumorigenesis, promoting tumor progression, cancer recurrence, metastasis, and resistance to conventional treatments (6-8). This leads to poor curative effects and prognostic outcomes for BC patients. Notably, a variety of molecular markers associated with BCSCs have been identified, for example, CD44, CXCR4, ALDH1, ABCG2, and CD24 (9). Among them, CD44, ALDH1 and CD24 can forecast the prognosis in triple-negative patients (10,11). Moreover, CD24 can promote immune escape in triple-negative patients by regulating programmed cell death ligand 1 (PD-L1) expression, which has the potential to be exploited as an immunotherapeutic target for triple-negative patients (12). Furthermore, ALDH1A1 enzymatic activity induces myeloid-derived suppressor cell expansion and immunosuppression, which in turn promotes BC progression (13). Therefore, the search for BC molecular markers from the new viewpoint of BCSCs will provide a new research direction.

The tumor microenvironment (TME) is complex and diverse. There is growing evidence that BCSCs interact with the TME and regulate multiple signaling pathways, affecting tumor evolution and immune evasion. The cytokine CXCL1 can maintain self-regrowth of BCSCs and promote tumor development and immune escape (14). Tumor-associated macrophages (TAMs) have an established place in the BC microenvironment: on the one hand, TAMs inhibit BC development by secreting a variety of inflammatory factors; on the other hand, IL6 secreted by TAMs activates the JAK/STAT3 signaling pathway and stimulates the self-renewal of BCSCs, thus promoting the development of BC (15). Tumor-associated fibroblasts secrete SDF-1, which binds to CXCR4 on the surface of BCSCs, activates the WNT/β-catenin signaling pathway, and promotes the generation of BCSCs to exert a durable tumorigenic and metastatic capability in BC (16). Triple-negative subtypes expressing high levels of PD-L1 also have robust stemness properties (17). These results suggest that a detailed investigation of the interaction between BCSCs and TME may reveal novel therapeutic possibilities. However, previous research on BCSCs has been constrained by two major limitations: firstly, it has largely been confined to exploring the function of single or a few markers, lacking a systematic approach to construct a molecular subtype system based on BCSCs heterogeneity as a classification criterion; secondly, it has failed to adequately elucidate the association between different BCSCs subtypes and clinical outcomes (such as patient prognosis and treatment resistance patterns), which has significantly limited its clinical translational value. Consequently, the establishment of a novel molecular subtype classification system based on BCSCs heterogeneity represents a significant clinical imperative.

Since CSCs are key factors that regulate cancer development and treatment resistance, existing studies have constructed a novel characterization of CSCs by machine learning using one-class logistic regression (OCLR). On this basis, a novel mRNA stemness index (mRNAsi), an indicator of tumor cell stemness, has been proposed to quantitatively reflect the differentiation patterns in tumor cells at different stages, as well as the distribution characteristics in tumor tissues (18,19). Therefore, this study utilised mRNAsi derived from the OCLR algorithm based on transcriptomic data to evaluate stem cell characteristics in BC patients. Bioinformatics analysis was employed to identify mRNAsi-associated genes, establishing a novel molecular subtyping framework grounded in tumour stemness heterogeneity. We systematically evaluated the characteristics of different stemness molecular subtypes across pathways, tumour mutations, immune infiltration, and immunotherapy response prediction. Ultimately, machine learning was employed to construct mRNAsi-associated prognostic models. Establishing novel molecular subtyping based on BCSCs heterogeneity holds promise for more precise prognostic risk assessment and personalised treatment strategy formulation, offering new avenues for overcoming tumour resistance and improving patient survival. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-837/rc).

Methods

Data acquisition and pre-treatment

The gene expression profiling and clinical information from BC patients (n=973) were acquired from The Cancer Genome Atlas (TCGA-BRCA). Among the clinical variables, they include the age at diagnosis, tumor size (T), lymph node status (N), distant metastasis (M), pathological tumor stage, estrogen receptor (ER), progesterone receptor (PR), human epidermal growth factor receptor 2 (HER2), and the molecular classification of BC. The primary prognostic outcome endpoint for analysis was overall survival (OS), defined as the time from diagnosis to death from any cause. The expression matrix and survival data for GSE20685 were taken from the Gene Expression Omnibus (GEO) database, a dataset with 327 BC patients. In addition, transcriptomic and survival information for 1,403 BC patients was extracted from the METABRIC database. Somatic mutation data from BC individuals were loaded from TCGA and displayed using the R package “maftools”. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Overall study design

The analytical workflow of this study is depicted in Figure S1. In summary, the TCGA-BRCA transcriptomic dataset was utilised to initially quantify tumour stemness by calculating the mRNAsi employing the OCLR algorithm. In order to identify genes that are central to the stemness phenotype, weighted gene co-expression network analysis (WGCNA) was performed in order to identify modules that are highly correlated with mRNAsi. Subsequently, genes from significant modules were subjected to univariate Cox analysis in order to filter for those with prognostic value. These were then used to establish a novel stemness-based molecular classification of BC patients. This classification was comprehensively delineated in terms of its clinicopathological features, tumour immune microenvironment, and tumor mutation burden (TMB). In order to translate these findings into a clinically applicable tool, a prognostic signature was constructed from key genes using least absolute shrinkage and selection operator (LASSO) regression, and its robustness was assessed in two independent validation cohorts (GSE20685 and METABRIC). Subsequently, an independent prognostic analysis was conducted to ascertain whether this signature could serve as an independent factor for predicting the prognosis of patients diagnosed with BC. Ultimately, by integrating the risk model with key clinical variables, a quantitative tool was developed to quantify individualised risk in BC patients, and its predictive accuracy was evaluated.

Calculating the mRNA stemness index (mRNAsi)

Mean RNA-seq data for pluripotent stem cells (PSCs) were acquired from the Progenitor Cell Biology Consortium (PCBC) database using the “synapser” R package. Stemness features were identified by the OCLR algorithm and tested by leave-one-out cross-validation (18). Then, the gene expression matrix of the transcriptome was subjected to Spearman’s correlation analysis with the stemness profile, and a linear transformation was used to obtain the mRNAsi for each sample. The mRNAsi ranges from 0 to 1. A higher mRNAsi indicates a more advanced degree of tumor dedifferentiation. The mRNAsi was connected with the degree of tumor dedifferentiation and the biological process of CSCs.

WGCNA

To identify mRNAsi-related, co-expressed gene modules in TCGA-BRCA, the “WGCNA” R package was utilized (20). For network construction, the network topology is analyzed using different soft-threshold powers (β). The minimum value of β that matches the distribution of the scale-free network is chosen to compute the adjacency matrix. Next, to detect the genetic connectivity of the network, the adjacency matrix is converted into a topological overlap matrix and a correlation dissimilarity matrix. Similar co-expressed genes were clustered into different gene modules (minModuleSize =70). Target modules associated with mRNAsi were selected by calculating the correlation between each module and the stemness index. Subsequently, genes within the target modules were submitted to Reactome analysis by applying the ‘clusterProfiler’ R package (21,22).

Identification of stemness molecular subgroups

After selection of the target modules that were specifically correlated with mRNAsi, we carried out univariate Cox regression analyses for genes in the target modules to search for those that were robustly related to overall survival in BC individuals (P<0.001). Consensus clustering of cohorts, in TCGA-BRCA, was performed using mRNAsi-associated genes linked to prognosis in BC patients via the “ConsensusClusterPlus” R package (23).

Differential expression analysis and functional enrichment analysis

The differential expression of stemness subgroups was analyzed using the “limma” R package [|log₂ fold change (FC)| >1, false discovery rate (FDR) <0.05] (24). At the same time, differential genes were analyzed via Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment using the “clusterProfiler” R package to identify functional distinctions between the two subtypes (21,25,26). In addition, hallmark gene sets were fetched from the Molecular Signatures Database (MSigDB) and gene set enrichment analysis (GSEA) was performed between stemness molecular subtypes (27). A normalized enrichment score (NES) was created to assess differences in cancer-associated classical pathways between molecular subtypes.

Analyzing immune infiltration and predicting reaction to immunotherapy

The “estimate” R package was utilized to assess differences in stromal, immune, ESTIMATE score and tumor purity among stemness molecular subtypes (28). ESTIMATE is a tool that uses gene expression profiles to derive stromal and immune cell abundance in tumor tissue. The MCPcounter algorithm can use transcriptome data to quantify the absolute abundance of stemness molecular subtypes for 8 immune cells (CD3⁺ T cells, CD8⁺ T cells, cytotoxic lymphocytes, NK cells, B lymphocytes, monocytes, myeloid dendritic cells, and neutrophils) and 2 stromal cells (fibroblasts and endothelial cells) (29). In addition, based on the principles of linear support vector regression, CIBERSORT is a framework for the deconvolution of known gene-expressed profiles characterizing immune cells (30). Immune cell infiltration between stemness molecular subtypes was assessed by the “CIBERSORT” R script with the gene expression signature set LM22 of 22 immune cell subtypes.

We obtained the information within 20 immune checkpoints from the HisgAtlas database and analyzed the variation in the levels of stemness molecular subtypes across the 20 immune checkpoints (31). At the same time, we predicted the immunotherapy response of each BC patient by the TIDE database (32). TIDE is a computational framework that combines the expression profiles of genes characterizing T cell dysfunction and exclusion in tumors. Increased TIDE scores indicate poorer response rates to immunotherapy.

Creation and testing of an mRNAsi-based risk-prognostic model

The cohort in TCGA-BRCA was taken as the training group, and LASSO regression analysis was applied to construct an mRNAsi-related risk model. GSE20685 and METABRIC were selected as validation sets. After computing the risk score for BC individuals, the R package “survminer” was utilized to separate the cohort into high-risk and low-risk groups under the ideal threshold. Survival analysis was carried out via the “Survival” R package and the “timeROC” R package, and ROC curves were drawn over time (33).

Nomogram development

To determine whether the mRNAsi-related risk score could be used as an independent factor to forecast the prognosis for BC patients, univariate and multivariate COX regression algorithms were performed. To avoid collinearity between composite staging and its components, our multivariate Cox model incorporated the risk score together with the individual T, N, and M categories rather than the composite pathological stage. This approach ensures a statistically robust evaluation of the risk score’s independent prognostic value by controlling for core cancer progression drivers while minimizing variance inflation.

Next, a nomogram was constructed by integrating the stemness-based risk score with established clinical prognostic factors (including age, T, N and M) using the R packages ‘rms’, ‘regplot’, ‘foreign’ and ’survival’. This integration was motivated by the imperative to enhance the clinical translatability of our novel prognostic signature: while the stemness-based risk score provides novel biological insight, clinical parameters such as age and T/N/M categories are universally available, well-validated, and immediately interpretable to clinicians. This hybrid model bridges the gap between discovery science and routine clinical practice, thereby providing a more holistic and actionable risk assessment for potential clinical use. The discriminative ability of the nomogram was assessed using Harrell’s concordance index (C-index). To correct for potential overfitting and obtain a robust estimate of model performance, bootstrapping with 1,000 resamples was performed using the “rms” package in R. Subsequently, calibration curves were drawn to assess the agreement between the nomogram-predicted and the actual observed probabilities of 3- and 5-year survival.

Drug sensitivity analysis

To determine whether risk profiles can guide clinical medication use, the “oncoPredict” R package has been implemented to calculate response to common therapeutic drugs for each BC patient using the gene expression matrix of cell lines from the GDSC database and the corresponding drug susceptibility data as a training set, and to compare between two risk groups (34).

Statistical analysis

The R programming language was chosen for the statistical analysis of all data. Student’s t-test was employed when comparisons between two groups were consistent with independence, normality and homogeneity of variance, and the Wilcoxon rank sum test when not. Similarly, multiple group comparisons were chosen using the analysis of variance or the Kruskal-Wallis test. All statistical tests conducted were two-sided. If no P value was assigned to indicate a statistical difference, the default value of P<0.05 was taken as statistically significant.

Results

mRNAsi-related clinical and immunological features

Firstly, an mRNAsi was computed for everyone in the TCGA-BRCA cohort. Figure 1A shows that BC tissue samples had significantly higher mRNAsi than paired normal samples. The mRNAsi was not statistically different between age groups (Figure 1A). On this basis, TCGA-BRCA patients were sorted by mRNAsi, and the correlation with clinicopathological and immunological characteristics of BC was investigated. The heatmap reflects changes in mRNAsi for BC patients who had altered clinical and immunological characteristics (Figure 1B). The beeswarm plot further shows a gradual increase in mRNAsi with increasing tumor size (Figure 1C). Moreover, the mRNAsi was considerably greater in stage II/III individuals than in stage I individuals (Figure 1C). However, the association of mRNAsi with lymph node metastasis and distant metastasis was not significant. (Figure 1C). In addition, mRNAsi was significantly elevated in ER-negative, PR-negative, and HER2-positive BC patients (Figure 1D). Triple-negative patients had the highest mRNAsi, while luminal subtype patients had the lowest mRNAsi (Figure 1D). Spearman’s correlation showed that mRNAsi had an inverse relationship with stromal, immune, and ESTIMATE scores, and a positive relationship with tumor purity (Figure 1E). These results suggest that clinical and immune characteristics were related to mRNAsi changes.

Figure 1 Relationship between the mRNAsi and clinical and immunological features of BC. (A) mRNAsi levels between BC tissue and paired normal tissue and between age groups. (B) Heatmap of clinical and immunological features of BC patients sorted by mRNAsi. (C) mRNAsi levels in relation to tumor size, lymph node status, distant metastasis and AJCC stage. (D) mRNAsi levels between ER, PR and HER2 molecular states and BC molecular subtypes. (E) Spearman correlation of mRNAsi with immune characteristics (stromal, immune, ESTIMATE score and tumor purity). AJCC, American Joint Committee on Cancer; BC, breast cancer; ER, estrogen receptor; HER2, human epidermal growth factor receptor 2; M, distant metastasis; mRNAsi, mRNA stemness index; N, lymph node status; PR, progesterone receptor; T, tumor size; TN, triple-negative.

Identification of mRNAsi-related gene modules

Since the mRNAsi is different in different BC patients, it is believed that there are genes that can regulate the stemness in tumor cells. WGCNA was run on the TCGA-BRCA dataset to detect gene modules associated with mRNAsi. The least β value found to be consistent with a scale-free network distribution was chosen to be 5 (Figure 2A). Then, after dynamic tree pruning and average hierarchical clustering, 11 gene modules were identified from the gene clustering tree (Figure 2B). Subsequently, exploring the link between gene modules and mRNAsi, the green module (r=0.71, P<0.001), brown module (r=−0.71, P<0.001), and magenta module (r=−0.73, P<0.001) exhibited a high correlation with mRNAsi (Figure 2C). Among these modules, the brown and magenta modules showed a strong negative correlation with mRNAsi, while the green modules showed a strong positive correlation with mRNAsi (Figure 2C).

Figure 2 Gene module identification associated with mRNAsi. (A) Each soft threshold, scale-free topological model fit index (left) and mean connectivity (right). (B) Gene clustering tree. (C) Genetic module-feature correlation mapping. *, P<0.05; **, P<0.01; ***, P<0.001; ns, no significance. mRNAsi, mRNA stemness index.

This will be used to further validate the efficiency of the three modules. Scatter plots showed that there was a good correlation between module membership and the gene significance of these modules (Figure 3A). The target modules for subsequent study were therefore the green, brown, and magenta modules. Pathway enrichment analysis found that the genes within the magenta module were mainly focused on the pathways of G protein-coupled receptor (GPCR)-ligand binding, nuclear receptor signaling, extracellular matrix (ECM) organization and NTRK signaling (Figure 3B). The genes in the brown module focus on ECM degradation, ECM proteoglycans and collagen formation in the context of the TME (Figure 3B). In contrast, the genes in the green module are focused on the cell cycle pathway associated with mitosis (Figure 3B). The genes in the targeted modules are strongly linked with tumor development.

Figure 3 Pathway analysis of target gene modules. (A) Relationship between MM and GS for green, brown, and magenta modules. (B) Reactome pathway for gene enrichment in green, brown, and magenta modules. GS, gene significance; MM, module membership; mRNAsi, mRNA stemness indices.

Identification of stemness subtypes using mRNAsi-related genes

All genes within the target module were further analyzed using a univariate Cox regression algorithm to investigate mRNAsi-related genes. Eleven mRNAsi-related genes were obtained that were linked to overall survival in BC patients (Figure 4A). Notably, only BEND5, TNN, and PDLIM4 were protective factors, while the other 8 genes were risk factors (Figure 4A). The above 11 genes were used as target genes for subsequent studies. Figure 4B shows that the majority of these target genes had a positive correlation. Considering that BC is a highly heterogeneous tumor, the expression of mRNAsi-related genes might be different in different BC patients. To this end, the current work was undertaken to provide a molecular typing system, in terms of stem cell properties, using consensus clustering techniques based on target genes. In the consensus clustering process, the clustering variable (k) was used to divide the samples into 2 to 10 classes. From the CDF and CDF delta area curves, the area below the curve increased at least at k=3 (Figure 4C). Since the CDF delta area curve reflects the relative change in area under the curve from k to k-1, 2 is the most adequate value of k (Figure 4C). Figure 4D shows that the TCGA-BRCA cohort was divided into 2 subgroups defined as C1 (476 cases, 48.92%) and C2 (497 cases, 51.08%). The mRNAsi in C2 patients was higher than in C1 patients, suggesting that C2 patients had higher tumor stemness (Figure 4E). The variations in clinicopathological characteristics between the two stemness subgroups were then compared. C2 patients had poorer outcomes, more patients with intermediate to advanced stages, and a major concentration of HER2-overexpressing and triple-negative patients (Figure 4F, Table S1). In contrast, early-stage patients and luminal subtype patients were mainly concentrated in C1 (Table S1). The above results suggest that stemness subtypes can be used as a marker for prognostic assessment and deserve further in-depth study.

Figure 4 Determination of the molecular stemness subtypes of BC in terms of mRNAsi-related genes. (A) HR and P values for 11 genes with significance in BC prognosis. (B) Heatmap reflects the 11 genes correlated with each other. (C) The area of relative variation under the CDF curve for the k range 2–10 (left) and CDF curves in the k range 2–10 (right). (D) Clustering heat map for k=2. (E) mRNAsi levels between the two stemness subgroups. C1 and C2 on the x-axis represent the two consensus stemness subgroups identified in this study. (F) Comparison of Kaplan-Meier curves for survival in the two stemness subgroups. BC, breast cancer; CDF, cumulative distribution function; HR, hazards ratio; mRNAsi, mRNA stemness indices.

Distinction in regulatory pathways and patterns of somatic mutations among stemness subtypes

Because of the heterogeneity of stemness characteristics, the molecular pathways and potential mechanisms of stemness subtypes were further explored. Firstly, the differential expression of the two subgroups was analyzed, and provided 352 differentially expressed genes comprised of 111 upregulated and 241 downregulated genes (Figure 5A). The differential genes were mainly involved in cell cycle-related biological activities and were associated with collagen-containing ECM, chromosomes, and other cellular components (Figure 5B). In addition, its major molecular functions include microtubule and associated protein binding, ECM composition, glycosaminoglycan binding, peptidase regulators and their inhibitor activity, and endopeptidase regulator activity (Figure 5B). Figure 5C illustrates that the differential genes were markedly enriched for the cell cycle pathway. Further exploration of the classical cancer pathway differences between stemness subtypes showed that both C1 and C2 were enriched for multiple cancer-related signaling pathways and metabolism-related pathways (Figure 5D). For example, C2 was enriched in the mTORC1, the PI3K/AKT/MTOR and glycolytic pathways (Figure 5D). Additionally, pathways that regulate the cell cycle and proliferation, such as the E2F-TARGETS, and G2M checkpoints, were markedly activated in C2 (Figure 5D).

Figure 5 Regulatory pathway differences between the two stemness subtypes. (A) Volcanic map of differential genes between stemness subtypes. Red dots represent upregulated genes in the C2 subgroup compared to the C1 subgroup, blue dots represent downregulated genes. (B) GO terms that were substantially enriched in the DEGs between C1 and C2. (C) KEGG pathway analysis with highly significant enrichment for DEGs between C1 and C2. (D) Ridge diagram shows the signaling pathway activity of C1 and C2. DEGs, differentially expressed genes; FC, fold change; FDR, false discovery rate; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Accumulation of somatic mutations results in the growth of cancerous tumors. Thus, the landscape of somatic mutations between stemness subtypes was analyzed. There is a relatively high frequency of somatic mutations in C2, characterized by high mutation frequencies in TP53, PIK3CA, and TTN. Of these, PIK3CA tended to be the most mutated gene in both typologies, and it was the most frequently mutated gene in C1 (Figure 6A,6B). TP53 was a common mutation in C2 (Figure 6B). The expression level of mRNAsi was found to be proportional to TMB, and TMB was notably greater in the C2 subgroup. (Figure 6C,6D). The results suggest that there may be potential for variation in immune infiltration and response to immunotherapy in the stemness subtype.

Figure 6 Differences in somatic mutations between the two stemness subgroups. (A) Waterfall diagrams show the top 10 most mutated genes in C1 subgroup. (B) Waterfall diagrams show the top 10 most mutated genes in C2 subgroup. (C) Relationship between TMB and mRNAsi. (D) Variation in TMB between the two stemness subtypes. C1 and C2 on the x-axis represent the two consensus stemness subgroups identified in this study. mRNAsi, mRNA stemness index; TMB, tumor mutational burden.

Stemness subtypes have different TME status and immunotherapy response

Genomic alterations have potential value in modulating tumor immunity and immune infiltration patterns. Differences in TME and immunotherapy response between stemness subtypes were therefore explored. C2 had a decreased stromal, immune, and ESTIMATE score, and increased tumor purity, suggestive of an immunosuppressive microenvironment (Figure 7A). MCP-counter showed that the amounts of CD8⁺ T cells, cytotoxic lymphocytes, B cells, NK cells, and neutrophils in C1 were considerably higher than those in C2 (Figure 7B). CIBERSORT further showed that naïve B cells, resting CD4 memory T cells, and resting mast cells were more abundant within C1 (Figure 7C). In contrast, M0 macrophages, M2 macrophages, and follicular helper T cells were mainly concentrated within C2 (Figure 7C). These results indicate a distinct infiltration abundance of tumor-resistant and tumor-promoting immune cells separating the stemness subtypes.

Figure 7 TME status and susceptibility to immunotherapy in two stemness subgroups. (A) Differences in stromal, immune, ESTIMATE scores and tumor purity within two stemness subgroups. C1 and C2 on the x-axis represent the two consensus stemness subgroups identified in this study. (B) Differences in the levels in 8 immune cells and 2 stromal cells across two stemness subgroups. (C) Differences in the proportion of 22 immune cells among two stemness subtypes. (D) Differential levels of immune checkpoints between C1 and C2. (E) T cell dysfunction, exclusion, and TIDE scores within two stemness subtypes. C1 and C2 on the x-axis represent the two consensus stemness subgroups identified in this study. (F) mRNAsi levels in responding and non-responding patients. (G) Percentage of responders between the two stemness subtypes. *, P<0.05; **, P<0.01; ***, P<0.001; ****, P<0.0001. mRNAsi, mRNA stemness index; TIDE, Tumor Immune Dysfunction and Exclusion; TME, tumor microenvironment.

We therefore focused on the response of stemness subtypes to immune checkpoint blockade (ICB). Figure 6D shows that there were changes in the levels of immune checkpoint expression between the different stemness subtypes. For example, CD80, LAG3, and ICOS were increased in C2, whereas CD27, BTLA, and PDCD1 were increased in C1 (Figure 7D). It is suggested that stemness subtypes may exhibit different reactions to immunotherapy. Further, TIDE showed that the TIDE value, T cell dysfunction and exclusion values were reduced in C2 versus C1 (Figure 7E). In addition, mRNAsi was substantially increased in responders than in non-responders in patients who received immunotherapy (Figure 7F), with a large proportion of immunotherapy responders in C2 patients (Figure 7G). These results suggest that C2 patients have stronger immunotherapy responses, suggesting ICB might have developmental promise in BC patients with high stemness.

Creation and testing of an mRNAsi-risk-prognostic model

To improve understanding of the link between mRNAsi-related genes and outcomes in BC individuals, a risk model was built from target genes using the TCGA-BRCA cohort. Seven signature genes (BEND5, TNN, PDLIM4, CD24, POP1, PRDX1, PGK1) were selected by the LASSO regression algorithm at λ=5 (Figure S2). A risk model was generated with the minimum error solution after screening (Table 1) and computed to be risk score = −0.0654*BEND5 − 0.1023*TNN − 0.0930* PDLIM4 + 0.0816*CD24 + 0.0995* POP1 + 0.1828* PRDX1 + 0.3399* PGK1. The model was subsequently used to compute risk scores for the TCGA-BRCA cohort and two external validation cohorts (GSE20685 and METABRIC). We observed in the TCGA-BRCA cohort that receiver operating characteristic (ROC) analysis revealed better survival prediction at 1, 3, and 5 years, with an area under the curve (AUC) ≥0.6 (Figure 8A). And based on the ideal threshold criteria, the TCGA-BRCA cohort was divided into high- and low-risk groups. The survival curves showed that the outcomes of the high-risk group were worse (Figures 8B). Furthermore, we also observed the same results in two external validation cohorts (Figure 8C-8F). These findings suggest the risk score model related to mRNAsi has good predictability for the prognosis of BC patients.

Figure 8 To establish and test the mRNAsi-related risk model. (A-F) ROC curves and survival curves in TCGA, GSE20685 and METABRIC cohorts by risk scores. AUC, area under curve; mRNAsi, mRNA stemness index; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas.

Building and assessing the nomogram

After analyzing risk scores and clinicopathological parameters of BC patients, it was noted that although risk scores did not correlate with lymph node metastasis, increased risk scores were found to be associated with higher tumor infiltration, distant metastasis, and advanced pathological stage (Figure 9A). PR-negative, ER-negative, and HER2-positive patients had higher risk scores compared to controls (Figure 9B). In addition, risk scores were much lower in luminal subtype patients compared to HER2-overexpressing and triple-negative patients (Figure 9B). However, the risk scores were not statistically meaningful between HER2-overexpressing and triple-negative patients (Figure 9B). These observations indicate that risk scores are closely connected to the clinicopathological parameters of BC patients, with those with elevated risk scores exhibiting poorer clinicopathological characteristics.

Figure 9 Nomogram development and validation. (A) Risk scores for tumor size, lymph node status, distant metastasis and AJCC stage. (B) Risk scores between ER, PR and HER2 molecular states and BC molecular subtypes. (C) Nomogram built from risk scores and pathological parameters to predict 3- and 5-year survival in TCGA breast cancer patients. (D) Calibration curves depicting the results of a nomogram predicting 3- and 5-year survival and true survival. *, P<0.05; **, P<0.01; ***, P<0.001. AJCC, American Joint Committee on Cancer; BC, breast cancer; ER, estrogen receptor; HER2, human epidermal growth factor receptor 2; M, distant metastasis; N, lymph node status; OS, overall survival; PR, progesterone receptor; T, tumor size.

Univariate Cox analysis showed that tumor infiltration, lymph node metastasis, distant metastasis, pathological stage, and risk score were all strongly linked with the outcome in BC patients, with an HR risk score of 3.38 (P<0.001) (Table 2). In order to ensure robust independence assessment while avoiding multicollinearity, the risk score and basic variables (including tumor infiltration, lymph node metastasis, and distant metastasis) were incorporated into the multivariate Cox regression analysis, rather than the composite American Joint Committee on Cancer (AJCC) stage. Multivariate Cox analysis demonstrated that the risk score could be an independent prognostic factor [hazard ratio (HR) =3.58, P<0.001] (Table 2). The nomogram, combined with the risk score, age, tumor infiltration, lymph node metastasis, and distant metastasis can be used to forecast the survival rate of BC individuals at different periods (Figure 9C). After verification by the Bootstrap method, the adjusted C-index of this model was 0.76, indicating that it has a good ability to distinguish. Moreover, the calibration curves in the 3- and 5-year periods were well-fitted (Figure 9D), showing that the nomogram has excellent predictability and can help in the clinical detection and management of BC patients.

Candidate drug screening in patients from different risk groups

At present, systemic chemotherapy and targeted therapy are the main therapeutic strategies for BC patients. Therefore, it is crucial to assess the susceptibility of BC patients to different risks of common therapeutic agents. These results showed that paclitaxel, lovastatin, lapatinib, docetaxel, decitabine, and olaparib exhibited higher sensitivity for high-risk individuals (Figure 10A-10F). Conversely, the half-maximal inhibitory concentration (IC50) of tamoxifen, cyclophosphamide and fulvestrant was decreased for low-risk individuals (Figure 10G-10I). These observations suggest different risk groups show different sensitivities to BC therapeutic agents.

Figure 10 Vulnerability of different risk groups to popular therapeutic agents. IC50 estimates for paclitaxel (A), lovastatin (B), lapatinib (C), docetaxel (D), decitabine (E), olaparib (F), tamoxifen (G), cyclophosphamide (H) and fulvestrant (I) in different risk groups. The x-axis represents the TCGA-BRCA cohort stratified into ‘High risk’ and ‘Low risk’ based on the optimal cut-off value of the risk score. IC50, half maximal inhibitory concentration; TCGA, The Cancer Genome Atlas.

Discussion

BC is a female malignancy with high incidences of morbidity and death. Current strategies to prevent and treat this malignant tumor still face great challenges of tumor recurrence, metastasis and drug resistance (3). In recent years, CSCs have received much attention in the cancer field. CSCs are at the root of unlimited tumor proliferation and susceptibility to recurrence (35). Previous studies have implicated CSCs in several malignant phenotypes, such as proliferation and invasion of BC cells, and are known to promote tumor recurrence metastasis and treatment resistance (7,8,36). BC is characterized by substantial molecular and clinical heterogeneity, wherein molecular subtyping is strongly associated with distinct clinicopathological features (37). In this context, the identification of stem cell-like subtypes through the lens of BCSCs and the construction of corresponding prognostic models represent a critical direction for developing novel therapeutic strategies. Supporting this approach, a systematic review and meta-analysis by Zhang et al. demonstrated that integrating multi-omics data—such as genomics and radiomics—with machine learning significantly improves the accuracy and robustness of molecular subtype classification (38). Their findings establish a solid methodological basis for advancing more refined subtyping frameworks rooted in stemness characteristics.

At present, the mRNAsi has been developed through computational biology as an indicator to measure the stemness in tumor patients, which can be effectively used to mine genes associated with tumor stemness (18,19). The mRNAsi is increasingly being used for cancer typing and prognosis. Taking hepatocellular carcinoma as an example, the mRNAsi from cancer tissue can be used to identify novel molecular subtypes of hepatocellular carcinoma with differing prognoses (39). The above reviews serve as key references for the identification of BC stemness subtypes and the construction of mRNAsi-related risk models.

This study first analyzed the correlation between mRNAsi and clinical and immune characteristics, and found that tumor tissues have higher stemness than normal tissues. As tumors progress, the tumor stemness shows an increasing trend. In addition, enhanced tumor stemness is characterized by an accumulation of tumor cells and a reduction of stromal and immune cells in the TME. This is consistent with previous findings (40). Several gene modules closely related to the mRNAsi were identified by WGCNA. Signaling pathway analysis suggested that these gene modules may play different functions in BC. For example, the GPCR ligand-binding and nuclear receptor signaling pathways were the most enriched in the magenta module. GPCRs are the largest known family of cell surface signaling receptors. They affect abnormal cell proliferation by activating several signaling pathways, including the MAPKs, AKT/mTOR, and Hippo pathways (41,42). GPCRs also promote tumor invasion and metastasis by altering the cytoskeleton and activating Rho GTPases. Similarly, members of the nuclear hormone receptor family (e.g., ER, PR), which act as transcription factors when bound to ligands, play a vital part in the progression of BC. Genes involved in ECM-related pathways were mostly concentrated in the brown module. The involvement of tumor-associated extracellular matrix promotes cancer cell growth, metastasis and thus resistance to cell death (43,44). In addition, extracellular matrix-derived mechanical forces can control the transition between stemness and quiescence in BC cells (45). Genes in cell cycle-related pathways were substantially accumulated in the green module. Cell cycle dysregulation is a hallmark of cancer. Abnormalities in this pathway may underlie the maintenance of rapid proliferation and invasion of tumor stem cells. Cyclins are involved not only in BC generation but also in the stem cell-like cellular behavior of BC cells. Inhibition of the cyclins CDK4 and CDK6 reduces the ability of stem cells to migrate in BC (46). In addition, BCSCs are thought to be in a quiescent state resulting from cell cycle dysregulation (47). Therefore, these pathways are closely linked to cancer occurrence. Eleven target genes were subsequently chosen by univariate Cox regression analysis to identify BC stemness subtypes. The subtypes involve BEND5, TNN, PDLIM4, CD24, SHCBP1, MTHFD2, POP1, STIP1, PRDX1, RAD54B and PGK1. Among them, BEND5, TNN, and PDLIM4 can be used as prognostic protective factors for BC. BEND5, as a tumor suppressor gene, inhibits Notch signaling by blocking RBPJ/CSL transcription factor interaction, thereby suppressing BC growth and metastasis (48). Similarly, PDLIM4 is a tumor suppressor. PDLIM4 is an intermediate adaptor protein in BC cells and is involved in the assembly of intracellular stress fibers, remodeling of the microfilament cytoskeleton and epithelial-mesenchymal transition (49). The other genes are risk factors related to tumor prognosis, mostly identified in previous research. For example, CD24 is a widely accepted cell surface marker for BC stem cells. Its low expression or absence and high CD44 levels correlate with stem cell properties. CD24 not only forecasts the prognosis of triple-negative BC but also promotes immune evasion from triple-negative BC by controlling PD-L1 levels (10,12). SHCBP1, an important protein-coding gene for intracellular signaling and cell division, can promote BC cell proliferation by inhibiting CXCL2 expression (50). MTHFD2 with high expression is closely connected to BC prognosis (51). Furthermore, it was demonstrated that MTHFD2 activates the AKT signaling pathway, thereby promoting BC growth (52). POP1 is upregulated in BC and correlates with a worse prognosis (53). It was found that POP1 binds to the RNA component of telomerase and activates the telomerase complex, thereby promoting the growth of BC (54). In addition, high expression of STIP1 (55), PRDX1 (56), RAD54B (57) and PGK1 (58) promotes tumor emergence and are related to poor prognosis in BC patients. Most target genes identified here to be connected to BC stemness strongly correlate with tumor progression or cellular metabolism. It is feasible to use the stemness subtypes, identified by the above genes and the prognostic model constructed by a further screening of the characteristic genes (BEND5, TNN, PDLIM4, CD24, POP1, PRDX1 and PGK1) through machine learning, for clinical prognostic guidance.

In the BC stemness subtype model constructed in this study, patients with higher stemness characteristics tended to exhibit poorer prognosis. Further exploration revealed that the difference in survival between C2 and C1 may be due to the activities of multiple cell cycle and cell proliferation-related signaling pathways in C2. For example, activation of the PI3K/AKT/MTOR pathway contributes to resistance to BC cell death, leading to disease progression (59). In addition, the difference in survival rate between stemness subtypes may also be related to TME status. Compared with C2, the C1 subgroup showed active TME status, such as higher CD8⁺ T cells, B cells, and NK cell infiltration, and the highest stromal score and immune score. Among them, CD8⁺ T cells can detect and eliminate cancer cells, activate cytotoxic T lymphocytes in the tumor immune circulatory system and mediate an anti-tumor response (60). At the same time, tumor-infiltrating B cells stimulate T cell defenses by secreting immunoglobulin, and directly kill cancer cells, thus inhibiting tumor progression (61). In clinical studies, NK cells can indirectly inhibit tumor growth by secreting various cytokines in addition to their killing function (62). Moreover, NK cells tend to work together with CD8⁺ T cells in anti-tumor immunity (62). Both have similar killing mechanisms. Therefore, the present study suggests that the upregulation of infiltration levels of various anti-tumor immune cells may be the cause of the better prognosis of C1 patients. Moreover, our study also supports the idea that BCSCs are highly associated with the TME. Immunotherapy is a new anti-cancer method developed in recent years. Normally, immune cells in the body can identify, kill and destroy abnormal cells, including cancer cells (63). ICB is one of the immunotherapeutic tools, mostly used in combination with chemotherapy drugs. Notably, one of the reasons for the poor efficacy of ICB and its combination with chemotherapy drugs in clinical treatment is the inability of existing means to accurately screen immunotherapy-sensitive populations. It is proposed that TMB triggers an effective anti-tumor immune response in the body, resulting in long-lasting clinical efficacy of immunotherapy (64). This study suggests that ICB may contribute to individualized precision therapy for patients with high stemness BC in clinical practice.

Inevitably, there are certain restrictions to this study. Firstly, it should be acknowledged that the present analysis was conducted exclusively using retrospective data from public databases (TCGA, GEO, METABRIC). Although these cohorts provide valuable resources, this reliance inevitably introduces certain limitations. Specifically, the investigators had no control over sample collection procedures, sequencing quality, or the consistency of clinical data annotation. In addition, the lack of detailed treatment information and the limited availability of long-term follow-up data in some datasets may introduce potential confounding biases in the analyses of subtypes and prognostic outcomes. Therefore, future prospective multicenter studies with standardized protocols and comprehensive clinical metadata will be essential to validate and further refine our findings. Secondly, the stemness-related genes identified in this study (e.g., BEND5 and CD24) were derived from comprehensive bioinformatics analyses. While this strategy provides a theoretical framework for hypothesis generation, further experimental validation is necessary to clarify the biological functions and causal roles of these genes in regulating BCSC properties. The absence of in vitro assays (such as spheroid formation and invasion assays following gene knockout or overexpression in BC cell lines) and in vivo experiments (such as tumorigenicity studies in xenograft models) represents an important limitation. Nevertheless, by systematically screening candidate genes and molecular subtypes, this study lays a foundation for future investigations. Subsequent mechanistic studies will be required to substantiate the computational predictions, which may ultimately uncover novel therapeutic targets for overcoming stemness-driven treatment resistance and disease recurrence. Thirdly, the focus on mRNAsi, whilst providing a reproducible and well-validated framework (18), inherently fails to encompass epigenetic regulation (such as the methylation DNA stemness index) or the spatial heterogeneity of intratumoural stemness. The future integration of multi-omics stemness measurements with single-cell or spatial transcriptomics studies will be crucial for validating our subtype classification and constructing more comprehensive, clinically precise models.

Conclusions

In sum, this study categorized BC in terms of mRNAsi-related genes. Differences in prognosis, biological pathways, TME, and immunotherapy response were found between the stemness subtypes. Additionally, the established mRNAsi-related risk model worked well in forecasting BC prognosis. This investigation provided new possibilities for diagnosis and clinical treatment strategies for BC.

Acknowledgments

We are thankful to Prof. Stanley Lin for his critical and careful editing and proofreading of the manuscript.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-837/rc

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-837/prf

Funding: This study was supported by the National Natural Science Foundation of China (No. 82273457), the Guangdong Basic and Applied Basic Research Foundation (No. 2023A1515012762), and Science and Technology Special Project of Guangdong Province (No. 210715216902829).

Conflics of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-2025-837/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Bray F, Laversanne M, Sung H, et al. Global cancer statistics 2022: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 2024;74:229-63. [Crossref] [PubMed]
2. Tsang JYS, Tse GM. Molecular Classification of Breast Cancer. Adv Anat Pathol 2020;27:27-35. [Crossref] [PubMed]
3. Trayes KP, Cokenakes SEH. Breast Cancer Treatment. Am Fam Physician 2021;104:171-8.
4. Lu J, Fang Q, Liu X, et al. Construction of molecular subtypes and prognostic model for breast cancer based on sulfur metabolism-related genes. Transl Cancer Res 2025;14:3452-70. [Crossref] [PubMed]
5. Fang Z, Liao SC, Guo YY, et al. Development of coenzyme Q10-related molecular subtypes and a prognostic signature for predicting breast cancer prognosis and response to immunotherapy. Transl Cancer Res 2025;14:2010-28. [Crossref] [PubMed]
6. Dey P, Rathod M, De A. Targeting stem cells in the realm of drug-resistant breast cancer. Breast Cancer (Dove Med Press) 2019;11:115-35. [Crossref] [PubMed]
7. Xu H, Zhang F, Gao X, et al. Fate decisions of breast cancer stem cells in cancer progression. Front Oncol 2022;12:968306. [Crossref] [PubMed]
8. Butti R, Gunasekaran VP, Kumar TVS, et al. Breast cancer stem cells: Biology and therapeutic implications. Int J Biochem Cell Biol 2019;107:38-52. [Crossref] [PubMed]
9. Ibragimova M, Tsyganov M, Litviakov N. Tumour Stem Cells in Breast Cancer. Int J Mol Sci 2022;23:5058. [Crossref] [PubMed]
10. Escudero Mendez L, Srinivasan M, Hamouda RK, et al. Evaluation of CD44+/CD24- and Aldehyde Dehydrogenase Enzyme Markers in Cancer Stem Cells as Prognostic Indicators for Triple-Negative Breast Cancer. Cureus 2022;14:e28056. [Crossref] [PubMed]
11. Panigoro SS, Kurnia D, Kurnia A, et al. ALDH1 Cancer Stem Cell Marker as a Prognostic Factor in Triple-Negative Breast Cancer. Int J Surg Oncol 2020;2020:7863243. [Crossref] [PubMed]
12. Zhu X, Yu J, Ai F, et al. CD24 May Serve as an Immunotherapy Target in Triple-Negative Breast Cancer by Regulating the Expression of PD-L1. Breast Cancer (Dove Med Press) 2023;15:967-84. [Crossref] [PubMed]
13. Liu C, Qiang J, Deng Q, et al. ALDH1A1 Activity in Tumor-Initiating Cells Remodels Myeloid-Derived Suppressor Cells to Promote Breast Cancer Progression. Cancer Res 2021;81:5919-34. [Crossref] [PubMed]
14. Ciummo SL, D'Antonio L, Sorrentino C, et al. The C-X-C Motif Chemokine Ligand 1 Sustains Breast Cancer Stem Cell Self-Renewal and Promotes Tumor Progression and Immune Escape Programs. Front Cell Dev Biol 2021;9:689286. [Crossref] [PubMed]
15. Radharani NNV, Yadav AS, Nimma R, et al. Tumor-associated macrophage derived IL-6 enriches cancer stem cell population and promotes breast tumor progression via Stat-3 pathway. Cancer Cell Int 2022;22:122. [Crossref] [PubMed]
16. Huang TX, Guan XY, Fu L. Therapeutic targeting of the crosstalk between cancer-associated fibroblasts and cancer stem cells. Am J Cancer Res 2019;9:1889-904.
17. Guha A, Goswami KK, Sultana J, et al. Cancer stem cell-immune cell crosstalk in breast tumor microenvironment: a determinant of therapeutic facet. Front Immunol 2023;14:1245421. [Crossref] [PubMed]
18. Malta TM, Sokolov A, Gentles AJ, et al. Machine Learning Identifies Stemness Features Associated with Oncogenic Dedifferentiation. Cell 2018;173:338-354.e15. [Crossref] [PubMed]
19. Zheng H, Song K, Fu Y, et al. An absolute human stemness index associated with oncogenic dedifferentiation. Brief Bioinform 2021;22:2151-60. [Crossref] [PubMed]
20. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 2008;9:559. [Crossref] [PubMed]
21. Wu T, Hu E, Xu S, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb) 2021;2:100141. [Crossref] [PubMed]
22. Rothfels K, Milacic M, Matthews L, et al. Using the Reactome Database. Curr Protoc 2023;3:e722. [Crossref] [PubMed]
23. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics 2010;26:1572-3. [Crossref] [PubMed]
24. Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015;43:e47. [Crossref] [PubMed]
25. Ashburner M, Ball CA, Blake JA, et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet 2000;25:25-9. [Crossref] [PubMed]
26. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res 2000;28:27-30. [Crossref] [PubMed]
27. Subramanian A, Tamayo P, Mootha VK, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A 2005;102:15545-50. [Crossref] [PubMed]
28. Yoshihara K, Shahmoradgoli M, Martínez E, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun 2013;4:2612. [Crossref] [PubMed]
29. Becht E, Giraldo NA, Lacroix L, et al. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol 2016;17:218. [Crossref] [PubMed]
30. Newman AM, Liu CL, Green MR, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods 2015;12:453-7. [Crossref] [PubMed]
31. Liu Y, He M, Wang D, et al. HisgAtlas 1.0: a human immunosuppression gene database. Database (Oxford) 2017;2017:bax094. [Crossref] [PubMed]
32. Jiang P, Gu S, Pan D, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med 2018;24:1550-8. [Crossref] [PubMed]
33. Blanche P, Dartigues JF, Jacqmin-Gadda H. Estimating and comparing time-dependent areas under receiver operating characteristic curves for censored event times with competing risks. Stat Med 2013;32:5381-97. [Crossref] [PubMed]
34. Maeser D, Gruener RF, Huang RS. oncoPredict: an R package for predicting in vivo or cancer patient drug response and biomarkers from cell line screening data. Brief Bioinform 2021;22:bbab260. [Crossref] [PubMed]
35. Nassar D, Blanpain C. Cancer Stem Cells: Basic Concepts and Therapeutic Implications. Annu Rev Pathol 2016;11:47-76. [Crossref] [PubMed]
36. Zhang L, Chen W, Liu S, et al. Targeting Breast Cancer Stem Cells. Int J Biol Sci 2023;19:552-70. [Crossref] [PubMed]
37. Abousahmeen A, Saud MAB, Elrais S, et al. Molecular subtypes and clinicopathological features of breast cancer in Libya: a first glance. Ther Radiol Oncol 2023;7:15.
38. Zhang Y, Li G, Bian W, et al. Value of genomics- and radiomics-based machine learning models in the identification of breast cancer molecular subtypes: a systematic review and meta-analysis. Ann Transl Med 2022;10:1394. [Crossref] [PubMed]
39. Chen D, Liu J, Zang L, et al. Integrated Machine Learning and Bioinformatic Analyses Constructed a Novel Stemness-Related Classifier to Predict Prognosis and Immunotherapy Responses for Hepatocellular Carcinoma Patients. Int J Biol Sci 2022;18:360-73. [Crossref] [PubMed]
40. Zhao X, Lin J. Construction and Validation of a Prognostic Model Based on mRNAsi-Related Genes in Breast Cancer. Comput Math Methods Med 2022;2022:6532591. [Crossref] [PubMed]
41. Chaudhary PK, Kim S. An Insight into GPCR and G-Proteins as Cancer Drivers. Cells 2021;10:3288. [Crossref] [PubMed]
42. Lappano R, Jacquot Y, Maggiolini M. GPCR Modulation in Breast Cancer. Int J Mol Sci 2018;19:3840. [Crossref] [PubMed]
43. Zhao Y, Zheng X, Zheng Y, et al. Extracellular Matrix: Emerging Roles and Potential Therapeutic Targets for Breast Cancer. Front Oncol 2021;11:650453. [Crossref] [PubMed]
44. Yu TY, Zhang G, Chai XX, et al. Recent progress on the effect of extracellular matrix on occurrence and progression of breast cancer. Life Sci 2023;332:122084. [Crossref] [PubMed]
45. Li C, Qiu S, Liu X, et al. Extracellular matrix-derived mechanical force governs breast cancer cell stemness and quiescence transition through integrin-DDR signaling. Signal Transduct Target Ther 2023;8:247. [Crossref] [PubMed]
46. Lamb R, Lehn S, Rogerson L, et al. Cell cycle regulators cyclin D1 and CDK4/6 have estrogen receptor-dependent divergent functions in breast cancer migration and stem cell-like activity. Cell Cycle 2013;12:2384-94. [Crossref] [PubMed]
47. De Angelis ML, Francescangeli F, Zeuner A. Breast Cancer Stem Cells as Drivers of Tumor Chemoresistance, Dormancy and Relapse: New Challenges and Therapeutic Opportunities. Cancers (Basel) 2019;11:1569. [Crossref] [PubMed]
48. Shi Y, Zhang D, Chen J, et al. Interaction between BEND5 and RBPJ suppresses breast cancer growth and metastasis via inhibiting Notch signaling. Int J Biol Sci 2022;18:4233-44. [Crossref] [PubMed]
49. Kravchenko DS, Ivanova AE, Podshivalova ES, et al. PDLIM4/RIL-mediated regulation of Src and malignant properties of breast cancer cells. Oncotarget 2020;11:22-30. [Crossref] [PubMed]
50. Yu X, Feng G, Nian R, et al. SHCBP1 Promotes the Proliferation of Breast Cancer Cells by Inhibiting CXCL2. J Cancer 2023;14:3444-56. [Crossref] [PubMed]
51. Liu F, Liu Y, He C, et al. Increased MTHFD2 expression is associated with poor prognosis in breast cancer. Tumour Biol 2014;35:8685-90. [Crossref] [PubMed]
52. Huang J, Qin Y, Lin C, et al. MTHFD2 facilitates breast cancer cell proliferation via the AKT signaling pathway. Exp Ther Med 2021;22:703. [Crossref] [PubMed]
53. He X, Wang J, Yu H, et al. Clinical significance for diagnosis and prognosis of POP1 and its potential role in breast cancer: a comprehensive analysis based on multiple databases. Aging (Albany NY) 2022;14:6936-56. [Crossref] [PubMed]
54. Zhu M, Wu C, Wu X, et al. POP1 promotes the progression of breast cancer through maintaining telomere integrity. Carcinogenesis 2023;44:252-62. [Crossref] [PubMed]
55. Wu R, Liu F, Peng P, et al. Tumor stress-induced phosphoprotein 1 as a prognostic biomarker for breast cancer. Ann Transl Med 2018;6:302. [Crossref] [PubMed]
56. Mei J, Hao L, Liu X, et al. Comprehensive analysis of peroxiredoxins expression profiles and prognostic values in breast cancer. Biomark Res 2019;7:16. [Crossref] [PubMed]
57. Zhang Z, Li X, Han Y, et al. RAD54B potentiates tumor growth and predicts poor prognosis of patients with luminal A breast cancer. Biomed Pharmacother 2019;118:109341. [Crossref] [PubMed]
58. Li Y, Wang S, Zhang X, et al. Expression Characteristics and Significant Prognostic Values of PGK1 in Breast Cancer. Front Mol Biosci 2021;8:695420. [Crossref] [PubMed]
59. Miricescu D, Totan A, Stanescu-Spinu II, et al. PI3K/AKT/mTOR Signaling Pathway in Breast Cancer: From Molecular Landscape to Clinical Aspects. Int J Mol Sci 2020;22:173. [Crossref] [PubMed]
60. Farhood B, Najafi M, Mortezaee K. CD8(+) cytotoxic T lymphocytes in cancer immunotherapy: A review. J Cell Physiol 2019;234:8509-21. [Crossref] [PubMed]
61. Li M, Quintana A, Alberts E, et al. B Cells in Breast Cancer Pathology. Cancers (Basel) 2023;15:1517. [Crossref] [PubMed]
62. Wolf NK, Kissiov DU, Raulet DH. Roles of natural killer cells in immunity to cancer, and applications to immunotherapy. Nat Rev Immunol 2023;23:90-105. [Crossref] [PubMed]
63. Barzaman K, Moradi-Kalbolandi S, Hosseinzadeh A, et al. Breast cancer immunotherapy: Current and novel approaches. Int Immunopharmacol 2021;98:107886. [Crossref] [PubMed]
64. Anagnostou V, Bardelli A, Chan TA, et al. The status of tumor mutational burden and immunotherapy. Nat Cancer 2022;3:652-6. [Crossref] [PubMed]