The BC related single-cell dataset GSE114724 was downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) through the GEOquery package. The data source is Homo sapiens, and its data platform is GPL16791. A total of three BC samples were selected for follow-up analysis: GSM3148575, GSM3148577 and GSM3148578. Sample GSM3148575 included 7,096 cells, sample GSM3148577 included 4,926 cells and sample GSM3148577 included 5,232 cells.

In the present study, R packet ‘Seurat’ (18) was used to create single-cell data as Seurat objects, and the proportion of mitochondria in each cell was calculated by the ‘PercentageFeatureSet’ function of Seurat package. Low-quality cells were excluded based on three quality measures: Mitochondrial gene content >5%, feature <500, and unique molecular identifier (UMI) count >20,000. After the aforementioned steps, 10,261 cells were obtained.

Subsequently, the samples were integrated and the Canonical Correlation Analysis (CCA) method was used to remove the batch effect. For Seurat objects, linear dimensionality reduction was performed and the most variable genes expressed in the dataset were used to calculate principal component (PC). Then, the ‘FindNeighbors’ and ‘FindClusters’ functions of Seurat were used to group the cells into the optimal number of cluster in order to identify the cell types, and Uniform Manifold Approximation and Projection (UMAP) was used to reduce the information captured in the selected important PCs to two dimensions, and the visual clustering of cells was realized based on graphics.

BC transcriptome data was downloaded from The Cancer Genome Atlas Genomic Data Commons (TCGA GDC) official website (https://portal.gdc.cancer.gov/), and the selected data type was Fragments Per Kilobase Million. And the clinical data of patients with BC were downloaded from the TCGA GDC website (n=1,285). After removing samples that lacked clinical information, TCGA-BC included 878 BC samples and 41 paracancerous tissue samples. Details of patients with BC are listed in Table I.

In addition, the BC dataset GSE42568 was downloaded from the GEO database through the GEOquery package to verify the accuracy of the prognostic model. The GSE42568(19) data set comes from Homo sapiens, and its data platform is GPL570. It contains a total of 121 samples, including 104 cases of BC and 17 cases of normal samples. BC samples were selected from these samples for survival analysis.

In the present study, NRGs were obtained from previous literature for follow-up analysis, including 101 NRGs (Table SI) (20). At the same time, the immune checkpoint genes (ICGs) were obtained from the previous literature for follow-up analysis; a total of 47 ICGs were obtained (Table SII) (21). To analyze somatic mutations in TCGA patients with BC, ‘TCGAbiolinks’ R package was used to download the patient's ‘Masked Copy Number Segment’ data (n=897).

A total of 14 clusters were visualized by UMAP, and 7 different cell types were revealed by artificial annotation of cell type marker genes, including B cells, CD4⁺ T cell, CD8⁺ T cell, progenitor cells, endothelial cells, natural killer (NK) T cells and regulatory T cells. Cell group marker genes were displayed in Table II. In order to understand the expression pattern of diagnostic markers in patients with BC, the differences of gene expression were compared in different cell groups.

For the seven cell groups that had been annotated, the function ‘FindAllMarkers’ was used to calculate the differential genes among all cell groups, and the genes selected according to the criteria of |log2FoldChange|>0.25 and P<0.05 were used as the single-cell differentially expressed genes (scDEGs) for further study. Then, NRGs and scDEGs were intersected, and the obtained genes were used as the key genes in the present study.

AUCell could identify cells with active gene sets in single-cell RNA-seq data. AUCell used ‘Area Under Curve’ (AUC) to calculate whether the key subset of the input gene set was enriched in the expressed genes of each cell. Key genes for AUCell scoring were selected and high-scoring cell populations were searched.

The differential genes in CD4⁺ T cell population were analyzed by GO annotation and KEGG pathway enrichment analysis using ‘clusterProfiler’ package of R, and the critical value of FDR <0.05 was considered to be statistically significant. The first 8 results with the lowest P-value of biological processes (BP), cellular components (CC) and molecular functions (MF) in GO and the top 10 results with the lowest P-value of KEGG were shown in bar chart and bubble chart, respectively.

Pseudotime analysis could arrange each cell on the corresponding track according to the time sequence gene expression, and divide the sample into several differentiated cell groups according to the gene expression status to generate an intuitive tree map of pedigree development, which could predict the differentiation and development trajectory of cells. For the subsets of CD4⁺ T cells, pseudotime analysis was used to predict the differentiation and development of subsets, and the changes of NRGs during the pseudotime process were analyzed.

Consistent clustering is a method that could determine the number and members of possible clusters in a dataset. The R-packet ‘ConsensusClusterPlus’ was used to cluster the TCGA-BC dataset consistently using the aforementioned key genes in order to improve distinguishing different subtypes of BC. In addition, PC analysis (PCA) was used to analyze different BC molecular subtypes.

In order to determine the effect of key genes on the prognosis of patients with BC, the Least Absolute Shrinkage and Selection Operator (LASSO) algorithm was used to construct the prognostic risk model of BC. The ‘glmnet’ package of R was used to select the LASSO feature genes. The selected features were screened in the process of model construction, and the genes in the best model were selected as the BC prognostic genes. Finally, according to the risk regression coefficients of all genes in the model, a formula for calculating the risk score was constructed:

In the present study, a prognostic risk model for BC samples was constructed from TCGA datasets, and R-packet ‘maxstat’ was used to calculate the best cut-off value of the ability to predict the survival time of patients with BC. Based on cut-off value, the patients were divided into high- and low-risk groups. The Kaplan-Meier (KM) method of R-packet ‘survival’ was used for survival analysis, and R-packet ‘survminer’ was used to visualize the results. In addition, GEO datasets were used to validate the prognostic model.

Based on the TCGA-BC dataset, univariate and multivariate COX regression analysis was used to evaluate the ability of risk scores combined with clinicopathological features of patients with BC to predict OS, and the results were shown through forest maps. In addition, according to the clinical information of the sample in the present study, R-packet ‘rms’ was used to construct a line chart to study the relationship between clinical factors and prognosis, and the calibration curves of 1-, 3- and 5-year, respectively, were drawn.

CIBERSORT (22) calculation method was used to analyze the data of high- and low-risk groups in TCGA-BC dataset to obtain immune infiltration information, and R ‘ggplot2’ package was used to draw a bar chart to show the distribution of 22 kinds of immune cell infiltration in each sample. Subsequently, a heat map was drawn to reflect the correlation between prognostic genes and immune cells. Finally, the scores of different immune cells were compared between high- and low-risk groups, so as to obtain immune cells with different infiltration levels between the two groups.

In the study of Ru et al (23), the data of 28 kinds of immune cells were collected, R-packet Gene Set Variation Analysis was used to explore the difference of immune infiltration among different BC molecular subtypes, and the results were shown by box chart.

In the present study, Estimate R package was used to estimate the proportion of immune components and stromal components of each sample in TME in the form of ImmuneScore, StromalScore and ESTIMATEScore. The higher the score, the greater the proportion of corresponding components in TME.

R-packet ‘maftools’ was used to visualize somatic mutations in high- and low-risk groups in TCGA-BC datasets, and the mutation differences between the two groups were compared; the results were shown in the waterfall map. At the same time, the differential expression of human leukocyte antigen (HLA) family genes was analyzed in high- and low-risk groups.

The MDA-MB-231 cell line (CL-0150, Procell Life Science & Technology Co., Ltd.) was cultured at 37˚C with 5% CO₂ in DMEM (cat. no. PM150210; Procell Life Science & Technology Co., Ltd.) with 10% fetal bovine serum (FBS) (cat. no. 10091148; Gibco; Thermo Fisher Scientific, Inc.), 100 µg/ml streptomycin and 100 U/ml penicillin (cat. no. SV30010; Hyclone; Cytiva). The MCF10A cell line (cat. no. CRL-10317; Shanghai Zhong Qiao Xin Zhou Biotechnology Co., Ltd.) was cultured in 37˚C, 5% CO₂, complete medium (cat. no. ZQ-1311; Shanghai Zhong Qiao Xin Zhou Biotechnology Co., Ltd.) containing 5% horse serum and 1% penicillin/streptomycin (Hyclone; Cytiva). MCF-7 cells were borrowed from Key Laboratory of Fertility maintenance, Ministry of Education (Professor Jing Chen). However, appropriate measures were received to ensure that the source and quality of these cells were reliable, and these cells were used for the present research to explore some of the characteristics of BC.

Total RNA was isolated using TRIZOL (Ambion; Thermo Fisher Scientific, Inc.). Afterwards, the RNA underwent phenol-chloroform extractions to further purify it. The quantity and quality (the acceptable ratio of A260/A280 was ≤1.8 and ≥2.2) of the purified RNA were assessed by measuring the absorbance at 260/280 nm (A₂₆₀/A₂₈₀) using Nanodrop One (NanoDrop Technologies; Thermo Fisher Scientific, Inc.). The level of gene expression was assessed using GAPDH as a control gene. cDNA synthesis [HiScript^® III RT SuperMix for qPCR (+gDNA wiper), Vazyme Biotech Co., Ltd.] was performed using standard procedures (37˚C for 15 min, 85˚C for 5 sec, and then maintained at 4˚C) and RT-qPCR analysis was carried out on the Bio-Rad S1000 instrument with HieffÔ qPCR SYBR^® Green Master Mix (Low Rox Plus) (Shanghai Yeasen Biotechnology Co., Ltd.). The procedure was as follows: 95˚C for 5 min, then 40 cycles at 95˚C for 10 sec, and 30 sec at 60˚C. Afterwards, the concentration of each transcript was normalized to the level of GAPDH mRNA using the 2^-ΔΔCq method (24). The primer sequences were as follows: GAPDH forward, 5'GGTCGGAGTCAACGGATTTG-3' and reverse, 5'-GGAAGATGGTGATGGGATTTC-3'; BCL2 forward, 5'-TCATGTGTGTGGAGAGCGTCAAC-3' and reverse, 5'-GTGTGCAGGTGCCGGTTCAG-3'; BIRC3 forward, 5'-TATCCACATCAGACAGCCCAGGAG-3' and reverse, 5'-TTCCACGGCAGCATTAATCACAGG-3'; AIFM1 forward, 5'-GGCGGCGGGTGCTTTGAAG-3' and reverse, 5'-CATGCCATCGCTGGAACAAGTTG-3'; IFNG forward, 5'-TGACTTGAATGTCCAACGCAAAGC-3' and reverse, 5'-CGACCTCGAAACAGCATCTGACTC-3'; and VDAC1 forward, 5'-GATTGACCCTGACGCCTGCTTC-3' and reverse, 5'-CTTGCCATCCAGAAGAGCTGACAG-3'.

All data calculation and statistical analysis were carried out by R software (https://www.r-project.org/; 4.1.2). For the comparison of two groups of continuous variables, the statistical significance of normal distribution variables was estimated by independent Student t-test, and the difference between non-normal distribution variables was analyzed by Wilcoxon rank sum test. For the comparison of several groups of continuous variables, ANOVA (one-way, parametric) was used when the variance was uniform, and Kruskal-Wallis test (non-parametric) was used when the variance was uneven, and then Dunnett's was used for the appropriate multiple comparison test after the occurrence of significant results. The survival package of R was used for survival analysis, the Kaplan-Meier survival curve was used to show the difference in survival, and the log-rank test was used to evaluate the significant difference in survival time between the two groups. Univariate and multivariate Cox regression analysis was based on survival package, and LASSO model was based on ‘glmnetR’ package. All the statistical P-values were bilateral, and P<0.05 was considered to indicate a statistically significant difference.

The flow chart of the overall analysis of the present study is demonstrated in Fig. 1. The BC samples of single-cell data were integrated and CCA method was used to remove the batch effect, and then ‘Seurat’ package was used to analyze the cellular heterogeneity of single-cell data. After removing the cells with mitochondrial gene content >5%, feature number <500 and UMI >20,000, the cluster diagrams of 7 cell types in BC samples were obtained by UMAP clustering (Fig. 2A). Then the expression of marker genes in different cell clusters of 7 single-cell types was analyzed and shown by a bubble diagram (Fig. 2B). It was found that there were significant differences in marker genes' expression among different cell clusters (Fig. 2C-K). At the same time, the expression differences among different cell populations were analyzed and the expression heat map of top 10 differential genes was constructed (Fig. 2L). Next, the differences in the proportion of different cell groups in patients with BC were also analyzed. It was found that CD8+ T cells accounted for the highest proportion in BC samples, followed by B cells and regulatory T cells, while progenitor cells accounted for the lowest, and other cell types had a similar proportion (Fig. 2M).

For the seven cell types that had been annotated, the function ‘FindAllMarkers’ was used to calculate the differential genes among cell types, and a total of 1,832 scDEGs were selected based on |log2FoldChange|>0.25 and P<0.05. A total of 101 NRGs were obtained in the literature (18). Then, the intersection of NRGs and scDEGs was received (Fig. 3A) and 15 key genes resulted for follow-up research (BAX, BCL2, PPIA, IFNG, AIFM1, IFNAR2, VDAC1, HMGB1, CHMP5, PYCARD, VDAC3, CHMP7, BIRC3, EIF2AK2 and PARP1).

Through Pearson correlation analysis, it was found that there was a positive correlation in expression of most key genes, which indicated that there might be similar expression patterns of key genes in patients with BC (Fig. 3B). Among them, the correlation between CHMP7 and BIRC3 was the highest, while BIRC3 and PYCARD had the largest negative correlation. Finally, a heat map was used to demonstrate the expression of key genes in the cell populations (Fig. 3C).

In order to verify the activity of key genes in different cell types, AUCell packages were used to demonstrate the activity of key genes in BC samples of single-cell data sets with UMAP maps (Fig. 4A). Through the UMAP map, it was found that CD4⁺ T cells had higher AUCell scores, therefore CD4⁺ T cells were defined as high-scoring cells with NRGs.

In order to explore whether there were significant differences in genes between high-scoring cell population and other cell populations, CD4⁺ T cells were selected as a reference, using the function ‘FindAllMarkers’, and 499 DEGs were screened for GO and KEGG enrichment according to the standard of |log2FC|>0.5 and P<0.01. GO analysis revealed that DEGs were mainly related to T cell activation and regulation of T cell activation (Fig. 4B). KEGG analysis demonstrated that DEGs were related to T cell receptor signaling pathway, Th17 cell differentiation and PD-L1 expression and PD-1 checkpoint pathways in cancer (Fig. 4C). These results suggested that NRGs were highly active in CD4⁺ T cells and might be involved in the activation of CD4⁺ T cells.

In order to explore new subtypes of CD4⁺ T cells, CD4⁺ T cell clusters were clustered and four CD4⁺ T cell subtypes were obtained, which were then annotated with significantly upregulated DEGs, resulting in ANXA1 CD4⁺ T, CCL4 CD4⁺ T, FASBP5 CD4⁺ T and TCF7 CD4⁺ T cells (Fig. 5A). For these four CD4⁺ T cell subtypes, the function ‘FindAllMarkers’ was used to calculate the marker genes between them, and a violin map was employed to show the expression of marker genes in each subtype (Fig. 5B). To verify whether there was a difference in the expression of key genes among subsets of CD4⁺ T cells, a heat was used map to demonstrate the results. It was found that there was no difference in the expression of the majority of key genes' subsets; only the expression of BAX2 was lower in CCL4 CD4⁺ T cells (Fig. 5C). Subsequently, pseudotime analysis was used to show the developmental trajectory of CD4⁺ T cells, and the results revealed that the developmental trajectories among the subtypes were CCL4 CD4⁺ T, ANXA1 CD4⁺ T, TCF7 CD4⁺ T and FASBP5 CD4⁺ T cells (Fig. 5D and E). To analyze the changes in the expression of NRGs in the developmental trajectory, the expression of key genes was illustrated by heat maps (Fig. 5F). The expression of PYCARD, IFNG, CHMP5, BAX, HMGB1, VDAC3, EIF2AK2 and other genes was found to be downregulated during the development of CD4⁺ T cell subtypes, while the expression of BCL2, PPIA, PARP1, VDAC1, AIFM1, BIR3 and other genes was upregulated. This suggested that NRGs could regulate the development of CD4⁺ T cell subtypes.

The aforementioned 15 key genes were employed as feature genes to carry out consistent clustering analysis of the samples in the TCGA-BC dataset. The process of selecting the k-value of the cluster is shown in Fig. 6A, and k=3 was selected as the result of clustering (Fig. 6B). The BC samples were divided into three disease subtypes: Subtype 1, 2 and 3. The clustering results were demonstrated by heat maps (Fig. 6C). PCA could distinguish different subtypes to a certain extent, indicating that there was a certain heterogeneity between different subtypes (Fig. 6D).

Single-sample Gene Set Enrichment Analysis was used to evaluate the differences in immune cells between BC molecular subtypes (Fig. 6E). The results revealed that the degree of immune cell infiltration in BC subtypes from high to low was subtypes 3, 1 and 2. Survival differences were then compared among the three molecular subtypes and illustrated with a KM survival curve (Fig. 6F). It was found that the survival rate of subtypes 1 and 3 was higher than that of subtype 2 (P<0.0001), which was consistent with a higher degree of immune cell infiltration and improved prognosis in these subtypes.

The expression of key genes among different subtypes is shown with a box chart. The results identified that BAX, IFNG, IFNAR2, HMGB1, CHMP7, BIRC3, EIF2AK2 and PARP1 were most highly expressed in subtype 1 (Fig. 7A, D, F, H, L, M, N and O, respectively). The expression of PPIA, AIFM1, VDAC1 and VDAC3 was the highest in subtype 2 (Fig. 7C, E, G and K). BCL2, CHMP5 and PYCARD were most highly expressed in subtype 3 (Fig. 7B, I and J). This suggested that different subtypes could be distinguished by the expression level of key genes, and the prognosis of patients could be predicted.

To determine the prognostic genes of BC and analyze their diagnostic ability for diseases, based on 15 key genes (BAX, BCL2, PPIA, IFNG, AIFM1, IFNAR2, VDAC1, HMGB1, CHMP5, PYCARD, VDAC3, CHMP7, BIRC3, EIF2AK2 and PARP1), in the TCGA-BC datasets, univariate Cox regression analysis was first used to screen prognostic genes, and 5 genes were screened according to the criterion of P<0.05 (Fig. 8A). These 5 genes (BCL2, BIRC3, AIFM1, IFNG and VDAC1) in the model were identified as potential prognostic genes of BC by LASSO regression (Fig. 8B). According to the Cox risk regression coefficient of the 5 genes, the riskScores of each sample were calculated. The samples were divided into two groups by the cut-off provided by the R package ‘maxstat’, including 440 samples in the high-risk group and 599 samples in the low-risk group. According to the aforementioned grouping, survival analysis was conducted and a survival curve was drawn. In the LASSO prognostic model, the survival rate of the high-risk group was lower than that of the low-risk group (P<0.0001) (Fig. 8C).

Next, the prognosis model was validated with GEO datasets and a survival curve was generated through survival analysis. In contrast to the aforementioned results, the survival rate in the high-risk group was higher than that in the low-risk group (P<0.042) (Fig. 8D). This contradictory result may be caused by the heterogeneity between datasets and the differences in sample characteristics. Finally, the Receiver Operating Characteristic (ROC) curve was drawn and the AUC was calculated to verify the model. The results revealed that the AUC of TCGA-BC dataset was 0.595 (Fig. S1A), which had a certain prediction effect, while the AUC of GEO datasets was 0.457 (Fig. S1B), indicating that its prediction effect was poor.

Univariate and multivariate Cox regression analysis was performed based on the TCGA-BC dataset according to the patients' TNM stage, sex, age and riskScore. In univariate Cox regression, age, T-stage, N-stage, M-stage and riskScore had a significant influence on OS (P<0.001) (Fig. 9A), while only age, N-stage and riskScore had a significant effect in multivariate Cox regression (P<0.001) (Fig. 9B). A nomogram was also drawn using the ‘rms’ package (Fig. 9C). The 1-, 3- and 5-year survival probabilities of patients with BC were predicted by drawing 1- (Fig. 9D), 3- (Fig. 9E) and 5-year calibration curves (Fig. 9F). The results demonstrated that the 5-year nomogram model was the most consistent with the ideal model, and other nomogram models were essentially equal to the ideal model, indicating that the present model was relatively accurate.

First, the correlation between 5 prognostic genes (BCL2, BIRC3, AIFM1, IFNG and VDAC1) and 22 types of immune cells was analyzed and visualized using CIBERSORT in the BC dataset (Fig. 10A). The positive correlation between IFNG and M1 macrophages was the highest (r=0.65, P=1.6x10^-125), while the negative correlation between IFNG and M2 macrophages was the highest (r=-0.29, P=5.1x10^-22). The infiltration panorama of 22 types of immune cells in each sample was next mapped (Fig. 10B). Finally, the ImmuneScore, StromalScore and EstimateScore in each sample were calculated by the ESTIMATE algorithm, and the correlation between each score and 22 types of immune cells was analyzed and visualized (Fig. 10C). The highest positive correlation was between ImmuneScore and M1 macrophages (r=0.5, P=7.7x10^-66). Finally, the differences in immune infiltration between high- and low-risk groups were shown as box plots (Fig. 10D). The results revealed that the infiltration of plasma cells, activated NK cells, M0 and M2 macrophages and eosinophils in the high-risk group was significantly higher than that in the low-risk group. However, compared with the high-risk group, the low-risk group had a higher degree of infiltration of naive B, CD8⁺ T cells, memory resting CD4⁺ T cells, memory activated CD4⁺ T cells and follicular helper T cells, as well as M1 macrophages and other immune cells.

According to the risk score of the aforementioned model, the total BC samples in the somatic mutation dataset were divided into high- and low-risk group, and the characteristics of various mutations in different groups were analyzed. The top 20 mutated genes in the high- and low-risk groups were analyzed (Fig. 11A and B), and it was found that there were differences in gene mutations in the two groups. The gene with the highest mutation frequency in the high-risk group was TP53, while in the low-risk group it was PIK3CA. The differential expression of HLA family genes in the high- and low-risk groups of TCGA-BC dataset was also analyzed. As revealed in Fig. 11C, various HLA molecules were upregulated in the low-risk group.

The constructed LASSO model produced 5 BC prognostic genes (BCL2, BIRC3, AIFM1, IFNG and VDAC1), and 47 ICGs were obtained by consulting the literature (19). By calculating the Pearson correlation coefficient between the two, correlations of prognostic genes with immune checkpoints were obtained (Fig. 12). Among them, IFNG had the highest positive correlation with IDO1 and ICOS (r=0.79 and P=2.4x10^-232; r=0.79 and P=4.9x10^-237, respectively), while BCL2 had the highest negative correlation with TNFSF9 (r=0.34, P=1.3x10^-30).

The mRNA expression of 5 prognostic genes (BCL2, BIRC3, AIFM1, IFNG and VDAC1) in the human BC cell lines MDA-MB-231 and MCF-7 and the human breast epithelial cell MCF-10A was detected by RT-qPCR. The results showed that, compared with MCF-10A, BCL2, BIRC3, AIFM1, IFNG and VDAC1 were significantly higher in the aforementioned BC cell lines (Fig. 13A-E) (P<0.05). Although the expression of BCL2 and BIRC3 in MCF-7 cells was higher than that in MCF-10A, there was no statistical significance (P>0.05). This suggested that the expression of BCL2 and BIRC3 genes may be related to BC types.