Tumor tissue samples were collected from eight patients with BC at The First Affiliated Hospital of Nanjing Medical University, with matched paracancerous tissues from eight healthy individuals serving as controls. All specimens underwent RNA-seq and stringent quality control (QC) to form the whole transcriptome RNA-seq dataset. The study was approved by the Ethics Committee of The First Affiliated Hospital of Nanjing Medical University, and informed consent was obtained from all participants.

To enhance analytical robustness, additional BC-related datasets were sourced from public repositories. Specifically, RNA-seq data comprising 1,082 BC tumor samples and 113 paracancerous samples were retrieved from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/) for expression validation. Concurrently, the scRNA-seq dataset GSE245601, based on the GPL18573 platform, was downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/gds), including 10 BC tumor samples and 2 normal tissue samples for cellular-level expression analysis (34).

Stem cell-associated gene data were curated by referencing established literature (35), and 26 stem cell gene sets were extracted from the StemChecker database (http://stemchecker.sysbiolab.eu/), yielding 4,419 SCRGs. Additionally, 218 TMMRGs were incorporated based on published literature (13). The overall analytical workflow of this study was illustrated in Figure 1 .

Total RNA from the 8 tumor and 8 paracancerous tissue samples was isolated using TRIzol reagent (Invitrogen, CA, USA). RNA integrity and concentration were evaluated using a NanoDrop ND-1000 spectrophotometer (Wilmington, DE, USA) and an Agilent Bioanalyzer 2100 system (Agilent, CA, USA). Sequencing libraries were prepared with the Hieff NGS Ultifaillumina Dual-mode mRNA Library Prep Kit to generate inserts of 300 ± 50 bp, followed by high-throughput paired-end sequencing (PE150) on the Illumina NovaSeq 6000 platform.

Raw sequencing data in FASTQ format underwent QC using fastp software (v0.23.4) (36) with default parameters, including adapter trimming, deduplication, and removal of low-quality reads. Cleaned reads were aligned to the human reference genome (GRCh38) using HISAT2 (v2.2.1) (37), and alignment outputs were stored in BAM format. Transcript assembly and quantification were performed using StringTie (v2.2.0) (38), and expression levels were normalized to FPKM, calculated as:

Differential expression analysis between tumor and normal tissues was conducted on the transcriptome-wide RNA-seq data using the DESeq2 package (v1.38.0) (39). Differentially expressed genes (DEGs) were defined by |log₂ fold change (FC)| > 0.5 and P < 0.05. Visualization of DEGs was performed using a volcano plot and heatmap, generated with ggplot2 (v3.3.6) (40) and ComplexHeatmap (v2.14.0) (41), respectively. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using the clusterProfiler package (v4.7.1.3) (42), with significance set at P < 0.05. GO terms were classified into three categories: biological processes (BPs), molecular functions (MFs), and cellular components (CCs).

The intersection of DEGs, SCRGs, and TMMRGs was identified using the VennDiagram package (v1.7.1) (43) to isolate candidate genes implicated in both stemness and telomere regulation within the context of BC. These candidate genes were subsequently analyzed using the Search Tool for the Retrieval of Interacting Genes (STRING) database (http://string-db.org) with a confidence threshold set at > 0.7, and a protein–protein interaction (PPI) network was generated to elucidate potential molecular interactions at the protein level. The Molecular Complex Detection (MCODE) plugin was used to select the highest-scoring cluster and identify hub genes. A PPI network was subsequently constructed for these hub genes. Both networks were generated using Cytoscape software (v 3.9.1) (44).

To assess the diagnostic potential of the hub genes, Least Absolute Shrinkage and Selection Operator (LASSO) regression analysis was performed using the glmnet package (v4.1.4) (45), configured with the ‘binomial’ family parameter. The optimal model was established at the point of minimal model error and the lowest lambda value, at which the genes as candidate biomarkers were selected for expression analysis. These genes were then subjected to expression analysis in both the in-house whole transcriptome RNA-seq dataset and the TCGA-BC dataset, comparing tumor and normal tissues. Expression differences were visualized via box plots generated with the ggplot2 package, and genes exhibiting consistent and statistically significant differential expression across both datasets (P < 0.05) were selected as biomarkers.

To investigate the functional relevance of these biomarkers, enrichment analyses were conducted on the whole transcriptome dataset. Using the Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb), the curated gene set collection c2.cp.v2023.2.Hs.symbols.gmt was employed. For each biomarker, Spearman correlation coefficients were computed with all other genes, and genes were ranked in descending order of correlation. GSEA was performed using the clusterProfiler package, with statistical thresholds set at P < 0.05 and False Discovery Rate (FDR) < 0.05. The top five significantly enriched pathways per biomarker were visualized using the enrichplot package (v1.18.4) (46), ranked by enrichment significance.

To further explore the differences in the immune microenvironment between patients with BC and normal individuals, a single-sample GSEA (ssGSEA) approach from the GSVA package (v1.42.0) (47) was applied to estimate the infiltration levels of 28 predefined immune cell types (48) within each sample of the whole transcriptome dataset. Comparative analysis revealed statistically significant differences in immune cell abundance between tumor and normal samples (P < 0.05). Additionally, a set of 40 immune checkpoints was compiled from previously published literature (49), and differences in the expression of these checkpoints between tumor and normal samples were analyzed within the same dataset (P < 0.05).

Drug–biomarker interactions were identified using the Drug Gene Interaction Database (DGIDB, https://dgidb.org), with the top five candidates ranked by interaction score. The resulting drug–biomarker network was visualized using Cytoscape.

To investigate the underlying molecular mechanisms of the biomarkers, the DESeq2 package was employed to analyze differentially expressed miRNAs and lncRNAs (DE-miRNAs and DE-lncRNAs) from a whole transcriptome RNA-seq dataset, comparing tumor and normal tissues under the criteria of P < 0.05 and |log₂FC| > 0.5. Volcano plots were generated to depict the distribution of DE-miRNAs and DE-lncRNAs.

Candidate miRNAs targeting the biomarkers were predicted via miRDB (https://mirdb.org/) and TargetScan (https://www.targetscan.org/) using a prediction score threshold of 500,000. Correspondingly, target lncRNAs for these miRNAs were identified through starBase (http://starbase.sysu.edu.cn/index.php), filtered by a clipExpNum > 4.

Subsequent analysis involved intersecting DE-miRNAs with the predicted miRNA targets and DE-lncRNAs with the predicted lncRNA targets to isolate key regulatory miRNAs and lncRNAs. A lncRNA–miRNA–mRNA regulatory network was then constructed in Cytoscape to elucidate the multifaceted regulatory landscape influencing biomarker expression.

To investigate the cellular mechanisms underlying BC and to characterize biomarker heterogeneity at the single-cell level, comprehensive analyses were conducted using the GSE245601 dataset. QC, clustering, and annotation were performed with the Seurat package (v5.0.1) (50). Key quality metrics—including gene count, cell count, and mitochondrial gene percentage—were calculated. Cells with fewer than 200 genes and genes covered by fewer than three cells were filtered out, as were cells expressing either less than 200 or more than 5000 genes. Genes expressed in fewer than three cells or with total counts below 500 or above 20,000 were filtered out. Mitochondrial content was capped at 10%. Distributions of nFeature_RNA, nCount_RNA, and percent.mt before and after filtering were visualized, and only cells meeting all criteria were retained for downstream analysis. Subsequently, data normalization was applied, and the top 2,000 highly variable genes were selected using the vst method. Principal component analysis (PCA) was conducted to evaluate variance across components, with results visualized via a PCA elbow plot. Statistically significant principal components, identified using the JackStrawPlot function based on gene-level P-values, were selected for further analysis. Unsupervised clustering was performed using the FindNeighbors and FindClusters functions, and clusters were visualized using Uniform Manifold Approximation and Projection (UMAP) at a resolution of 0.4. Cell type annotation was based on canonical marker genes reported in prior research (34), and a bubble plot was generated to illustrate the specificity of these markers across identified cell populations. In addition, the NF-κB activation, stemness-related signatures, and ERK-JUN signaling pathway scores of the annotated cell types were evaluated using the AddModuleScore method.

In the scRNA-seq analysis, the AddModuleScore function from the Seurat package was utilized to compute enrichment scores of specific gene sets for each sample in the GSE245601 dataset, thereby assessing metabolic pathway activity across cell populations. A curated background gene set, encompassing candidate genes associated with both stem cells and TMM in BC, served as the basis for this scoring. Scores were calculated across annotated cell types in tumor and normal tissues, and cell populations demonstrating significant differences (P < 0.05) were designated as candidate key cells. Subsequent single-cell-level analysis examined the expression profiles of biomarkers within these candidate key cells across tumor and normal tissues. Cells exhibiting consistent and significant differential expression across all biomarkers (P < 0.05) were defined as key cells relevant to both stemness and telomere maintenance in BC. Additionally, heterogeneity within these key cells was examined. UMAP-based reclustering was performed using FindNeighbors and FindClusters to delineate cellular subtypes, followed by assessment of biomarker expression across the resulting subclusters (P < 0.05).

To infer intercellular interactions, cell–cell communication analysis was conducted by evaluating receptor–ligand expression and pairing patterns among cell types. The distribution and proportions of key cell subtypes relative to other populations in tumor and normal tissues were quantified. Communication networks were constructed using the CellChat package (v1.6.1) (51), generating aggregated intercellular signaling maps and visualizing the contribution of each cell cluster to the overall communication landscape. Developmental trajectories and lineage dynamics of key cells and their subtypes were further explored through pseudo-time analysis using the Monocle package (v2.26.0) (52).

Tumor tissue samples were obtained from five patients diagnosed with BC, alongside five adjacent normal tissue samples from healthy individuals, all collected at The First Affiliated Hospital of Nanjing Medical University. RT-qPCR was performed for experimental validation of biomarker expression. Total RNA was extracted from the ten samples using TRIzol reagent (Ambion, Austin, USA), and cDNA synthesis was carried out using the SureScript First-Strand cDNA Synthesis Kit (Servicebio, Wuhan, China), following manufacturer protocols. qPCR amplification was performed with the 2 × Universal Blue SYBR Green qPCR Master Mix (Servicebio, Wuhan, China). GAPDH served as the internal control, and gene expression was quantified using the 2^−ΔΔCt method (53). Primer sequences are provided in Supplementary Table 1 .

All statistical analyses were conducted in R software (v4.2.2). Group comparisons were performed using the Wilcoxon test, with P < 0.05 considered statistically significant.

Differential expression analysis of the whole transcriptome RNA-seq dataset identified 8,193 DEGs between tumor and normal samples, including 7,014 upregulated and 1,179 downregulated genes in BC ( Figures 2A, B ). These DEGs underwent enrichment analysis, yielding 649 significantly associated GO terms: 404 BPs such as “cell cycle,” “DNA damage response,” and “signal transduction” ( Figure 2C ); 112 CCs, including “plasma membrane,” “extracellular exosome,” and “cell junction” ( Figure 2D ); and 133 MFs, including “protein binding,” “protein kinase activity,” and “DNA-binding transcription factor activity” ( Figure 2E ). Additionally, 68 KEGG pathways enriched in DEGs were identified, including “pathways in cancer,” “PI3K-Akt signaling pathway,” “proteoglycans in cancer,” “MAPK signaling pathway,” and “focal adhesion” ( Figure 2F ). These pathways reflect essential mechanisms contributing to BC pathogenesis, reinforcing the reliability of differential expression results and providing functional insight into molecular alterations associated with BC. By intersecting the 8,193 DEGs with 4,419 SCRGs and 218 TMMRGs, 33 candidate genes related to both stemness and TMM in BC were identified ( Figure 2G ).

A PPI network encompassing 26 nodes and 77 edges was constructed based on these 33 candidate genes using a publicly available database ( Figure 3A ). Further analysis with the MCODE plugin identified the most prominent functional module, comprising 9 hub genes: JUN, HDAC1, HSP90AB1, SP1, ESR1, NFKB1, MYC, STAT3, and TP53. An additional PPI network was generated to visualize the interactions among these hub genes ( Figure 3B ).

Subsequently, the nine hub genes were subjected to LASSO regression analysis. At log(lambda.min) = −7.49352 and log(lambda.1se) = −5.725879, the optimal model retained five genes: JUN, SP1, NFKB1, STAT3, and TP53 ( Figures 3C, D ). Expression profiling across both the whole transcriptome RNA-seq dataset and the TCGA-BC cohort identified JUN, NFKB1, and SP1 as exhibiting statistically significant and consistent expression patterns (P < 0.01) ( Figures 3E, F ). These three genes were subsequently designated as biomarkers associated with both stem cell function and TMM in BC. Notably, all three demonstrated markedly reduced expression in tumor tissues compared to normal counterparts (P < 0.01).

GSEA was performed to elucidate biological pathways associated with the identified biomarkers. Due to the lack of significant pathway enrichment for NFKB1, downstream analysis focused on pathways enriched by JUN and SP1. Enrichment analysis revealed that JUN was predominantly associated with pathways such as “ECM organization,” “RhoA regulation pathway,” “Nrf2 pathway,” “S1P–S1P3 pathway,” and “ECM proteoglycans” ( Figure 4A ). In contrast, SP1 enriched pathways included “selenoamino acid metabolism,” “SRP-dependent cotranslational protein targeting to membrane,” “translation initiation via medicus reference,” “ribosome,” and “eukaryotic translation elongation” ( Figure 4B ). These pathways are involved in key processes including ECM structure and tissue reconstruction, cell signaling regulation, oxidative stress responses, and protein synthesis and transport, all of which may contribute to BC pathogenesis through their interaction with these biomarkers.

Comprehensive analysis of the immune microenvironment in BC revealed distinct immune infiltration patterns, as visualized in a stacked bar chart illustrating the relative abundance of 28 immune cell types across individual samples, computed via the ssGSEA algorithm ( Figure 4C ). Notably, activated B cells and NK cells demonstrated statistically significant differences between tumor and normal tissues (P < 0.05), with both cell types showing elevated infiltration in tumor samples ( Figure 4D ). This heightened immune presence likely reflects immune activation within the tumor microenvironment. Specifically, activated B cells may contribute to antibody-mediated responses, while NK cells are known for their direct cytotoxic activity against malignant cells, indicating that differential immune cell infiltration may be associated with immune escape, tumor growth, and the broader pathophysiology of BC.

Expression profiles were available for 34 out of the 40 immune checkpoints identified in the literature. An analysis of these immune checkpoints in tumor and normal samples revealed significant differences in the expression of CD200, CD274, TIGIT, TNFRSF25, and TNFSF15 between the two groups (P < 0.05) ( Figure 4E ). These five immune checkpoints also exhibited elevated expression levels in tumor samples (P < 0.05). This suggests that these checkpoints may be activated within the tumor microenvironment, potentially aiding tumor cells in evading immune surveillance and clearance, thereby positioning them as promising therapeutic targets. Targeting these immune checkpoints could inhibit tumor growth and metastasis.

Utilizing DGIdb, drug–biomarker interactions were predicted and ranked by interaction scores ( Supplementary Table 2 ). Based on the selected candidates, a drug–biomarker interaction network was constructed, comprising 14 nodes and 11 edges ( Figure 5A ).

Parallel differential expression analysis identified 272 DE-miRNAs (248 upregulated, 24 downregulated) and 12,287 DE-lncRNAs (11,899 upregulated, 388 downregulated) between tumor and normal samples ( Figures 5B, C ). Predictions of miRNAs targeting the biomarkers, along with their corresponding target lncRNAs, were obtained from public databases. Intersections between DE-miRNAs and predicted target miRNAs, as well as between DE-lncRNAs and predicted target lncRNAs, enabled the identification of key miRNAs and lncRNAs. These regulatory elements were integrated into a lncRNA–miRNA–mRNA network consisting of 18 nodes and 31 edges ( Figure 5D ), including 7 miRNAs, 9 lncRNAs, and 2 biomarkers (SP1 and JUN). Within this network, complex regulatory relationships were delineated; for example, XIST may regulate JUN via hsa-miR-139-5p, and potentially modulate SP1 through hsa-miR-199a-3p, hsa-miR-7-5p, hsa-miR-324-5p, and hsa-miR-299-5p.

Following transcriptomic analysis, single-cell studies provided insights into the unique expression patterns of biomarkers across various cell types and investigated the underlying cellular mechanisms of BC. Data QC was initially performed on the GSE245601 dataset, with metrics for nFeature_RNA, nCount_RNA, and percent.mt visualized before and after QC ( Supplementary Figures 1A, B ). Initially, the dataset contained 69,275 cells and 25,558 genes, which were reduced to 55,204 cells post-filtering, retaining the same number of genes. The top 2,000 highly variable genes were selected for further analysis ( Supplementary Figure 1C ). PCA demonstrated that cells from both tumor and normal samples were concentrated, with minimal variation beyond the 30th principal component, leading to the selection of the first 30 components (dims = 30) and a clustering resolution of 0.4 ( Supplementary Figures 1D, E ). UMAP clustering revealed 21 distinct clusters ( Figures 6A, B ). Marker gene analysis identified nine cell types: plasmablasts, mast cells, B cells, myeloid cells, basal cells, endothelial cells, mesenchymal cells, T cells, and luminal cells ( Figure 6C ). The specificity of these marker genes was validated using a bubble plot, enabling precise annotation of cell types ( Figure 6D ). The distribution of these nine cell types in tumor and normal samples was also visualized, revealing similar clustering patterns between groups but with higher cell counts in tumors across all types compared to normal samples ( Supplementary Figure 2 ). These differences are likely attributed to the tumor microenvironment and the activation of immune responses in patients with BC. Additionally, basal cells had higher scores in the ERK-JUN signaling and stemness-related signatures pathways, while myeloid cells had the highest score in the NF-κB activation pathway ( Supplementary Figure 3 ). This suggested that the proliferation, differentiation, and stemness regulation of basal cells, as well as the inflammatory response of myeloid cells, might play important roles in the biological processes of BC.

Subsequent analyses employed the AddModuleScore function, using the previously identified 33 candidate genes associated with stem cells and TMM in BC as the background gene set. The scores for the nine annotated cell types in tumor and normal samples were calculated and compared. Noteworthy differences in AddModuleScore were observed among endothelial cells, luminal cells, mesenchymal cells, and myeloid cells (P < 0.01), indicating alterations in cell function and state within the tumor microenvironment. These cells were thus classified as candidate key cells ( Figure 7A ). Specifically, endothelial, luminal, and mesenchymal cells exhibited significantly lower scores in tumor samples (P < 0.0001), suggesting suppression of their stem cell properties and TMM within the tumor microenvironment. This suppression may be linked to tumor aggressiveness, metastatic potential, or treatment resistance. In contrast, myeloid cells displayed significantly higher scores in tumor samples (P < 0.01), potentially reflecting their active role in supporting tumor growth and maintaining the tumor microenvironment, such as by promoting inflammation, enhancing tumor cell survival, or suppressing immune responses.

Further investigation of the expression patterns of JUN, NFKB1, and SP1 in these key cell types revealed significant differences in both mesenchymal and luminal cells (P < 0.01) ( Figures 7B–D ).

Subsequent investigations examined the heterogeneity of mesenchymal cells, with UMAP clustering categorizing these cells into three subtypes: MSC1, MSC2, and MSC3 ( Figure 8a ). The distribution of these subtypes in tumor and normal tissues revealed that, compared to normal samples, MSC1 and MSC2 were more abundant in BC, while MSC3 was less prevalent ( Figure 8b ). A bubble plot further emphasized the specificity of marker genes across the three subtypes ( Figure 8c ). Biomarker expression was then evaluated across these subtypes. In tumor samples, JUN exhibited significantly lower expression in MSC1 and MSC2 (P < 0.05), while NFKB1 was significantly reduced in MSC2 and MSC3 (P < 0.01) ( Figures 8d, e ). SP1, however, showed no significant expression differences across the three subtypes (P > 0.05) ( Figure 8f ).

The study further assessed the distribution and abundance of three mesenchymal cell subtypes, alongside eight other cell types, in both tumor and normal tissues, revealing significant variations in cell type proportions between these environments ( Figure 9A ). These differences suggested that BC may influence the interaction dynamics among various cell populations. Subsequently, communication analysis was conducted to evaluate interactions between MSC1, MSC2, and MSC3, as well as the other eight cell types, under both normal and tumor conditions, focusing on both the frequency and strength of these interactions. Under BC conditions, the number of interactions between the mesenchymal subtypes and other cell types was generally reduced compared to normal tissues ( Figure 9B ), implying that the tumor microenvironment may alter cell functionality and behavior by diminishing communication between mesenchymal cells and other cell types.

In terms of interaction strength, compared to normal conditions, communication between MSC1 and endothelial cells was weakened in BC, while interactions with B cells and myeloid cells were enhanced. Similarly, MSC2 showed reduced interaction with endothelial cells but increased communication with basal and myeloid cells. MSC3 exhibited similar patterns, with decreased communication with endothelial cells and increased interaction with myeloid cells ( Figure 9C ). These findings may reflect alterations in immune regulation within the tumor microenvironment, potentially contributing to tumor invasiveness and environmental changes. The results suggest that BC could reorganize functional and interaction patterns among cells, particularly between mesenchymal and immune cells, which may be critical for understanding tumor progression and developing effective therapeutic strategies. Additionally, ligand-receptor interactions among the mesenchymal cell subtypes and other cell types under both normal and BC conditions were also examined ( Figure 10 ), highlighting the importance of these interactions in regulating cellular behavior, maintaining tissue function, and responding to environmental changes.

In the final stages of the study, pseudo-time analysis was performed on mesenchymal cells to trace their developmental trajectory. The analysis identified a clear starting point, from which cells underwent maturation as they progressed along the trajectory. Four differentiation nodes were detected, signaling points where individual cell clusters began to diverge into distinct cellular states or fates, underscoring the heterogeneity within mesenchymal cells that guides their progression through various developmental paths ( Figure 11A ). The mesenchymal cell trajectory was divided into nine stages, with cells evenly distributed across early, middle, and late developmental stages, highlighting the continuous and complex nature of mesenchymal cell development ( Figure 11B ).

Notably, under BC conditions, the number of mesenchymal cells across all developmental stages significantly increased, corroborating previous findings ( Figures 11C, D ). This suggests that the BC microenvironment promotes an expansion of mesenchymal cells, enhancing their differentiation potential and possibly contributing to tumor progression and the remodeling of the tumor microenvironment. The study also depicted the distribution of the three mesenchymal subtypes—MSC1, MSC2, and MSC3—along the mesenchymal cell differentiation trajectory ( Figure 11E ). These subtypes were present throughout all stages, indicating the intrinsic heterogeneity and complex interactions among mesenchymal cells within the tissue. Interestingly, MSC1 was predominantly more abundant at every stage, possibly reflecting its central role in mesenchymal cell development and function. This suggests that MSC1 may possess greater biological activity or functionality, influencing both neighboring cells and the broader cellular community.

Previous analyses revealed significantly lower expression levels of JUN, NFKB1, and SP1 in BC samples from both the whole transcriptome RNA-seq dataset and the TCGA-BC dataset (P < 0.01) ( Figures 3E, F ). To further validate these findings at the clinical level, RT-qPCR was employed to assess the expression of these biomarkers in patients with BC ( Figures 12A–C ). RT-qPCR results confirmed that both NFKB1 and SP1 exhibited significantly reduced expression in BC samples (P < 0.05), aligning with initial observations. Although JUN demonstrated a trend toward decreased expression in BC, the difference was not statistically significant (P = 0.0139).