We downloaded mRNA sequencing data and clinical information of BC patients from the Cancer Genome Atlas (TCGA) database as a training dataset (TCGA-BRCA; https://portal.gdc.cancer.gov/) Data integration and processing were as follows :1) normal samples were removed; 2) samples with a survival time less than 30 days were screened; 3) only one sample was retained for each patient, and 4) genes expressed in more than half of the samples were retained. After careful screening, 937 BC samples were used for follow-up studies. Gene expression profiles and clinical information of BC patients in the GSE21653, GSE20685, and GSE58812 datasets were used as test data sets from Gene Expression Omnibus (GEO) data (https://www.ncbi.nlm.nih.gov/geo/) downloaded and processed according to the same filtering criteria. Single-cell RNA sequencing (scRNA-seq) data (GSE188600) from triple-negative BC patients were also downloaded from the GEO database for verification. The specific data download and processing process is included in Supplementary Table S1. Moreover, a total of 1369 CRRGs about Homo sapiens were obtained from the Circadian Gene Data Base (GCDB, http://cgdb.biocuckoo.org/). These genes have been used as the basis for further research. The detailed flow chart of the research design is shown in Figure 1.

In this study, we first performed a univariate COX regression analysis on genes expressed in TCGA-BRCA, then identified genes with p < 0.05 as overall survival (OS)-related genes. To find candidate genes, these genes were crossed with 1369 CRRGs. The least absolute contraction and selection operator (LASSO) regression was then used to determine non-zero coefficients, allowing us to eliminate potential predictors and select the best OS-related genes while avoiding model overfitting. Finally, multivariate Cox regression was used to identify and calculate correlation coefficients for candidate genes involved in the final modeling. The associated risk score of the prognostic CRRGs was equal to the product of the prognostic rhythm-related gene’s expression and its coefficient.

As previously stated, the training dataset was TCGA-BRCA, while the test dataset was GSE20685 for OS of BC. Based on the median calculated risk scores associated with rhythm genes, samples in each dataset were divided into high-risk and low-risk groups. The Kaplan-Meier survival analyses were then plotted to investigate significant differences in patient survival. Furthermore, receiver operating characteristic (ROC) curves for 1, 2, and 3 years were generated in R software to calculate the area under the ROC curve (AUC) values for each prediction model to further evaluate the model’s efficiency and accuracy. Simultaneously, disease-free survival (DFS) and metastasis-free survival (MFS) of BC were validated using outcome events and time in GSE21653 and GSE58812, respectively, and related graphs were created. In addition, validation of scRNA-seq for BC was performed in the GSE188600 dataset.

Clinical information was first extracted from the TCGA-BRCA dataset, including age, stage, and risk group. Univariate and multivariate COX regression analyses were then performed to identify independent prognostic factors affecting OS in BC patients. The corresponding operations were then performed in the validation set databases (GSE21653, GSE20685, and GSE58812) to make our extracted information universally convincing (due to missing data, the Stage was only validated in GSE20685). To visualize the Cox regression results, we used the “rms” R package to generate a Normogram that included CRRGs and clinical variables. Age, TNM stage, and CRRGs risk groups were three independent predictors. The scoring criteria were developed based on the magnitude of their regression coefficients, allowing 1-year, 2-year, and 3-year OS probabilities to be calculated. Calibration curves were used to visually demonstrate the consistency of the Normogram model at 1, 2, and 3 years.

Differentially expressed genes (DEGs) between high- and low-risk groups were first identified by the “limma” package with screening criteria of |log FC| > 1 and p-value_ t < 0.05. This DEG list was then used for gene ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis for TCGA-BRCA to explore possible biological functions and signaling pathways. GO analysis included biological process (BP), cellular composition (CC), and molecular function (MF) (p < 0.05 was statistically significant). The results of the GO analysis and KEGG pathway were visualized using the R package “ggplot2”.

To further explore the potential molecular mechanisms of rhythm-related genes in breast cancer formation and progression, we performed a gene set enrichment analysis (GSEA) on the BC dataset in the TCGA database to explore the enrichment pathways associated with rhythm genes in high-risk and low-risk groups. “c2. cp.kegg.v7.4. symbols” was chosen to be used for our analysis.

CIBERSOFT was used to explore the relationship between high- and low-risk groups based on candidate gene characteristics and the abundance of each immune cell type. Information on the 22 immune cell types can be downloaded from the attachment to the previous article (https://www.nature.com/articles/nmeth.3337#MOESM207). Box plots were then plotted using the ggplot2 R package to visually represent differences in abundance. Additionally, we investigated the relationship between CCGRs and immunological checkpoints (PD1, PDL2, and CTLA4) as well as the connection between risk scores for CCGRs and immune checkpoints.

To investigate the response of BC patients in high- and low-risk groups to ICIs treatment, we obtained immunophenotype score (IPS) data for BC from The Cancer Immunome Atlas (TCIA, https://tcia.at/), which quantifies tumor immunogenicity scores ranging from 0 to 10. IPS values are associated with tumor immunogenicity and can be used to predict patient response to ICIs treatment (Charoentong et al., 2017). Furthermore, we compared the expression of immune checkpoints in high- and low-risk groups.

Statistical analyses were carried out using R (version 4.2.0) and SPSS 26.0. The independent prognostic value of clinical characteristics and risk groups associated with CRRGs was described using univariate and multifactorial COX. Pearson correlation coefficients were used to examine the relationship between CRRGs and immune checkpoints. The presence of two-sided p values less than 0.05 was considered statistically significant.

The flow chart in Figure 1 shows that we finally constructed a prognostic model associated with 18 CRRGs, which include STXBP5, CYP27A1, TAGLN2, SIPA1L1, ZNF485, IFNG, FOXJ1, TP53I11, CEACAM1, TCF7, PPA2, MAK, EMP1, NDRG2, DGAT1, RGL3, TULP4, and CABYR. Figures 2A, B shows the LASSO regression process. The relationship between these 18 candidate genes and the core clock genes is shown in Supplementary Figure S3. A prognostic risk score formula was developed based on the linear combination of CRRGs expression levels and the weighted regression coefficients of multiple Cox regression analysis: Risk score = 0.557*STXBP5-0.189*CYP27A1+0.70*TAGLN2-0.814*SIPA1L1+0.747*ZNF485–1.192*IFNG-0.139*FOXJ1+0.314*TP53I11 0.238*CEACAM1+0.373*TCF7+0.577*PPA2-0.253*MAK+0.333*EMP1-0.220*NDRG2+0.366*DGAT1-0.152*RGL3+0.350*TULP4-0.237*CABYR. These candidate genes were categorized into risk types (STXBP5, TAGLN2, ZNF485, TP53I11, TCF7, PPA2, EMP1, DGAT1, TULP4) and protective types (CYP27A1, SIPA1L1, IFNG, CEACAM1, RGL3), with HR > 1 (p < 0.05) being associated with poor prognosis and HR < 1 (p < 0.05) being associated with better prognosis (Figure 2C). Interestingly, NDRG2 and CABYR were present and expressed at higher levels in the normal group despite FOXJ1, MAK, and EMP1 not being in the human protein atlas (HPA; https://www.proteinatlas.org/), proving that these two genes are also protective. In addition, the expression differences of the other 15 CRRGs in normal and BC tissues were indexed in the HAP database and shown in Supplementary Figure S1. Clinical stage was found to be positively correlated with risk score (Figure 2E), with a statistically significant difference (p = 0.01). Based on the median risk score, patients were divided into high-risk and low-risk groups, with the expression of 18 candidate genes in TCGA in both groups shown in Figure 2D and the expression in the three validation sets shown in Supplementary Figure S6.

To determine the prediction accuracy of the 18 candidate genes, OS validation was performed using GSE20685 datasets, and DFS and MFS validation were performed using two other external datasets (GSE21653 and GSE58812), respectively. The AUCs for 1-year, 2-year, and 3-year survival in the TCGA test set were 0.83, 0.82, and 0.82 (Figure3A), whereas in GSE20685, GSE21653, and GSE58812, the AUCs for 1-year, 2-year, and 3-year survival were 0.80, 0.69, and 0.65; 0.63, 0.65, and 0.70; 0.75, 0.65, and 0.65, respectively (Figure 3B, Figures 4A, B).

Kaplan-Meier survival analysis showed that patients in the high-risk group had a significantly shorter survival time than those in the low-risk group (TCGA: p < 0.0001; GSE20685: p = 0.028; GSE21653: p = 0.00061; GSE58812: p = 0.0026) (Figure 3C, D; Figures 4C, D). Moreover, risk scores were negatively correlated with survival time in all four data sets (Figures 3E, F, Figures 4E, F). In addition, the results of scRNA-seq are shown in Supplementary Figure S4, 5.

First, the risk score and other clinical parameters were used as covariates in univariate and multivariate Cox proportional risk regression models. The risk group was found to perform well as a prognostic factor affecting the survival cycle of BC patients in both the TCGA and validation set databases. Table 1; Table 2 also showed the results for age and stage. As a result, we selected the three variables listed above to create follow-up nomograms. The findings revealed that as the total score increased, the survival time decreased (Figure 5A). The calibration curve revealed that the nomogram model had high predictive accuracy for OS at 1, 2, and 3 years and significantly overlapped with the 45-degree angle line (Figures 5B–D). Additionally, Supplementary Table S1 displayed the survival times (OS, MFS, and DFS) of the three external validation sets (GSE21653, GSE20685, and GSE58812) in relation to CRRGs.

Based on the results of the “limma” variance analysis, we identified 354 differentially expressed genes (DEGs), including 193 up-regulated DEGs and 161 down-regulated DEGs. The bubble diagram shows the terms of these DEGs in GO-BP, GO-CC and GO-MF (Figure 6A). GO functional annotation shows that these DEGs are mainly involved in cytoplasmic translation, ribosomes, and ribosome structural composition. A two-way bar graph showing the top five pathways of KEGG upregulation and downregulation (Figure 6B). GSEA showed that the significantly enriched pathway in the high-risk group was systemic lupus erythematosus, whereas in the low-risk group it was the adipocytokine signaling pathway, ascorbate and aldehyde metabolism, pentose and glucuronide interconversion, phenylalanine metabolism, and ribosomes (Figures 6C–H). Furthermore, using a risk-score model, we investigated the genetic alterations of these CRRGs and discovered that they are relatively conserved evolutionarily (mutation rates of 1% or less) (Supplementary Figure S2), which is consistent with the current study (Cederroth et al., 2019).

We investigated the level of infiltration of 22 immune cells in various risk groups, and box plots revealed that the level of immune infiltration of T-cell CD8 and T-cell CD4 memory activated was significantly lower in the high-risk group, while immunosuppressive cells such as macrophages M0 and M2 were significantly higher (Figure 7A). The high-risk group appeared to have an immunosuppressive tumor microenvironment filled with a large number of immunosuppressive cells, which was consistent with a poor prognosis. The correlation heat map between CRRGs and immune checkpoints is shown in Figure 7B, which shows that IFNG and TCF7 are strongly correlated with immune checkpoints (Figures 7D, E). These two genes were also shown in scRNA analysis to be distributed mainly in T-cell (Supplementary Figure S4E). The chord plot of the correlation between risk scores and immune checkpoints is represented in Figure 7C.

Because the TCGA-BRCA dataset lacked information on ICIs treatment, we used two IPS-valued subtypes (IPS- PD-L1/PD-1/PD-L2 blocker and IPS-CTLA-4 blocker) as proxies for response to anti-PD-1/PD-L1 and anti-CTLA-4 treatment in BC patients. In the CRRGs prediction model, the low-risk group had a higher relative probability of receiving anti-PD-1/PD-L1 and anti-CTLA-4 treatment (Figures 8A, B). The findings imply that patients with low CRRG risk scores may be candidates for ICIs treatment. Furthermore, we compared immune checkpoint expression levels between the high-risk and low-risk groups. Patients with low risk had significantly higher levels of PD-1, PD-L2, and CTLA-4 (Figure 8C).