Skip to main content
Download PDF

Prognostic value of circadian rhythm-associated genes in breast cancer

World Journal of Surgical Oncology volume 23, Article number: 186 (2025) Cite this article

Abstract

Objective

Breast cancer (BC) remains the most prevalent malignancy among women. Clinical evidence indicates that genetic variations related to circadian rhythms, as well as the timing of therapeutic interventions, influence the response to radiation therapy and the toxicity of pharmacological treatments in women with BC. This study aimed to identify key circadian rhythm-related genes (CRGs) using bioinformatics and machine learning, and construct a prognostic model to predict clinical outcomes.

Methods

Transcriptome data for BC were retrieved from The Cancer Genome Atlas database. Univariate Cox regression and least absolute shrinkage and selection operator regression analyses were used to develop a prognostic model based on CRGs. The predictive performance of the risk score model was evaluated. Univariate and multivariate Cox regression analyses were applied to construct the prognostic model and stratify patients into high-risk and low-risk groups. Additionally, differences in immune microenvironment, immunotherapy efficacy, and tumor mutation burden were assessed between risk groups.

Results

A prognostic risk score model comprising 17 CRGs was developed. The areas under the receiver operating characteristic curve for overall survival at 1, 3, 5, and 7 years exceeded 0.6, indicating acceptable predictive performance. Calibration plots and decision curve analyses demonstrated the use of the model in prognostic prediction. Significant differences in immune microenvironment, immunotherapy efficacy, and tumor mutation burden were identified between the low-risk and high-risk groups.

Conclusion

The circadian rhythm-based gene model, effectively predicted the prognosis of individuals with BC, highlighting its potential to inform personalized therapeutic strategies and improve patient outcomes.

Introduction

Breast cancer (BC) is among the most frequently diagnosed cancers in women [1]. Improvements in diagnostic techniques and treatment strategies have led to a notable reduction in BC-related mortality [2]. However, challenges persist in accurately predicting individual patient responses to specific treatments or identifying those at risk for adverse effects [3, 4]. Additionally, determining prognostic outcomes and selecting appropriate salvage therapies for patients with advanced BC remain significant clinical hurdles. As a result, there is an ongoing need to identify diagnostic biomarkers and therapeutic targets to develop novel strategies for the management and treatment of BC.

Circadian rhythm refers to biological processes that follow an approximately 24-hour cycle. These rhythms are regulated by a central oscillator located in the suprachiasmatic nucleus of the hypothalamus and several peripheral oscillators in organs such as the liver, heart, and kidneys [5]. The central and peripheral rhythms are synchronized, enabling the body to anticipate and adapt to daily behavioral and environmental changes under normal physiological conditions. Disruptions in circadian rhythms, however, have been associated with an increased risk of various diseases, including asthma, epilepsy, and Parkinson’s disease [6,7,8]. Notably, circadian rhythm disorders are considered a significant risk factor for several types of BC development [9,10,11].

Circadian rhythms have been shown to influence both the efficacy of chemotherapy and the severity of treatment-related side effects in patients with BC [12, 13]. Xue et al. demonstrated that several clock genes, including DEC1, DEC2, PER1, PER2, PER3, CRY2, and RORC, exhibit a negative correlation with proliferation-associated genes involved in cell division. Furthermore, certain clock genes, such as CRY2, PER2, PER3, and DEC1, exhibit strong positive correlations with estrogen receptor (ER) signaling-dependent genes, particularly DEC1 [14]. These findings indicate that circadian rhythms modulate breast tumor microenvironment cellular components involving alterations in cellular metabolism, gene expression, and activation of signaling pathways [15].

The study aimed to identify key circadian rhythm-associated genes (CRGs) in BC by combining bioinformatics and machine learning, and to construct a prognostic model to assess its utility in predicting clinical outcomes.

Materials and methods

Data source

A total of 1,147 CRGs were identified from the Circadian Gene Database [16]. The TCGA-BRCA gene expression dataset, including 1,104 tumor samples and 113 control samples, was obtained from The Cancer Genome Atlas (TCGA, https://www.cancer.gov/tcga) database [17]. Clinical data were available for 1,104 patient samples, while survival data were accessible for 1,082 patient samples. Samples lacking follow-up survival data were excluded from the analysis. Additionally, BC mutation data were also retrieved from the TCGA database.

Identification of differentially expressed genes in BC

The “DESeq2” package (v 1.34.0) was used to identify differentially expressed genes (DEGs) between the BC and control samples [18]. Thresholds of an adjusted p-value < 0.05 and|log2fold change (FC)| > 1 were applied to filter DEGs for further analysis. A volcano plot and heatmap were generated to visualize these genes. Furthermore, differentially expressed circadian rhythm-related genes (DECRGs) were identified by intersecting the DEGs with the circadian rhythm-related gene set.

Functional enrichment analysis of DECRGs

Functional enrichment analyses were performed to explore the functional roles of DECRGs. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were conducted using the “clusterProfiler” R package [19]. GO terms and KEGG pathways with an adjusted p-value < 0.05 were considered statistically significant and results were visualized for interpretation.

Protein–protein interaction network

The protein–protein interaction (PPI) network for the DECRGs was constructed using data from the Search Tool for the Retrieval of Interacting Genes (STRING, https://string-db.org/) database to explore the protein interaction relationships of the candidate genes [20]. A minimum interaction score threshold of > 0.4 was applied, and isolated nodes within the network were excluded.

Identification of prognostic genes and establishment of the risk model

The TCGA-BRCA dataset was randomly divided into a training set (n = 758) and a testing set (n = 324) in a 7:3 ratio. Univariate Cox regression analysis was conducted on the DECRGs in the training set to identify genes significantly associated with patient survival, with a p-value < 0.05 used as the selection criterion for candidate genes. The least absolute shrinkage and selection operator (LASSO) regression analysis was subsequently applied to develop the circadian rhythm-related prognostic gene signature. Ten-fold cross-validation was used to select the optimal lambda value (lambda.min), minimizing the model error rate. Genes whose regression coefficients were not penalized to zero were included as prognostic genes in the final model.

A prognostic risk score model was constructed using the formula: risk score = expression level of Gene1 × β1 + expression level of Gene2 × β2 +… + expression level of Genen × βn where β represents the regression coefficient obtained through LASSO regression analysis.

Subsequently, the expression levels of prognostic genes in tumor tissues and control tissues were subsequently analyzed using the Wilcoxon test. Based on the optimal expression threshold of the prognostic genes, BC samples were divided into high-expression and low-expression groups, and Kaplan-Meier (K-M) curves were generated to assess the relationship between biomarker expression and survival outcomes. In order to further evaluate the effectiveness of the risk model, using the median risk score as a cutoff, patients were classified into high-risk and low-risk groups. Kaplan–Meier (KM) analysis was conducted to compare overall survival (OS) between the two groups. Receiver operating characteristic (ROC) curves for 1-, 3-, 5-, and 7-year OS were generated to assess the predictive performance of the model. The prognostic value of the testing set was evaluated following the same procedure applied to the training set.

Gene set variation analysis (GSVA)

GSVA was conducted to estimate variations in gene set enrichment across the groups within the expression dataset. The “GSVA” package (v 1.46.0) [21] was used to analyze biological functions and pathways associated with the high-risk and low-risk groups. Pathways and functions with an adjusted p-value < 0.05 were deemed statistically significant.

Immune microenvironment analysis in high- and low-risk groups

Immune and stromal scores were calculated to evaluate the proportions of immune and stromal components within the tumor microenvironment, using the ESTIMATE algorithm, available in the “ESTIMATE.” package (v 1.0.13) [22]. Additionally, the CIBERSORT algorithm was applied to assess the proportions of 22 immune cell subtypes in the high-risk and low-risk groups [23], and to analyze the correlations between different immune cells. Subsequently, the Wilcoxon test was used to analyze the differences in immune cell infiltration between the two risk groups, with p-values being adjusted using the FDR method.

The expression of 34 immune checkpoint genes and 23 HLA family genes between the two risk groups was analyzed using the Wilcoxon test (p-value < 0.05), and boxplots were generated to display the top 15 differentially expressed immune checkpoint molecules. Subsequently, the TIDE immune dysfunction and immune rejection scores for each cancer sample were obtained using the online database TIDE (http://tide.dfci.harvard.edu/), and the correlation with risk scores was analyzed to assess the potential response to immunotherapy.

Prediction of therapeutic sensitivity

Additionally, the predictive ability of the risk score to assess responses to immunotherapy and seven drugs: 5-fluorouracil, cisplatin, docetaxel, doxorubicin, paclitaxel, shikonin, and vinorelbine was evaluated. The half-maximal inhibitory concentrations (IC50) of these drugs were estimated using the pRRophetic algorithm, with the values normalized for analysis. Subsequently, the Wilcoxon method was applied to analyze the differences in immunogenicity prediction scores (IPS) between the high-risk and low-risk groups.

Analysis of tumor mutation status in patients with varying risk scores

To further investigate the mutation status of prognostic genes in the training set, the “maftools” package (v 2.10.5) [24] was used to analyze the mutation data of patients in both groups. The top 20 frequently mutated genes were visualized in a waterfall plot. Additionally, the tumor mutation burden (TMB) score for each patient sample in the TCGA dataset was calculated, and the differences in TMB scores between the high-risk and low-risk groups were compared.

Predictive nomogram construction and evaluation

To assess the applicability of the risk model to different clinicopathological features, the Kruskal-Wallis test was used to analyze the distribution of risk scores in BC patients with different clinical characteristics (gender, M, N, T stages). Subsequently, the risk scores were combined with other clinical indicators to further identify independent prognostic factors associated with BC and construct a prognostic model. The risk score, patient age, gender, and tumor TNM stage were included in the risk model, with survival time and survival status as dependent variables. Independent prognostic factors identified through univariate and multivariate Cox regression analyses (p-value < 0.05) using the survival package (v 3.3-1) [25]. Then, a nomogram was constructed using the rms package (v 6.2-0) [26] for predicting overall survival (OS) at 1-, 3-, 5-, and 7-years for patients with BC. Calibration curves were then generated to assess the concordance between predicted and actual survival outcomes. Decision curve analysis (DCA) was conducted to assess the clinical utility of the nomogram.

Tissue samples and cell lines

The MCF10A, SK-BR-3, MCF7, and MDA-MB-231 cell lines were cultured in DMEM supplemented with 10% FBS and 1% penicillin-streptomycin. BC tissue samples were collected from patients with BC who underwent surgical procedures, while healthy breast tissue samples were obtained from patients with benign breast tumors. Approximately 50 µg of each tissue sample was used for western blot analysis, with the remaining pathological tissue processed for routine diagnostic evaluation in the pathology department.

Real-time PCR

Real-time PCR was conducted to determine the target gene copy numbers relative to glyceraldehyde-3-phosphate dehydrogenase (GAPDH) transcripts using specific primers. Gene expression levels were quantified with an SYBR Green kit (Biosharp, China), and data analysis was carried out using MxPro quantitative PCR software (Agilent Technologies). The oligonucleotide primer sequences for GAPDH, ARAGAP26, BACH2, and KCNH8 are provided in Table 1.

Table 1 Primer sequences for GAPDH, ARHGAP26, BACH2, and KCNH8 genes
Full size table

The PCR reactions were conducted in a 20 µl reaction mixture containing 10 µl of 2× QuantiFast SYBR Green PCR Master Mix, 2 µM of each primer, and 2 µl of cDNA. The reaction conditions were as follows: an initial denaturation step at 95 °C for 5 min, followed by 40 cycles of denaturation at 95 °C for 10 s, and a combined annealing and extension step at 60 °C for 30 s. All reagents were provided by Biosharp, China. Each experiment was performed in duplicate to ensure reproducibility and accuracy of the results.

Western blotting

The tissue samples were weighed and homogenized in RIPA lysis buffer containing a protease inhibitor cocktail at a ratio of 1:10. The samples were minced until no visible large particles remained, sonicated using a hand-held ultrasonic device for 10 s, and lysed on ice for 30 min. The lysate was then centrifuged at 12,000 rpm, and the supernatant was collected. Protein concentration was determined using the bicinchoninic acid method and adjusted to 5 µg/µl.

Subsequently, 5× loading buffer was added, and the samples were boiled in water for 10 min. Sodium dodecyl sulfate–polyacrylamide gel electrophoresis was conducted at 100 V for 90 min, with 50 µg of protein loaded per lane. The separated proteins were transferred onto a nitrocellulose membrane at 200 mA for 90 min. The membranes were blocked with 5% skim milk at room temperature for 1 h, followed by overnight incubation at 4 °C with primary antibodies (anti-BACH2, anti-GRAF, and anti-LY6G5C, all diluted 1:1,000; Proteintech).

After washing with Tris-buffered saline with Tween®, the membranes were incubated with secondary antibodies (Goat anti-rabbit IgG (H + L), Proteintech, USA, CAS: SA00001-2) diluted 1:10,000 for 1 h at room temperature. Signals were detected using an ImageQuant™ LAS 4000 Mini Imager. Protein expression levels were quantified using ImageJ software (version 180), with the gray value of target proteins normalized to that of GAPDH.

Immunohistochemistry

The tissue samples were dehydrated, embedded in paraffin, and sectioned into 4-µm thick slices. Immunohistochemical staining was conducted to assess the formation of mammary granulomas following modeling and to analyze the expression levels of BACH2, LY6G5C, and GRAF.

The tissue sections were deparaffinized and treated with 3% hydrogen peroxide for 20 min to inhibit endogenous peroxidase activity. Antigen retrieval was conducted using EDTA at 99 °C for 30 min, followed by natural cooling. The sections were then rinsed three times with PBS, each rinse lasting 3 min, and subsequently blocked with BSA for 30 min to prevent non-specific binding.

Primary antibodies (ER at 1:200 and PR at 1:800) were diluted appropriately and applied to the sections, which were incubated at room temperature for 1 h. Following incubation, the sections were washed three more times with PBS for 3 min per wash. Enhanced enzyme-labeled goat anti-rabbit IgG polymer was then added, and the sections were incubated at room temperature for 20 min before being washed again with PBS to remove any unbound secondary antibodies.

A freshly prepared DAB (3,3’-diaminobenzidine) color-developing solution, mixed at a ratio of 1:19, was applied dropwise to the sections, which were incubated at room temperature for 2–10 min. The sections were rinsed under running tap water and counterstained with Kazurkin stain for 50 s to enhance contrast. Finally, the sections were washed under tap water for an additional 10 min to remove any excess stain.

Statistical analysis

Statistical analysis was conducted using SPSS™ Statistics version 26.0. Continuous data are expressed as mean ± standard deviation (x ± s), while count data are presented as percentages (%). Comparisons between two groups were conducted using the t-test, and comparisons among multiple groups were conducted using one-way analysis of variance. A p-value of < 0.05 was considered indicative of statistically significant differences.

Results

Identification of DECRGs in BC

A total of 7,322 DEGs were identified between BC and control samples based on the criteria of an adjusted p-value < 0.05 and an absolute|log2FC| > 1. Among these, 4,592 genes were upregulated, while 2,730 genes were downregulated (Fig. 1a and Supplementary Table S1). The top 15 upregulated and downregulated genes, ranked by adjusted p-value, were visualized in a heatmap (Fig. 1b). Subsequently, 295 genes overlapping between the DEGs and CRGs were identified as DECRGs for further analysis (Fig. 1c and Supplementary Table S2). Of these, 175 DECRGs were upregulated and 120 were downregulated in tumor tissues (Fig. 1d).

Fig. 1
figure 1

Differential expression analysis. (a) A volcano plot depicting 7,322 DEGs in BC samples compared with normal samples. (b) A heatmap showing the expression levels of the top 12 upregulated and downregulated DEGs in BC samples. (c) A Venn diagram illustrating the overlap between DEGs and circadian rhythm-related genes, identifying DECRGs. (D) A bar chart depicting the upregulated and downregulated DECRGs in BC samples. Abbreviations: BC: breast cancer; DEGs: differentially expressed genes; DECRGs: differentially expressed circadian rhythm-related genes

Full size image

Functional enrichment analysis and PPI network of DECRGs

Subsequently, enrichment analysis was performed on the DECRGs, which identified 99 GO biological process terms and 7 KEGG pathways (Supplementary Table S3). The top 10 GO biological process terms were illustrated in Fig. 2a. The DECRGs were predominantly enriched in pathways associated with protein-DNA complexes, including “protein–DNA complex assembly” and “protein–DNA complex subunit organization,” as well as chromatin-related pathways such as “chromatin remodeling,” “chromatin organization,” “chromatin assembly,” and “chromatin assembly or disassembly.” The 7 KEGG pathways identified are illustrated in Fig. 2b. The DECRGs were primarily enriched in pathways related to “alcoholism,” “systemic lupus erythematosus,” and “neutrophil extracellular trap formation.” Additionally, Fig. 2c presented the PPI network of DECRGs, constructed using the STRING database.

Fig. 2
figure 2

The functional enrichment and PPI network construction of candidate genes. (a) Bubble plot representing the top 10 enriched GO terms for the DECRGs. (b) Bubble plot revealing the top 7 KEGG pathways enriched in DECRGs. (c) PPI network constructed for the DECRGs. Abbreviations: GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; PPI: protein–protein interaction; DECRGs: differentially expressed circadian rhythm-related genes

Full size image

Prognosis prediction model of BC based on DECRGs

Based on this, univariate Cox regression was conducted on the training dataset to identify DECRGs significantly associated with prognosis. This analysis led to the selection of 26 DECRGs with a p-value < 0.05 (Fig. 3a). LASSO regression analysis was subsequently used to refine the selection, identifying 17 signature genes for constructing the risk score model (Fig. 3b, c). The regression coefficients of these signature genes were detailed in Table 2. The expression analysis results of the prognostic genes in the training set are shown in Additional file 1 and Supplementary Table S4. Survival analysis showed that the expression levels of these prognostic genes are correlated with the prognosis of patients (Additional file S2).

Fig. 3
figure 3

The construction and evaluation of the risk model in the training set. (a) Univariate Cox regression analysis of 26 DECRGs to identify prognostic genes. (b) The error plot of LASSO cross-validation. The dashed line on the left represents lambda.min, which is the position of the minimum cross-validation error, and the number of feature genes is shown above it. The dashed line on the right represents the line with fewer features. (c) The LASSO coefficient path plot. The lines in different colors correspond to different genes and their coefficients; as lambda increases, the coefficients of the feature variables approach zero, which is the optimal condition for selecting the target genes. (d-e) Risk score distribution and associated survival status of samples in the training set. (f) KM survival curve comparing overall survival between the high-risk and low-risk groups. (g) ROC curve evaluating the predictive performance of the prognostic model in the training set. Abbreviations: DECRGs: differentially expressed circadian rhythm-related genes; LASSO: least absolute shrinkage and selection operator; KM: Kaplan–Meier; ROC: receiver-operating characteristic

Full size image
Table 2 Coefficients of the signature genes
Full size table

Patients were then stratified into high-risk and low-risk groups based on the median risk score. An analysis of the risk score distribution and corresponding survival status revealed that mortality rates increased with higher risk scores (Fig. 3d, e). Survival analysis demonstrated a significant difference in overall survival between the two groups (p-value < 0.05; Fig. 3f), with the high-risk group exhibiting a notably poorer prognosis compared to the low-risk group.

ROC curve analysis indicated that the AUC values for 1-, 3-, 5-, and 7-year survival predictions were 0.677, 0.733, 0.763, and 0.794, respectively (Fig. 3g). These results suggest that the prognostic model demonstrated moderately effective performance in predicting long-term survival outcomes.

Verification of the prognostic significance of the risk score within the testing cohort

A parallel analysis, following the methodology outlined in Section “Prognosis prediction model of BC based on DECRGs”, was conducted on the testing set to validate the predictive performance of the 17-gene signature. Patients were stratified into high-risk and low-risk groups using the median risk score as the cutoff point (Fig. 4a, b). KM survival curve analysis demonstrated that patients in the high-risk group exhibited significantly poorer survival outcomes compared to those in the low-risk group (p-value < 0.05; Fig. 4c).

Fig. 4
figure 4

The evaluation of the risk model in the validation set. (a-b) Risk score distribution and associated survival status of samples in the testing set. (c) KM survival curve showing significant survival differences between high-risk and low-risk groups. (d) ROC curve assessing the prognostic accuracy of the model in the testing set. Abbreviations: KM: Kaplan Meier; ROC: receiver-operating characteristic

Full size image

Additionally, the time-dependent ROC curve analysis revealed that the AUC values for the risk score signature were 0.692, 0.692, 0.634, and 0.656 for the 1-year, 3-year, 5-year, and 7-year survival periods, respectively (Fig. 4d).

GSVA between high- and low-risk groups

GSVA was conducted on the BC dataset to examine the functional and pathway differences between the high-risk and low-risk groups. The high-risk group revealed enrichment in biological process terms such as “MAINTENANCE_OF_PROTEIN_LOCALIZATION_IN_ENDOPLASMIC_RETICULUM,” “RETROGRADE_VESICLE_MEDIATED_TRANSPORT_GOLGI_TO_ENDOPLASMIC_RETICULUM,” and “PROTEIN_REFOLDING,” whereas the low-risk group was enriched in terms like “RESPONSE_TO_INTERLEUKIN_9,” “RESPONSE_TO_INTERLEUKIN_2,” and “GAMMA_DELTA_T_CELL_DIFFERENTIATION” (Fig. 5a).

Fig. 5
figure 5

GSVA between the high- and low-risk groups. (a) Top five pathways enriched in BP terms between the high- and low-risk groups. (b) Top five pathways enriched in CC terms. (c) Top five pathways enriched in MF terms. (d) KEGG pathways enriched in the high- and low-risk groups, highlighting distinct biological mechanisms in each group. Abbreviations: GSVA: gene set variation analysis; BP: biological process; CC: cellular component; MF: molecular function; KEGG: Kyoto Encyclopedia of Genes and Genomes

Full size image

In cellular component analysis, the high-risk group demonstrated a prominence in “COPI_COATED_VESICLE_MEMBRANE,” “COPI_COATED_VESICLE,” and “VESICLE_COAT,” while the low-risk group exhibited higher levels of “ALPHA_BETA_T_CELL_RECEPTOR_COMPLEX,” “NLRP3_INFLAMMASOME_COMPLEX,” and “INFLAMMASOME_COMPLEX” (Fig. 5b).

For molecular function, enrichment in “CATION_TRANSMEMBRANE_TRANSPORTER_ACTIVITY” and “ION_TRANSMEMBRANE_TRANSPORTER_ACTIVITY” was observed in the high-risk group. Conversely, the low-risk group was enriched in terms such as “DEATH_RECEPTOR_ACTIVITY,” “TRANSMEMBRANE_RECEPTOR_PROTEIN_TYROSINE_KINASE_ADAPTOR_ACTIVITY,” and “PHOSPHATIDYLINOSITOL_3_KINASE_REGULATOR_ACTIVITY” (Fig. 5c).

Regarding KEGG gene sets, the high-risk group was enriched in pathways including “STEROID_BIOSYNTHESIS,” “CARDIAC_MUSCLE_CONTRACTION,” and “VIBRIO_CHOLERAE_INFECTION.” In contrast, the low-risk group demonstrated enrichment in “PRIMARY_IMMUNODEFICIENCY,” “HEMATOPOIETIC_CELL_LINEAGE,” and “T_CELL_RECEPTOR_SIGNALING_PATHWAY” (Fig. 5d).

Differences in immune microenvironment between high and low-risk groups

Afterwards, the immune and stromal scores for BC samples were calculated to assess their correlation with risk scores. The immune score was significantly lower in the high-risk group compared to the low-risk group, whereas the stromal score revealed no significant correlation with risk scores (Fig. 6a).

Fig. 6
figure 6

The immune microenvironment analysis of the high-risk and low-risk groups was conducted. (a) Correlation between immune and stromal scores with the risk scores, highlighting associations between the tumor microenvironment and prognostic risk. (b-c) Abundance levels of 22 immune cell types in samples from the high-risk and low-risk groups, illustrating variations in immune cell infiltration. (d) Relationships among the 22 immune cell types, showing their interactions and relevance within the tumor microenvironment. (e) Differential analysis identifying 14 immune cell types that significantly vary between the high-risk and low-risk groups

Full size image

Furthermore, the CIBERSORT algorithm was used to assess the infiltration of 22 immune cell types in samples from the high-risk and low-risk groups (Fig. 6b-d). Significant differences in the infiltration levels of 14 immune cell types were observed between the two groups (Fig. 6e). Among them, 7 cell types, including B_cells_naive, T_cells_CD4_memory_activated, T_cells_follicular_helper, and T_cells_CD8, were expressed at higher levels in the low-risk group, while 7 cell types, including Macrophages_M0, Macrophages_M2, T_cells_CD4_memory_resting, and Neutrophils, were expressed at higher levels in the high-risk group.

Prediction of therapeutic sensitivity in patients with different risk scores

Further analysis related to immunotherapy showed that, the expression levels of 29 immune checkpoint genes and 23 HLA family genes were found to differ significantly between the high-risk and low-risk groups (Table 3). The top 15 immune checkpoint genes and the top 15 HLA family genes are depicted in Fig. 7a, b. In terms of immunotherapy, patients in the low-risk group exhibited a higher immune dysfunction score (Fig. 7c), a lower immune exclusion score (Fig. 7d).

Table 3 Different HLA family genes
Full size table
Fig. 7
figure 7

Differences in response to immunotherapy between the high- and low-risk groups: (a) Expression levels of the top 15 immune checkpoint genes in the two risk score groups. (b) Expression levels of the top 15 HLA family genes in the two risk score groups. (c) Correlation between the TIDE dysfunction score and the risk score. (d) Correlation between the TIDE exclusion score and the risk score. (e) Drugs exhibiting differential sensitivity between the high-risk and low-risk groups. (f) Differences in the IPS between the high-risk and low-risk groups, reflecting varying responses to immunotherapy. Abbreviations: HLA: human leukocyte antigen; TIDE: tumor immune dysfunction and exclusion; IPS: immunophenoscore

Full size image

To assess the potential of the risk score as a biomarker for drug response, the IC50 of 5-fluorouracil, cisplatin, docetaxel, doxorubicin, paclitaxel, shikonin, and vinorelbine were estimated for the TCGA-BRCA group. The results showed that compared to patients in the high-risk group, those in the low-risk group demonstrated higher sensitivity to these drugs (Fig. 7e), and had higher IPS (Fig. 7f) than those in the high-risk group.

Tumor mutation profiles in high-risk and low-risk groups

To examine the mechanisms underlying the risk stratification in BC, somatic mutations across all genes were analyzed using TCGA data. The waterfall plot identified TP53, PIK3CA, and TTN as the three most frequently mutated genes in BC samples (Fig. 8a). Similarly, TP53, PIK3CA, and TTN were the most commonly mutated genes in the high-risk group (Fig. 8b), whereas PIK3CA, TP53, and CDH1 were the top three mutated genes in the low-risk group (Fig. 8c). The high-risk group exhibited a higher overall rate of somatic mutations compared to the low-risk group (Fig. 8d). Additionally, 19 genes showed significantly different mutation rates between the two groups (Fig. 8e).

Fig. 8
figure 8

Analysis of tumor mutation status in the high- and low-risk groups (a) Waterfall plot displaying the mutation landscape of genes in the BC samples. (b-c) Waterfall plots of mutated genes in the high-risk and low-risk groups, respectively. (d) Comparison of TMB between the high-risk and low-risk groups. (e) Mutation frequency comparison for 19 genes between the two risk score groups. (f) Mutation analysis of the 17 signature genes in breast cancer samples, identifying key genetic changes. Abbreviations: BC: breast cancer; TMB: tumor mutational burden

Full size image

The mutation profiles of the 17 signature genes in the model were analyzed, identifying BACH2, KCNH8, and ARHGAP26 as the most frequently mutated genes in BC samples. The predominant mutation type observed for these genes was missense mutations (Fig. 8f).

Development and evaluation of a predictive nomogram

To further explore the relationship between risk scores, clinical indicators, and prognosis, differences in risk scores between various clinical subgroups were first analyzed. The results showed that patients with higher risk scores demonstrated more advanced T and N stages compared to those with lower risk scores. No significant relationship was found between risk scores and M stage or gender (Fig. 9a). Univariate and multivariate Cox regression analyses identified the risk score, patient age, and TNM stage as significant prognostic factors (Fig. 9b, c).

Fig. 9
figure 9

Independent prognostic analysis of clinical information. (a) The correlation analysis between the risk score and clinical information. (b-c) Identification of significant predictive factors through univariate and multivariate Cox regression analyses. These factors were selected based on their association with prognosis in breast cancer patients

Full size image

Based on these findings, a nomogram was developed to estimate 1-year, 3-year, 5-year, and 7-year survival probabilities for each patient (Fig. 10a). Calibration plots demonstrated a high degree of concordance between the survival probabilities predicted by the nomogram and the actual observed outcomes (Fig. 10b). Furthermore, decision curve analysis (DCA) demonstrated a substantial net benefit of the nomogram in predicting patient survival (Fig. 10c).

Fig. 10
figure 10

Development and validation of a nomogram. (a) Construction of a nomogram that integrates the identified predictive factors to estimate 1-, 3-, 5-, and 7-year survival probabilities for BC patients. (b) Calibration curve demonstrating the concordance between the nomogram-predicted survival probabilities and actual outcomes, indicating the accuracy of the model. (c) DCA evaluating the clinical utility of the nomogram by assessing the net benefit across different threshold probabilities

Full size image

Differential expression of CRGs in BC cell lines and tissues

To assess the expression of CRGs in BC, qRT-PCR was used to assess the levels of BACH2, KCNH8, and ARHGAP25 across various cell lines. The results indicated that the expression of the three genes was higher in BC cell lines (MCF7, MDA-MB-231, and SK-BR-3) compared to the normal breast cell line MCF10A (Fig. 11a).

Fig. 11
figure 11

Validation of the Expression of Circadian Rhythm-Related Genes in BC Cell Lines and Tissues: (a) qRT-PCR Analysis: Expression levels of BACH2, KCNH8, and ARHGAP25 in various cell lines, highlighting differential expression in BC cell lines compared to normal breast cell lines. (b) Western Blotting Results: Protein expression levels of GRAF, LY6G5C, and BACH2 in BC tissues compared to normal tissues, demonstrating significant differences in expression. (c) Immunohistochemistry Results: Visualization of GRAF, LY6G5C, and BACH2 expression in BC and normal tissues, providing further validation of their differential expression. BC: breast cancer

Full size image

Additionally, western blotting and immunohistochemical assays were conducted to examine the effects of these genes on their corresponding proteins and the immune environment. The expressions of GRAF, LY6G5C, and BACH2 were analyzed in BC and healthy breast tissues. Western blotting and immunohistochemistry revealed that GRAF expression was lower in BC tissue compared to healthy breast tissue, whereas LY6G5C and BACH2 exhibited higher expression in BC tissue when compared to healthy breast tissue (Fig. 11b-c).

Discussion

This study aimed to identify DCERGs and explore their potential roles in BC. Previous research has indicated that circadian rhythm-related genes, such as DEC1, DEC2, PER1, PER2, PER3, CRY2, and others, influence various aspects of BC, including cellular proliferation, apoptosis, and the immune microenvironment [14]. The results of this study provide robust support for the involvement of circadian genes in BC, as indicated by bioinformatic analyses and functional annotations.

The differential expression of DCERGs in BC underscores their crucial role in the pathophysiology of the disease. GO enrichment analysis indicated that the DECRGs were primarily associated with processes such as the “formation of protein-DNA complexes,” “organization of protein-DNA complex subunits,” and “modification of chromatin structure.” KEGG enrichment analysis further indicated that the candidate genes were predominantly enriched in pathways related to “Alcohol Dependence,” “Systemic Lupus Erythematosus,” and “Neutrophil Extracellular Trap (NET) Formation.”

NETs are structures composed of chromatin fibers intertwined with histones, granule-derived antimicrobial peptides, and various enzymes. They are released by neutrophils in response to pathogenic or antigenic stimuli and play a catalytic role in tumor progression and metastasis [27]. NETs facilitate the breakdown of the extracellular matrix through proteases, such as matrix metalloproteinase 9 and neutrophil elastase, and they contribute to angiogenesis by releasing vascular endothelial growth factor [28]. Additionally, NETs promote the epithelial-to-mesenchymal transition of cancer cells, thereby enhancing tumor invasion and angiogenesis [29]. By capturing tumor cells and adhering to the vessel wall through the von Willebrand factor, NETs disrupt normal vascular endothelial cell connections, making it easier for tumor cells to penetrate the vessel wall and enter the circulatory system.

In this study, 17 characteristic genes were identified as influential in BC progression through univariate Cox and LASSO regression analyses. A risk model based on these genes was developed to stratify patients with BC into high-risk and low-risk groups, followed by survival analysis. This risk model has the potential to significantly contribute to the understanding and management of BC. These key genes play key roles in BC regulation, either by modulating the tumor immune microenvironment or by preserving the normal epithelial phenotype.

Among these genes, OTOR is notably overexpressed in BC tissue compared to normal samples [30]. In primary ER + BC, elevated expression levels of the triadic ENTPD1/NT5E/ADORA3 signature have been associated with diminished OS, reduced distant metastasis-free survival, and lower progression-free survival rates [31].

In hypoxic conditions, the expression of PADI1 and PADI3 is upregulated, leading to the citrullination of PKM2 and subsequently increasing glycolysis in BC cells [32]. However, this process bypasses typical physiological regulatory mechanisms due to the low availability of serine, resulting in excessive glycolytic activity. This dysregulated metabolic pathway ultimately reduces the proliferation rate of cancer cells.

SUSD3 knockdown impairs actin-rich protrusion formation, reduces stress fibers, and inhibits the formation of large basal adherent spots. This also leads to increased cortical actin levels. This is accompanied by decreased activity of Rho and adherent spot kinases. Cells deficient in SUSD3 exhibit impaired cell spreading, diminished intercellular adhesion, and reduced motility [33]. SUSD3 functions as a novel promoter of estrogen-dependent cell proliferation in BC and plays a key role in regulating cell–cell and cell–matrix interactions, as well as cellular migration.

FOXJ1, a member of the forkhead box (FOX) family of transcription factors, undergoes changes in its subcellular localization during the cell cycle transition in HepG2 cells [34]. ECRG4 acts as a tumor suppressor, and its overexpression inhibits the proliferation and migration of BC cells, causing cell cycle arrest at the G0/G1 phase. Additionally, ECRG4 promotes apoptosis in BC cells by inducing the formation of apoptotic vesicles containing the Cytc/Apaf-1/caspase-9 complex [35].

BCL2A1, a member of the Bcl-2-related protein family, serves as an anti-apoptotic protein by inhibiting cytochrome c release and preventing protease activation [36]. Its expression is significantly higher in breast cancer tissues compared to normal breast tissues and increases with the malignancy of various molecular subtypes of BC.

Expression of LEF1 in invasive micropapillary BC is associated with increased lymphovascular invasion, lymph node metastasis, and disease recurrence. Furthermore, β-catenin, a key mediator of the Wnt signaling pathway, accumulates in the cytoplasm and translocates to the nucleus, where it forms a transcription-regulating complex with LEF1 [37].

Although not yet reported in BC, GOLGA8B has been found to be highly expressed in metastatic tissues of lung squamous carcinoma, where it shows a significant correlation with distant metastasis-free survival in patients. Loss-of-function studies indicate that silencing GOLGA8B inhibits tumorigenesis in vivo and reduces the invasion and migration capabilities of cancer cells in vitro [38].

Mutations in cytokine receptors, such as Flt-3, impair apoptotic mechanisms and affect the Ras/Raf/MEK/ERK signaling pathway, thereby influencing tumor progression [39].

STX11, a protein localized on endosomal membranes and lysosomes in macrophages, is involved in vesicle transport and secretion [40]. Silencing STX11 has been shown to enhance the phagocytosis of apoptotic and antibody-dependent target cells, while also increasing TNF-α secretion, contributing to antitumor activity [41].

EIF4A3/DDX48, an ATP-dependent RNA decapping enzyme of the eIF4A family, exhibits reduced expression of eIF4E3 compared to eIF4E1 in BC cells [42]. Inhibition of EIF4A3 activates p53, thereby promoting tumor suppression [43].

BACH2 regulates CD8+ and CD4+ effector T cells by limiting their effector cytokine expression within tumors [44]. Additionally, elevated expression of circBACH2 has been shown to enhance cell proliferation, invasion, and migration in triple-negative BC cells [45].

The PASK gene regulates the expression of PAS kinase, which interacts with the mTOR pathway by phosphorylating S6, a key target of mTOR [46]. The downstream mTOR pathway represents a critical target for cancer therapy.

IFNG, which encodes IFN-γ, plays a key role in coordinating innate immune and inflammatory responses. IFN-γ secreted by M1 macrophages, CD8+ T cells, and dendritic cells inhibits tumor cell proliferation, thereby contributing to antitumor immunity [47].

ARHGAP26 (also known as GRAF1) is located on human chromosome 12q24.31 and is a member of the RhoGAP family. It negatively regulates the activity of RhoA and Rac1 via its RhoGAP domain, thereby participating in tumorigenesis and tumor progression [48]. Additionally, GRAF1 is the protein encoded by ARHGAP26. Our bioinformatics analysis revealed that ARHGAP26 undergoes genetic mutations, which lead to downregulation of GRAF1 expression [49]. In this study, we found that compared with normal breast cells, ARHGAP26 is highly expressed in breast cancer cells. In contrast, GRAF1 protein is expressed at lower levels in human breast cancer tissue samples than in adjacent non-cancerous tissues, indicating that mutated ARHGAP26 regulates GRAF1 expression. Thus, ARHGAP26 and GRAF1 are closely related.

KCNH8 is a member of the voltage-gated potassium channel family (Kv12.3) and is closely associated with the electrophysiological activities of tumors [50]. LY6G5C belongs to the LY6/uPAR protein family, which typically encodes cell surface antigens or signaling molecules and may play a role in the tumor microenvironment [51]. Studies have shown that co-immunoprecipitation (Co-IP) experiments reveal a weak interaction between the LY6 domain of LY6G5C and the C-terminal PDZ motif (SVKI) of KCNH8 (Kd = 8.7 µM) [52]. Therefore, we hypothesize that LY6G5C is a downstream protein of KCNH8. This hypothesis was validated in cellular and tissue samples, revealing that compared with normal breast cells, KCNH8 is highly expressed in breast cancer cells. Meanwhile, LY6G5C protein is expressed at higher levels in human breast cancer tissue samples than in adjacent non-cancerous tissues, suggesting a certain correlation between KCNH8 and LY6G5C.

GSVA indicated distinct biological activities between the high-risk and low-risk groups. The high-risk group exhibited significant involvement in processes related to endoplasmic reticulum stress, such as the maintenance of protein localization within the ER. In contrast, the low-risk group was primarily engaged in immune microenvironment-related processes, including the formation of T-cell receptor complexes. Genetic changes within cancer cells can exacerbate ER stress and sustain the activation of the unfolded protein response pathway [53]. For instance, the deletion of tumor suppressor genes and the hyperactivation of oncogenes often lead to increased protein synthesis to meet the heightened metabolic demands of tumor development. Furthermore, rapidly proliferating cancer cells require extensive ER expansion to facilitate cell division and distribute ER content to daughter cells [54].

Cancer cells experiencing ER stress secrete factors that recruit or modulate the function of immune cells in the tumor microenvironment, including myeloid cells and natural killer (NK) cells [55, 56]. Additionally, CD3+, CD4+, and CD8+ T cells are key predictors of BC prognosis, as changes in T-cell activity influence the PD-1/PD-L1 pathway, which can either promote or inhibit tumor immune escape [57].

Further investigation into the potential mechanisms of the identified signature genes regulating BC indicated that risk scores were not associated with the M stage or gender but were significantly higher with advancing T and N stages. Moreover, pronounced differences in immune scores were observed between high-risk and low-risk groups. The composition of infiltrating immune cells, including dendritic cells, eosinophils, macrophages, neutrophils, plasma cells, and CD4+/CD8+ T cells, varied significantly between the two groups. These findings indicate that changes in the immune microenvironment contribute to differences in patient prognoses.

Neutrophils, for instance, can inhibit lymphocyte and NK cell cytolytic activity, suppress apoptosis of cancer cells, and enhance tumor cell adhesion to the extracellular matrix. Conversely, lymphocytes augment anticancer immune responses. Blood neutrophils secrete soluble factors that promote interactions between circulating tumor cells and endothelial cells, thereby facilitating metastasis [58, 59]. As a result, circadian-rhythm regulated genes with differential expression influence BC prognosis by modulating NETs and T-cell receptor complexes.

Disparities in immunotherapy responses between the high-risk and low-risk groups were further evaluated. Significant differences were observed in the expression levels of immune checkpoint molecules, such as CD40L, CD244, and BTLA, as well as HLA family genes, including HLA-DPB2, HLA-DQB1, and HLA-E. These findings indicate that patients in different risk groups may require distinct types of immune checkpoint inhibitors.

CD40 stimulation on antigen-presenting cells facilitates the direct activation of CD8+ cytotoxic T cells, bypassing the need for CD4+ T cells. CD244, an NK cell receptor, regulates NK cell cytotoxicity [60, 61]. Blocking the BTLA/HVEM pathway has been shown to enhance IFN-γ production in circulating CD4+ and CD8+ T cells, highlighting its role in peripheral T-cell suppression, as observed in hepatocellular carcinoma [62]. These immune checkpoints affect CD4+ and CD8+ T cells, which are key prognostic indicators in BC.

Patients in the low-risk group demonstrated greater sensitivity to chemotherapeutic agents such as 5-fluorouracil, cisplatin, docetaxel, doxorubicin, paclitaxel, shikonin, and vinorelbine. Chemotherapy exerts its effects by disrupting DNA synthesis in tumor cells and stimulating anticancer immunity through two mechanisms: activating dying tumor cells to release immunostimulatory molecules and directly targeting immune cells to enhance the immune response [63].

Higher immune dysfunction scores, lower immune rejection scores, and higher IPS were observed in the low-risk group. These findings indicate notable differences in immunotherapy responses among patients with BC.

TMB serves as a biomarker for predicting the response to immunotherapy, with higher TMB levels typically associated with greater therapeutic benefits. In this study, the mutation frequencies of 19 genes were found to differ between the high-risk and low-risk groups, which may help explain the elevated TMB observed in the high-risk group. Among the 17 signature genes identified, BACH2, KCNH8, and ARHGAP26 emerged as the most frequently mutated genes in BC samples.

In summary, 17 signature genes were identified, and a risk model was developed to predict BC prognosis. Additionally, insights into the mechanisms by which circadian rhythm influences BC regulation were provided.

Despite the comprehensive nature of our study, there are several limitations that should be acknowledged. First, the sample size, although substantial, is still limited, particularly in terms of representing diverse patient populations and clinical stages. This may affect the generalizability of our findings. Second, while our prognostic model shows promise, the AUC values are currently in the moderate range, indicating room for improvement in predictive accuracy. Finally, the validation of our model relies primarily on retrospective data analysis, and prospective validation in a larger, independent cohort is needed to further confirm its clinical utility. Future work will focus on addressing these limitations through expanded sample collection, refined model development, and prospective validation studies.

In conclusion, this study successfully developed and validated a novel prognostic model based on circadian rhythm-related genes, highlighting its potential to predict outcomes in breast cancer patients. Through comprehensive bioinformatics analysis and multi-level experimental validation, we identified key genes that significantly influence prognosis and demonstrated their impact on the tumor immune microenvironment. Our findings underscore the importance of circadian regulation in breast cancer and provide a foundation for future personalized treatment strategies.

Data availability

No datasets were generated or analysed during the current study.

Abbreviations

BC:

Breast cancer

ER:

Estrogen receptor

DEGs:

Differentially expressed genes

DECRGs:

Differentially expressed circadian rhythm-related genes

GO:

Gene Ontology

KEGG:

Kyoto Encyclopedia of Genes and Genomes

PPI:

Protein–protein interaction

STRING:

Search Tool for the Retrieval of Interacting Genes

LASSO:

Least absolute shrinkage and selection operator

OS:

Overall survival

IC50:

Half-maximal inhibitory concentrations

TIDE:

Tumor immune dysfunction and exclusion

IPS:

Immunophenoscore

TMB:

Tumor mutation burden

DCA:

Decision curve analysis

ROC:

Receiver operating characteristic

AUC:

Area under the curve

NET:

Neutrophil Extracellular Trap

FOX:

Forkhead box

NETs:

Neutrophil Extracellular Traps

References

  1. Fahad Ullah MB. Cancer: current perspectives on the disease status. Adv Exp Med Biol. 2019;1152:51–64.

    Article  CAS  PubMed  Google Scholar 

  2. Zhou Y, et al. Tumor biomarkers for diagnosis, prognosis and targeted therapy. Signal Transduct Target Ther. 2024;9(1):132.

    Article  PubMed  PubMed Central  Google Scholar 

  3. Zhang ZY, et al. Machine learning applications in breast cancer survival and therapeutic outcome prediction based on multi-omic analysis. Yi Chuan. 2024;46(10):820–32.

    PubMed  Google Scholar 

  4. Ye F, et al. Advancements in clinical aspects of targeted therapy and immunotherapy in breast cancer. Mol Cancer. 2023;22(1):105.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Klerman EB, Brager A, Carskadon MA, et al. Keeping an eye on circadian time in clinical research and medicine. Clin Transl Med. 2022;12(12):e1131.

    Article  PubMed  PubMed Central  Google Scholar 

  6. Litinski M, Scheer FA, Shea SA. Influence of the circadian system on disease severity. Sleep Med Clin. 2009;4(2):143–63.

    Article  PubMed  PubMed Central  Google Scholar 

  7. Mason IC, Qian J, Adler GK, et al. Impact of circadian disruption on glucose metabolism: implications for type 2 diabetes. Diabetologia. 2020;63(3):462–72.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Scheer FA, Hilton MF, Mantzoros CS, et al. Adverse metabolic and cardiovascular consequences of circadian misalignment. Proc Natl Acad Sci U S A. 2009;106(11):4453–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Straif K, Baan R, Grosse Y, et al. Carcinogenicity of shift-work, painting, and fire-fighting. Lancet Oncol. 2007;8(12):1065–6.

    Article  PubMed  Google Scholar 

  10. D’Cunha K, Park Y, Protani MM et al. Circadian rhythm disrupting behaviours and cancer outcomes in breast cancer survivors: a systematic review. Breast Cancer Res Treat, 2022.

  11. Briguglio G, Costa C, Teodoro M, et al. Women’s health and night shift work: potential targets for future strategies in breast cancer (Review). Biomed Rep. 2021;15(6):98.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Walker WN, Sprowls SA, Bumgarner JR, et al. Circadian influences on chemotherapy efficacy in a mouse model of brain metastases of breast Cancer. Front Oncol. 2021;11:752331.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Koritala B, Porter KI, Sarkar S, et al. Circadian disruption and cisplatin chronotherapy for mammary carcinoma. Toxicol Appl Pharmacol. 2022;436:115863.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Xue J, Dai Y, Li G, et al. DEC1 directly interacts with Estrogen receptor (ER) alpha to suppress proliferation of ER-positive breast cancer cells. Biochem Biophys Res Commun. 2020;528(4):740–5.

    Article  CAS  PubMed  Google Scholar 

  15. Malla RR, Padmaraju V, Amajala KC, et al. Association between the circadian clock and the tumor microenvironment in breast Cancer. Crit Rev Oncog. 2021;26(3):43–51.

    Article  PubMed  Google Scholar 

  16. Li S, Shui K, Zhang Y, et al. CGDB: a database of circadian genes in eukaryotes. Nucleic Acids Res. 2017;45(D1):D397–403.

    CAS  PubMed  Google Scholar 

  17. Weinstein JN, Collisson EA, Mills GB, et al. The Cancer genome atlas Pan-Cancer analysis project. Nat Genet. 2013;45(10):1113–20.

    Article  PubMed  PubMed Central  Google Scholar 

  18. Love MI, Huber W, Anders S. Moderated Estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.

    Article  PubMed  PubMed Central  Google Scholar 

  19. Yu G, Wang LG, Han Y, et al. ClusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. von Mering C, Huynen M, Jaeggi D, et al. STRING: a database of predicted functional associations between proteins. Nucleic Acids Res. 2003;31(1):258–61.

    Article  Google Scholar 

  21. Liang L, et al. Integration of scRNA-Seq and bulk RNA-Seq to analyse the heterogeneity of ovarian Cancer immune cells and Establish a molecular risk model. Front Oncol. 2021;11:711020.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Yoshihara K, Shahmoradgoli M, Martinez E, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013;4:2612.

    Article  PubMed  Google Scholar 

  23. Newman AM, Liu CL, Green MR, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015;12(5):453–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Mayakonda A, et al. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 2018;28(11):1747–56.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Li C, et al. Identification of a four-gene panel predicting overall survival for lung adenocarcinoma. BMC Cancer. 2020;20(1):1198.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Li M, et al. Recognition of refractory Mycoplasma pneumoniae pneumonia among Myocoplasma pneumoniae pneumonia in hospitalized children: development and validation of a predictive nomogram model. BMC Pulm Med. 2023;23(1):383.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Kos K, de Visser KE. Neutrophils create a fertile soil for metastasis. Cancer Cell. 2021;39(3):301–3.

    Article  CAS  PubMed  Google Scholar 

  28. Poto R, Cristinziano L, Modestino L et al. Neutrophil Extracell Traps Angiogenesis Cancer Biomedicines, 2022,10(2).

  29. Langiu M, Palacios-Acedo AL, Crescence L et al. Neutrophils, Cancer and thrombosis: the new Bermuda triangle in Cancer research. Int J Mol Sci, 2022,23(3).

  30. Wang YF, Han J. OTOR in breast carcinoma as a potent prognostic predictor correlates with cell proliferation, migration, and invasiveness. Biochem Cell Biol. 2019;97(6):750–7.

    Article  CAS  PubMed  Google Scholar 

  31. Shropshire DB, Acosta FM, Fang K, et al. Association of adenosine signaling gene signature with Estrogen receptor-positive breast and prostate cancer bone metastasis. Front Med (Lausanne). 2022;9:965429.

    Article  PubMed  Google Scholar 

  32. Coassolo S, Davidson G, Negroni L, et al. Citrullination of pyruvate kinase M2 by PADI1 and PADI3 regulates Glycolysis and cancer cell proliferation. Nat Commun. 2021;12(1):1718.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Zhao S, Chen SS, Gu Y, et al. Expression and clinical significance of Sushi Domain- containing protein 3 (SUSD3) and Insulin-like growth Factor-I receptor (IGF-IR) in breast Cancer. Asian Pac J Cancer Prev. 2015;16(18):8633–6.

    Article  PubMed  Google Scholar 

  34. Chen HW, Huang XD, Li HC, et al. Expression of FOXJ1 in hepatocellular carcinoma: correlation with patients’ prognosis and tumor cell proliferation. Mol Carcinog. 2013;52(8):647–59.

    Article  CAS  PubMed  Google Scholar 

  35. Tang GY, Tang GJ, Yin L et al. ECRG4 acts as a tumor suppressor gene frequently hypermethylated in human breast cancer. Biosci Rep, 2019,39(5).

  36. Dai Y, Yang L, Sakandar A, et al. Vemurafenib inhibits immune escape biomarker BCL2A1 by targeting PI3K/AKT signaling pathway to suppress breast cancer. Front Oncol. 2022;12:906197.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Dolezal D, Zhang X, Harigopal M. Increased expression of LEF1 and beta-Catenin in invasive micropapillary carcinoma of the breast is associated with lymphovascular invasion and lymph node metastasis. Appl Immunohistochem Mol Morphol. 2022;30(8):557–65.

    Article  CAS  PubMed  Google Scholar 

  38. Li Z, Li Y, Li N, et al. Silencing GOLGA8B inhibits cell invasion and metastasis by suppressing STAT3 signaling pathway in lung squamous cell carcinoma. Clin Sci (Lond). 2022;136(11):895–909.

    Article  CAS  PubMed  Google Scholar 

  39. McCubrey JA, Steelman LS, Chappell WH, et al. Roles of the Raf/MEK/ERK pathway in cell growth, malignant transformation and drug resistance. Biochim Biophys Acta. 2007;1773(8):1263–84.

    Article  CAS  PubMed  Google Scholar 

  40. Li Y, Zhao X, Liu Q, et al. Bioinformatics reveal macrophages marker genes signature in breast cancer to predict prognosis. Ann Med. 2021;53(1):1019–31.

    Article  PubMed  Google Scholar 

  41. Zhang S, Ma D, Wang X, et al. Syntaxin-11 is expressed in primary human monocytes/macrophages and acts as a negative regulator of macrophage engulfment of apoptotic cells and IgG-opsonized target cells. Br J Haematol. 2008;142(3):469–79.

    Article  CAS  PubMed  Google Scholar 

  42. Yi T, Papadopoulos E, Hagner PR, et al. Hypoxia-inducible factor-1alpha (HIF-1alpha) promotes cap-dependent translation of selective mRNAs through up-regulating initiation factor eIF4E1 in breast cancer cells under hypoxia conditions. J Biol Chem. 2013;288(26):18732–42.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Kanellis DC, Espinoza JA, Zisi A et al. The exon-junction complex helicase eIF4A3 controls cell fate via coordinated regulation of ribosome biogenesis and translational output. Sci Adv, 2021,7(32).

  44. Roychoudhuri R, Eil RL, Clever D, et al. The transcription factor BACH2 promotes tumor immunosuppression. J Clin Invest. 2016;126(2):599–604.

    Article  PubMed  PubMed Central  Google Scholar 

  45. Wang X, Xue B, Zhang Y, et al. Up-regulated circBACH2 contributes to cell proliferation, invasion, and migration of triple-negative breast cancer. Cell Death Dis. 2021;12(5):412.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. DeMille D, Grose JH. PAS kinase: a nutrient sensing regulator of glucose homeostasis. IUBMB Life. 2013;65(11):921–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Zhou Y, Tian Q, Gao H, et al. Correlation between immune-Related genes and Tumor-Infiltrating immune cells with the efficacy of neoadjuvant chemotherapy for breast Cancer. Front Genet. 2022;13:905617.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Zhang L, Zhou A, Zhu S, et al. The role of GTPase-activating protein ARHGAP26 in human cancers. Mol Cell Biochem. 2022;477(1):319–26.

    Article  CAS  PubMed  Google Scholar 

  49. Katoh M, Katoh M. Characterization of human ARHGAP10 gene in Silico. Int J Oncol. 2004;25(4):1201–6.

    Article  CAS  PubMed  Google Scholar 

  50. Zou A, Lin Z, Humble M, et al. Distribution and functional properties of human KCNH8 (Elk1) potassium channels. Am J Physiol Cell Physiol. 2003;285(6):C1356–66.

    Article  CAS  PubMed  Google Scholar 

  51. Mallya M, Campbell RD, Aguado B. Characterization of the five novel Ly-6 superfamily members encoded in the MHC, and detection of cells expressing their potential ligands. Protein Sci. 2006;15(10):2244–56.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Fang H, Guo Z, Chen J, et al. Combination of epigenetic regulation with gene therapy-mediated immune checkpoint Blockade induces anti-tumour effects and immune response in vivo. Nat Commun. 2021;12(1):6742.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Chen X, Cubillos-Ruiz JR. Endoplasmic reticulum stress signals in the tumour and its microenvironment. Nat Rev Cancer. 2021;21(2):71–88.

    Article  CAS  PubMed  Google Scholar 

  54. Xu D, Liu Z, Liang MX, et al. Endoplasmic reticulum stress targeted therapy for breast cancer. Cell Commun Signal. 2022;20(1):174.

    Article  PubMed  PubMed Central  Google Scholar 

  55. Liu J, Fan L, Yu H, et al. Endoplasmic reticulum stress causes liver Cancer cells to release Exosomal miR-23a-3p and Up-regulate programmed death ligand 1 expression in macrophages. Hepatology. 2019;70(1):241–58.

    Article  CAS  PubMed  Google Scholar 

  56. Obiedat A, Seidel E, Mahameed M, et al. Transcription of the NKG2D ligand MICA is suppressed by the IRE1/XBP1 pathway of the unfolded protein response through the regulation of E2F1. FASEB J. 2019;33(3):3481–95.

    Article  CAS  PubMed  Google Scholar 

  57. Gil DAC, Huh SJ, Ekram MB, et al. Immune escape in breast Cancer during in situ to invasive carcinoma transition. Cancer Discov. 2017;7(10):1098–115.

    Article  Google Scholar 

  58. De Larco JE, Wuertz BR, Furcht LT. The potential role of neutrophils in promoting the metastatic phenotype of tumors releasing interleukin-8. Clin Cancer Res. 2004;10(15):4895–900.

    Article  PubMed  Google Scholar 

  59. Wculek SK, Malanchi I. Neutrophils support lung colonization of metastasis-initiating breast cancer cells. Nature. 2015;528(7582):413–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Buller CW, Mathew PA, Mathew SO. Roles of NK cell receptors 2B4 (CD244), CS1 (CD319), and LLT1 (CLEC2D) in Cancer. Cancers (Basel), 2020,12(7).

  61. Zippelius A, Schreiner J, Herzig P, et al. Induced PD-L1 expression mediates acquired resistance to agonistic anti-CD40 treatment. Cancer Immunol Res. 2015;3(3):236–44.

    Article  CAS  PubMed  Google Scholar 

  62. Carter JM, Polley MC, Leon-Ferre RA, et al. Characteristics and spatially defined immune (micro)landscapes of Early-stage PD-L1-positive Triple-negative breast Cancer. Clin Cancer Res. 2021;27(20):5628–37.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. Galluzzi L, Humeau J, Buque A, et al. Immunostimulation with chemotherapy in the era of immune checkpoint inhibitors. Nat Rev Clin Oncol. 2020;17(12):725–41.

    Article  PubMed  Google Scholar 

Download references

Acknowledgements

This work was supported by Key Laboratory of Dongzhimen Hospital of Beijing University of Traditional Chinese Medicine.

Funding

the Young Faculty Project of Beijing University of Traditional Chinese Medicine, Funding ID: 2022-BUCMXJKY009.

National Natural Science Foundation of China, ID: 82474510.

Author information

Author notes

  1. Ling Wang and Xiang Gao contributed equally to this work.

Authors and Affiliations

  1. Department of Breast Surgery, Dongzhimen Hospital, Beijing University of Chinese Medicine, Haiyuncang 5th, Dongcheng District, Beijing, 100700, China

    Ling Wang, Xiang Gao, Ximeng Zuo, Tangshun Wang & Xiaoguang Shi

Authors
  1. Ling Wang
    View author publications

    You can also search for this author inPubMed Google Scholar

  2. Xiang Gao
    View author publications

    You can also search for this author inPubMed Google Scholar

  3. Ximeng Zuo
    View author publications

    You can also search for this author inPubMed Google Scholar

  4. Tangshun Wang
    View author publications

    You can also search for this author inPubMed Google Scholar

  5. Xiaoguang Shi
    View author publications

    You can also search for this author inPubMed Google Scholar

Contributions

Conception and design of the research and statistical analysis: Ling WangStatistical analysis and Writing of the manuscript: Xiang Gao, Ximeng Zuo Acquisition of data and analysis and interpretation of the data: Tangshun WangCritical revision of the manuscript for intellectual content and Obtaining financing: Xiaoguang ShiAll authors read and approved the final draft.

Corresponding author

Correspondence to Xiaoguang Shi.

Ethics declarations

Ethics approval and consent to participate

This study was conducted with approval from the Ethics Committee of Dongzhimen Hospital, Beijing University of Chinese Medicine. This study was conducted in accordance with the declaration of Helsinki.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Electronic supplementary material

Below is the link to the electronic supplementary material.

12957_2025_3829_MOESM1_ESM.tif

Supplementary Material 1: Additional file 1: The box plot of the differential expression of prognostic genes in the training set was generated.

12957_2025_3829_MOESM2_ESM.tif

Supplementary Material 2: Additional file 2: The KM survival curve of prognostic genes was plotted. Abbreviations: KM: Kaplan Meier.

Supplementary Material 3

Supplementary Material 4

Supplementary Material 5

Supplementary Material 6

Supplementary Material 7

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Wang, L., Gao, X., Zuo, X. et al. Prognostic value of circadian rhythm-associated genes in breast cancer. World J Surg Onc 23, 186 (2025). https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12957-025-03829-8

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s12957-025-03829-8

Keywords

World Journal of Surgical Oncology

ISSN: 1477-7819

Contact us