The circadian rhythm gene set (comprising 210 genes) was retrieved from the Circadian Genes Database (CGDB) (http://cgdb.biocuckoo.org) on 27 June 2023. Concurrently, the BC RNA-Seq data (TPM values) and clinicopathological data (1,086 cancer samples and 99 adjacent normal tissue samples) were acquired from the TCGA database (https://portal.gdc.cancer.gov). The expression levels of the circadian rhythm genes in BC patients were screened and subjected to log2 transformation. The external validation dataset GSE42568 (104 cancer samples and 17 adjacent normal tissue samples) and GSE58812 (107 cancer samples)was obtained from the GEO (Gene Expression Omnibus) (https://www.ncbi.nlm.nih.gov/geo) database. All the datasets in this study originated from public databases and thus did not require ethical approval.

The gene mutation data and copy number variation data of BRCA were downloaded from the TCGA database and corresponded to the circadian rhythm genes. They were processed and visualized respectively by the “maftools” and “Rcircos” packages of R-4.3.0 software. The “DESeq2″ package was employed to screen the differentially expressed genes (DEGs) of circadian rhythm with the criteria of P < 0.05 and | log2(FoldChange)| > 1, and the “pheatmap”, “ggplot2″ and “psych” packages were utilized to draw the volcano plot and the heatmap of the correlation between genes and clinical traits.

The “enrichplot”, “org.Hs.e.g.,.db” and “clusterProfiler” packages were employed to conduct Gene Ontology (GO) functional enrichment analysis, Kyoto Encyclopedia of Genes and Genomes pathway (KEGG) analysis and Gene Set Enrichment Analysis (GSEA) on the DEGs. A P value <0.05 was regarded as statistically significant, and the results were visualized using the “ggplot2″ package. The protein-protein interaction network (PPI) of the differential genes was generated via The STRING website (www.string-db.org) with a minimum correlation score of 0.400 as the threshold. The data was downloaded from the database in TSV format and visualized using Cytoscape_v3.9.0 software.

The expression levels of circadian rhythm DEGs and survival information of 1069 BC samples in the training group were consolidated. Patients without overall survival (OS) and survival status (Censor) were excluded. The remaining 1049 BC samples were successively incorporated into Cox analysis and three machine learning algorithms (LASSO regression, Random Forest, SVM-RFE). For the assessment of the prognostic value in this study, the intersection of DEGs calculated by the four algorithms was selected, and ultimately a risk score (Risk Score) prediction model was constructed. The formula for the risk score prediction model is: Risk Score = Exp1 × C1 + Exp2 × C2 + … + Expn × Cn (Exp: the expression level of the gene, C: the regression coefficient obtained by LASSO regression analysis, n: the number of intersection genes). The Risk Score of each sample was calculated based on the above formula, and the BC samples were classified into low-risk and high-risk groups using the median Risk Score as the cut-off value. A P value <0.05 was considered statistically significant.

Firstly, packages such as “survival” and “survminer” were exploited to draw the Kaplan-Meier (K-M) survival curve for comparing the overall survival (OS) of the low-risk group and the high-risk group. Secondly, the time-dependent receiver operating characteristic (ROC) curve was plotted through the “time ROC” package, and the area under the curve (AUC) values of the OS of the training group samples at 1, 3, and 5 years were calculated to appraise the accuracy of the predictive risk model. Finally, the clinicopathological factors of the two groups and the Risk Score were incorporated into the Cox regression analysis. Packages like “rms” were loaded to construct the nomogram with the variables presenting independent risk factors. The upper part represented the scoring system, and the lower part constituted the prediction system. The 1-year, 3-year, and 5-year survival rates of patients could all be predicted via the scoring system, and the calibration curve and ROC curve were employed to evaluate the accuracy of survival prediction. The GSE42568 + GSE58812 dataset was utilized as the validation set to verify the constructed prediction model.

The infiltration status of 22 distinct immune cells was evaluated by applying packages such as “CIBERSORT” to explore the relationship between the risk score and immune cell infiltration. Additionally, 46 common immune checkpoints were retrieved from relevant literature to investigate the association between the risk score and the 46 immune checkpoint genes. A P value <0.05 indicated a significant threshold.

The expression and drug data of relevant genes were retrieved from the CellMiner database (https://discover.nci.nih.gov/cellminer/home.do). The drug data were filtered through a combined method involving clinical laboratory validation and Food and Drug Administration (FDA) standard certification. Then, the expression data of the prognostic circadian rhythm genes were merged with the drug data, and a Pearson correlation test was conducted to determine the correlation and drug sensitivity between them.

After establishing a one-to-one correspondence between the genetic variation data and the circadian rhythm genes, the data of gene mutations and copy number variations obtained (199 genes, 534 samples) revealed a somatic mutation frequency of 55.43% (296/534 cases) in Breast cancer (BC) (Supplementary Figure S1). The circadian rhythm genes of BC patients collected from the TCGA database were screened for 62 differentially expressed genes (DEGs) with a threshold of P < 0.05 and |log2(Fold Change) | > 1, encompassing 30 upregulated genes and 32 downregulated genes. The visualization results of the volcano plot and the correlation heatmap between the differential genes and clinical traits are presented (Figures 1A, B). Meanwhile, the chromosomal locations of the circadian rhythm DEGs were depicted (Supplementary Figure S2).

To explore the pertinent biological functions and pathways of the circadian rhythm genes in BC, GO and KEGG analyses were implemented on the aforementioned 62 differentially expressed genes (DEGs), and the GO classifications and KEGG pathways with the highest enrichment degrees were presented (Figures 1C, D). The GO analysis indicated that the DEGs were primarily enriched in biological functions related to the circadian rhythm process, chromosome, cell junction, and molecular transducer activity, among others. The KEGG analysis mainly manifested enrichment in biological functions associated with neuroactive ligand-receptor interaction pathways, cAMP signaling pathways, Rap1 signaling pathways, and the like. Circadian rhythm genes with P < 0.05 were incorporated into the Gene Set Enrichment Analysis (GSEA) (Figure 1E), and the analysis outcomes were mainly enriched in relevant biological processes such as non-alcoholic fatty liver disease, prion diseases, blood lipids and atherosclerosis, and transcriptional dysregulation in cancer, further validating the intimate relationship between the regulation of genetic material and the metabolism of energy substances such as fatty acids and tumor progression.

Subsequently, after eliminating the DEGs that did not interact with other nodes, a protein-protein interaction (PPI) network of 53 nodes was constructed through the Degree algorithm, revealing the interactions among the circadian rhythm DEGs proteins. The darker the color, the greater the number of nodes associated with the gene (Figure 1F). The results of the above analyses elucidate the relationships at the gene level and also furnish distinct research concepts for identifying potential targets of biomarkers related to prognosis in the future.

The aforementioned 62 differentially expressed genes (DEGs) were respectively incorporated into Cox analysis, LASSO regression, Random Forest, and SVM-RFE. Cox analysis yielded 11 genes that were significantly associated with independent prognosis (P < 0.05) (Figure 2A). Subsequently, the dimensionality reduction analysis of LASSO regression was carried out, and the minimum value of logλ was indicated by the dotted line, obtaining 46 (seed 10,000) feature genes (Figures 2B, C). Then, after ranking the importance scores of the DEGs by applying Random Forest, 30 candidate genes were identified (Figure 2D), and 30 genes with the highest accuracy were identified by the SVM-RFE method based on their importance scores (Figure 2E). Finally, the DEGs screened by each method (Cox analysis, n = 11; LASSO regression, n = 46; Random Forest, n = 30; SVM-RFE, n = 30) were intersected, and the 5 prognosis-related feature genes (SUV39H2, OPN4, RORB, FBXL6, SIAH2) after taking the intersection were visualized as a Venn diagram (Figure 2F). The RiskScore of each sample was calculated based on the expression levels and risk coefficients of the 5 genes: RiskScore = (Expression level of SUV39H2 × 0.0436) + (Expression level of OPN4 × 1.4270) + (Expression level of RORB × 0.1917) + (Expression level of FBXL6 × 0.3190) + (Expression level of SIAH2 × −0.1984).

The median RiskScore of the training group BC samples was employed as the cut-off value to categorize them into the low-risk group (L = 525) and the high-risk group (H = 524). The relationship between the RiskScore and the clinicopathological characteristics was depicted through box plots and heat maps, and patients in different risk groups manifested distinct clinicopathological characteristics (Figures 3A–E). The pathological stage, T stage, and N stage were all correlated with the RiskScore, whereas age and gender were not associated with the RiskScore. Furthermore, we also presented the correlation scatter plots between the RiskScore and the 5 prognostic feature genes (Figures 3F–I, Supplementary Figure S3). SUV39H2, OPN4, RORB, and FBXL6 were positively correlated with the RiskScore, while SIAH2 was negatively correlated with the RiskScore.

The median RiskScore categorized the training group samples into the low-risk group and the high-risk group, and the Kaplan-Meier (K-M) survival curve was plotted. The blue color denoted the low-risk group, and the red color represented the high-risk group. The results demonstrated a statistically significant difference in overall survival (OS) between the two groups (P < 0.0001), suggesting that the high-risk group was more inclined to have an unfavorable prognosis than the low-risk group (Figure 4A). To guarantee the stability of the risk score prediction model, the receiver operating characteristic (ROC) curve and the time-dependent ROC curve were constructed, and the area under the curve (AUC) value was calculated to offer predictions of patient survival (Figures 4B, C). Utilizing the same approach, the samples of the validation group dataset were classified into the low-risk group (L = 106) and the high-risk group (H = 105), followed by the plotting of the K-M survival curve, the ROC curve, and the time-dependent ROC curve. A statistically significant difference in OS persisted between the two groups (p < 0.0001) (Figures 4D–F). The aforementioned results indicated that the model maintained favorable predictive performance in the validation group. By assessing the efficacy of the risk model based on circadian rhythm genes in clinical practice, RiskScore and clinicopathological characteristics were incorporated into the Cox regression analysis to identify variables with independent risk factors, and a Nomogram was developed (Figure 5A). This individualized prediction model could estimate the survival rates of BC patients (1-year, 3-year, and 5-year). It was conspicuously evident from the calibration curve that the outcomes of the Nomogram and the actual observed results exhibited precise overlap in both the training group and the validation group (Figure 5B). The ROC curve also attested that the prediction of the Nomogram was adequately reliable (Figure 5C).

The CIBERSORT algorithm was utilized to compute the correlations between the 5 prognostic feature genes and 22 distinct types of immune cells in order to assess the linkage between the risk score prediction model and immune cell infiltration. In this research, the expression of SUV39H2 was predominantly positively correlated with the infiltration of neutrophils, activated dendritic cells, resting dendritic cells, and M2 macrophages; the expression of SIAH2 was mainly negatively correlated with the infiltration of neutrophils, activated dendritic cells, and activated natural killer (NK) cells; the expression of RORB was mainly positively correlated with the infiltration of resting mast cells and M0 macrophages; the expression of OPN4 was mainly positively correlated with the infiltration of M0 macrophages and resting mast cells; the expression of FBXL6 was mainly positively correlated with the infiltration of CD8 T cells, resting CD4 memory T cells, resting NK cells, M1 macrophages, and γδ T cells (Figure 6A). Furthermore, the levels of immune cell infiltration between the low-risk group and the high-risk group were compared (Figure 6B). Among the infiltration proportions of 22 types of immune cells, 20 immune cells manifested significant discrepancies between the low-risk group and the high-risk group. Particularly, macrophages occupied a considerable proportion in both groups and exhibited significant differences, suggesting that controlling the behavior of this cell is indispensable for intervening in tumor progression in BC patients.

Within the BC dataset, we delineated the interactions between the risk score and 46 common immune checkpoint genes. The expressions of immune checkpoint genes in the low-risk group and the high-risk group are depicted as follows (Figure 6C). Based on the correlation analysis, genes such as SUV39H2, OPN4, RORB, and FBXL6 were positively correlated with the RiskScore, while genes such as SIAH2 were negatively correlated with the RiskScore (Figure 6D). The aforementioned results might offer potential directions for seeking biological targets related to immunotherapy.

Pearson correlation analysis was conducted between the five prognostic feature genes and the drug data. The drug sensitivity based on the high and low expression levels of the five prognostic feature genes was calculated (Figures 7A–L, Supplementary Figures S4, 5). We found that the high or low expression of SIAH2 was significantly associated with sensitivity to drugs such as Ifosfamide (p < 0.001). Through analysis, no significant correlations between drug sensitivity and expression were found for the other four genes.