Integration of scRNA-Seq and Bulk RNA-Seq Identifies Circadian Rhythm Disruption-Related Genes Associated with Prognosis and Drug Resistance in Colorectal Cancer Patients

Introduction

Colorectal cancer (CRC) is one of the most common malignancies worldwide, with its incidence continuously rising. This increase is partly attributed to the emergence of new favorable subsets of cancers of unknown primary (CUP), including CRC CUP, which are treated as CRC and contribute to the current increased incidence of CRC.¹ With the continuous rise in its incidence, early diagnosis and treatment of CRC have become increasingly crucial. Notably, the population of elderly patients with CRC is growing, and while older patients are more prone to severe postoperative complications, there is no consensus that age itself affects survival outcomes. The prognosis of older patients may be confounded by differences in stage at presentation, tumor site, preexisting comorbidities, and type of treatment received.² Current standard treatments for CRC include surgery, chemotherapy, radiotherapy, and targeted therapies. Identifying effective prognostic biomarkers is essential for the development of new CRC drug targets, optimizing their use, and enabling personalized treatment approaches, thereby improving patient outcomes and increasing survival rates. In this context, it is important to note that immune cell PD-L1 expression is significantly higher in mismatch repair (MMR)-deficient (MSI-H) CRC compared to MMR-proficient (MSI-L) tumors, with no differences among the different MSI-H molecular subtypes.³ The recommended screening for defective DNA mismatch repair includes immunohistochemistry (IHC) and/or MSI testing. However, there are challenges in distilling the biological and technical heterogeneity of MSI testing down to usable data. It has been reported in the literature that IHC testing of the MMR machinery may give different results for a given germline mutation, and this discrepancy has been suggested to be due to somatic mutations.

Circadian rhythm is a natural biological clock that maintains a 24-hour cycle of activity by regulating the molecular clock mechanisms in various organs and tissues.^4–6 It involves cycles such as sleep-wake, feeding-fasting, and activity-rest, ensuring the coordination of behavior and physiological functions across the body, thereby maintaining homeostasis. However, when circadian rhythms are disrupted, it can lead to various functional disorders, such as abnormal cell proliferation, increased gene mutations, and resistance to apoptosis, all of which can elevate the risk of various diseases, including cancer.⁷ Additionally, circadian rhythm disruption, through mechanisms such as apoptosis, cell cycle regulation, and metabolic control, directly or indirectly contributes to the initiation and progression of various cancers.^8–10 For instance, Zhu et al reviewed the significant functions of the circadian rhythm system coordinated by interlocking CLOCK genes in regulating processes from sleep-wake cycles to immune responses in endocrine tumors.¹¹ Zeng et al found an epidemiological link between cancer risk and circadian rhythm gene expression, associated with circadian components in tumor cells and their microenvironment.¹²

Previous studies have utilized circadian rhythm-related genes (CRDGs) to construct prognostic models for various cancers, achieving high predictive accuracy. For example, Hao et al developed a prognostic signature based on six CRDGs that could precisely evaluate the immune therapy response and prognosis in patients with head and neck squamous cell carcinoma.¹³ Li et al identified and validated the clinical relevance of circadian rhythm characteristics and tumor microenvironment in non-small cell lung cancer (NSCLC), revealing that circadian rhythm may play a critical role in NSCLC.¹⁴ However, the role of CRDGs in CRC remains unclear. Therefore, this study aimed to infer the CRD status from single-cell RNA sequencing (scRNA-seq) and bulk RNA sequencing data in CRC and to construct a CRC prognostic model based on CRD scores. This model effectively predicts patient prognosis following immunotherapy and offers new strategies for the development of novel CRC targets.

Methods

Data Acquisition and Preprocessing

This study retrieved two scRNA-seq datasets of CRC samples from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/), specifically GSE178318 and GSE205506. The GSE178318 dataset comprises transcriptomic data from 111,292 single cells, including six primary CRC samples, six matched liver metastasis samples, and three peripheral blood mononuclear cell samples. The GSE205506 dataset includes 40 samples from 19 CRC patients treated with mismatch repair-deficient and microsatellite instability-high (d-MMR/MSI-H) therapy. The scRNA-seq data were processed and analyzed using the “Seurat” package in R. Cells included in the analysis met the following criteria: gene count per cell exceeding 200, gene count per cell less than the 98th percentile across all cells, mitochondrial gene expression percentage below 15%, and total RNA count per cell less than the 98th percentile across all cells. Subsequently, data normalization was performed using the built-in functions of the Seurat package. Batch effects were removed using the “harmony” package in R, and clustering was performed based on dimensionality reduction values.

For bulk RNA-seq data, we downloaded level 3 data from the UCSC Xena (http://xena.ucsc.edu/), including mRNA expression and clinical data from 22 cancer types (eg, LUAD, BRCA, and PRAD). Additionally, we downloaded several CRC datasets from the GEO database, specifically GSE39582 and GSE72970, which include expression profiles from 566 and 124 CRC samples, respectively. Furthermore, bulk RNA-seq data from 83 neuroblastoma samples (GSE221103), 76 osteosarcoma samples (GSE221173), four samples from PD-1 (nivolumab)-responsive patients, seven non-responsive renal cell carcinoma samples (GSE6750), 26 ovarian cancer samples (GSE188249), and 191 advanced clear cell renal cell carcinoma samples (PMC7499153) were collected from the GEO database. For microarray data, the robust multi-array average (RMA) algorithm was applied to normalize the microarray signals, and the data were imported into R (http://www.r-project.org) for further analysis.

Weighted Gene Co-Expression Network Analysis (WGCNA) for scRNA-Seq

The scRNA-seq data from primary and liver metastatic CRC samples were grouped and analyzed using the “WGCNA” package in R. This study aimed to identify signature networks positively correlated with liver metastasis from CRDGs. Optimal gene filtering for single-cell data was performed using default parameters, followed by quantification of the correlation between co-expression gene modules and clinical features, as well as the identification of gene modules associated with liver metastasis using the “moduleEigengenes” function.

CRDG Scoring and Cell-Cell Communication Analysis

CRDG scoring for each cell was performed according to a previously described method.¹⁵ For cell-cell communication analysis, the “CellphoneDB” package in R was employed to assess the mean expression levels of annotated ligand-receptor (L-R) pairs from the STRING database (https://cn.string-db.org/). Significant interactions between different cell types were analyzed (p<0.05). Cytokine signaling activity was evaluated using the CytoSig tool.

CRDG-Based Model for Immunotherapy Response

This study extracted CRDG expression levels from four representative cohorts receiving immunotherapy. Prognosis-related genes were selected using the Lasso algorithm, and a logistic regression model was developed to predict patient outcomes. The model provides a risk score for each patient, and the differences in post-immunotherapy scores between non-responders and responders were further evaluated.

Statistical Analysis

Statistical analyses were conducted using R software (v4.0). Differences in continuous variables between two groups were analyzed using the Wilcoxon rank-sum test (non-parametric). Frequency and categorical variables were compared using the chi-square test. Correlations were assessed using the Spearman correlation coefficient. A p-value of less than 0.05 was considered statistically significant.

Results

The Role of CRDGs in Cancer Cells and Their Relationship with CRC

In this section, we performed an analysis based on the scRNA-seq data from the GSE178318 dataset. Specifically, this dataset includes six primary CRC samples, six matched liver metastasis samples, and three peripheral blood mononuclear cell samples. After initial processing, including quality control, dimensionality reduction, and clustering, we grouped the samples by type to quantify and compare the number and proportion of cells in different sample types (Figure 1A). Based on gene markers, all cell clusters were categorized into the following known cell lineages: immune cells (T cells, natural killer cells, B cells, macrophages, plasma cells, conventional dendritic cells, monocytes, mast cells, and plasmacytoid dendritic cells), fibroblasts, smooth muscle cells, endothelial cells, and cancer cells. The t-SNE plot shows the distribution of different cell types forming distinct clusters (Figure 1B, left) and their corresponding sample origins (Figure 1B, right). We further subclustered the cancer cell population, identifying 17 clusters (Figure 1C, left), and calculated the proportion of CRC and liver metastasis (LM) samples within each cluster (Figure 1C, right).

Figure 1 The landscape of CRDG expression in cancer cells and its correlation with CRC outcomes. (A) shows the clinical characteristics of CRC samples from the GSE178318 dataset (top), the number of cells in different samples (middle), and the proportion of cells in different samples (bottom). (B) illustrates the distribution of cells in a two-dimensional space, grouped by cell type (left) and by sample (right). (C) depicts the distribution of subclusters within the cancer cell population after re-clustering (left) and the distribution of primary CRC and liver metastasis CRC samples across different clusters (right). (D) presents the heatmap of the gene co-expression network constructed using WGCNA for malignant cells, with deeper colors indicating stronger interactions between modules. (E) shows the protein-protein interaction network among genes in the brown module. (F) illustrates the distribution of cancer cells into high and low CRD score groups (left), along with the proportions of cells in each score group and a summary of the total number of malignant cells (right). (G) provides a bar chart of significantly activated pathways in the high and low CRD score groups. (H) displays box plots showing the differences in pathway scores and CRD scores between high and low expression groups (*** represents a p-value less than 0.001; **** represents a p-value less than 0.0001).

Next, we performed WGCNA using the expression profiles of 2,091 circadian rhythm genes from cancer cells to identify CRC liver metastasis-related genes. We selected brown and blue modules containing 141 CRDGs (Figure 1D) and constructed a protein-protein interaction network for these modules (Figure 1E). Subsequently, we classified cancer cells into high and low expression groups based on the expression levels of CRDGs using a scoring method (Figure 1F, left). The low-scoring cell population was more prevalent in CRC samples, confirming that CRDGs are more highly expressed in LM (Figure 1F, right).

Finally, GSVA was performed on the high and low expression groups, identifying multiple pathways associated with CRC development and progression (Figure 1G and H). Notably, dysregulation of pyrimidine and purine metabolism supports rapid tumor cell proliferation, while alterations in the proteasome pathway assist tumor survival by enhancing protein degradation. Metabolic pathways such as linoleic acid metabolism and nitrogen metabolism contribute to lipid signaling and amino acid synthesis, promoting tumor growth and metastasis. Additionally, pathways related to neuroactive ligand-receptor interaction affect cancer cell proliferation and immune evasion. Pathways involved in base excision repair, focal adhesion, and mismatch repair are crucial for maintaining DNA integrity and cell adhesion; defects in these pathways lead to increased mutation burden and metastasis. Immune-related pathways play key roles in tumor immune surveillance and evasion, further driving CRC progression. The results indicate that CRDGs play a significant role in CRC development and metastasis. High expression of CRDGs is associated with liver metastasis, and these genes are involved in various pathways critical for tumor cell proliferation, survival, and immune evasion. These findings provide insights into the molecular mechanisms underlying CRC progression and suggest potential therapeutic targets for future research.

The Impact of CRD on Response to Immune Checkpoint Inhibitor Therapy

In this study, we investigated the influence of CRD on the response to immune checkpoint inhibitor therapy by analyzing scRNA-seq data from 27 patients with CRC who underwent treatment with immune checkpoint inhibitors. The dataset included patients with pathological complete response (pCR) and non-pathological complete response (non-pCR), sourced from the GSE205506 repository. Following quality control, dimensionality reduction, clustering, and annotation based on marker genes, we identified eight distinct cell types in two-dimensional space, including immune cells (B cells, T cells, myeloid cells), cancer cells, fibroblasts, smooth muscle cells, endothelial cells, and plasma cells (Figure 2A).

Figure 2 Correlation between CRD and pCR treatment outcomes. (A) shows the distribution of scRNA-seq data from 27 patients treated with immune checkpoint inhibitors, including quality control, dimensionality reduction, clustering, and cell type annotation in two-dimensional space. (B) presents the CRD scores for each cell type, with statistical comparisons of the proportions of different cell types in non-pCR and pCR samples (top) and CRD scores for different cell types (bottom). **** represents a p-value less than 0.0001. (C) illustrates UMAP plots of cancer cells (left) and immune cells (right). (D) displays the distribution of cancer cells based on treatment response (left) and CRD scores (right) in two-dimensional space. (E) shows boxplots of CRD score differences between non-pCR and pCR groups. * represents a p-value less than 0.05. (F) depicts the ranking of patients based on cancer cell CRD scores, categorized into CRD high and CRD low groups (top). The volcano plot shows differentially expressed genes between the two cell groups. (G) presents a heatmap of pathway enrichment in malignant cells from the CRD high score group.

CRD scores were calculated for each cell type, revealing significant differences in CRD scores across various cell types (Figure 2B, bottom). Statistical comparisons indicated that the proportions of different cell types varied between non-pCR and pCR samples (Figure 2B, top). UMAP plots were used to visualize the distribution of cancer cells and immune cells separately (Figure 2C). The distribution of cancer cells based on treatment response (non-pCR vs pCR) and CRD scores was depicted in two-dimensional space (Figure 2D). Boxplots of CRD score differences between non-pCR and pCR groups were presented (Figure 2E), demonstrating that non-pCR samples generally exhibited higher CRD scores. Patients were ranked based on cancer cell CRD scores, categorizing them into CRD high and CRD low groups (Figure 2F, top). A volcano plot illustrated the differentially expressed genes between these two cell groups (Figure 2F, bottom). Pathway enrichment analysis in malignant cells from the CRD high score group was performed using a heatmap, highlighting significant biological processes (Figure 2G).

The results indicated that the oxidative phosphorylation pathway was significantly upregulated in the pCR group, suggesting a close association between enhanced metabolic activity and pathological complete response in CRC. In contrast, T helper cell 17 differentiation was significantly downregulated in the pCR group, while other T-cell-related pathways such as T-cell receptor signaling and Toll-like receptor signaling did not show significant changes. Fatty acid metabolism and biosynthesis processes exhibited a slight upward trend, whereas biological processes like DNA replication, circadian rhythm, cellular senescence, and glycolysis showed minimal differences between the two groups. These findings suggest that metabolic pathways, particularly enhanced oxidative phosphorylation, may play a crucial role in pathological response in CRC, while T-cell-mediated immune responses might be suppressed in this context. The differential expression of genes and the distinct biological processes identified between CRD high and CRD low groups provide valuable insights into the mechanisms underlying the response to immune checkpoint inhibitor therapy in CRC patients.

Cellular Interactions and Cytokine Signaling Activity

To elucidate the interactions between cancer cells and T cells, we conducted a comprehensive cellular communication analysis using the CellPhoneDB package. Our analysis focused on chemokines (Figure 3A), co-stimulatory factors (Figure 3B), and co-inhibitory factors (Figure 3C). The size of the circles in these figures represents the p-values, while the color indicates the average expression levels of ligand-receptor pairs. We observed no significant differences in the overall expression levels of chemokines and co-stimulatory signals between high and low CRD score cell populations derived from cancer cells. This suggests that these interactions are relatively stable across different CRD scores.

Figure 3 Cellular Communication Analysis Results. (A) shows interactions between cancer cells and T cells through chemokines. The size of the circles represents the p-values, while the color indicates the average expression levels of LRs (ligand-receptor pairs). (B) illustrates interactions between cancer cells and T cells through co-stimulatory factors. The size of the circles represents the p-values, while the color indicates the average expression levels of LRs. (C) depicts interactions between cancer cells and T cells through co-inhibitory factors. The size of the circles represents the p-values, while the color indicates the average expression levels of LRs. (D) shows interactions between macrophages and T cells through chemokines. The size of the circles represents the p-values, while the color indicates the average expression levels of LRs. (E) illustrates interactions between macrophages and T cells through co-stimulatory factors. The size of the circles represents the p-values, while the color indicates the average expression levels of LRs. (F) depicts interactions between macrophages and T cells through co-inhibitory factors. The size of the circles represents the p-values, while the color indicates the average expression levels of LRs. (G) presents the ligand-receptor interaction network between cancer cells and T cell subtypes in the high CRD score cell population. (H) presents the ligand-receptor interaction network between cancer cells and T cell subtypes in the low CRD score cell population. (I) displays the activity of chemokine pathways in high and low CRD score cancer cells (top) and the interaction networks among CRC samples, macrophages, cancer cells, and T cell populations (bottom). (J) depicts cytokine signaling activity across various cell types using a bubble plot, with the size of bubbles representing the activity levels.

In contrast, co-inhibitory factors such as FAM3C-PDCD1 were significantly elevated in the high CRD score cell population (Figure 3C). This indicates a potential role of co-inhibitory factors in modulating T cell responses in a CRD-dependent manner. Figure 3D–F illustrate the interactions between macrophages and T cells. We found that the interactions between T cell subtypes and macrophages were predominantly mediated by CCL20-CXCR3 (Figure 3D). Additionally, M1 and M2 macrophages exhibited similar co-inhibitory and co-stimulatory signaling effects on T cells in both high and low CRD score cell populations (Figure 3E and F). This suggests that macrophages play a consistent role in modulating T cell responses, regardless of the CRD score.

Further network analysis of high and low CRD score cell populations revealed distinct patterns in the communication processes between cancer cells/macrophages and T cells (Figure 3G and H). These networks highlight the complexity of cellular interactions and the potential for differential signaling pathways in different CRD score populations. Figure 3I displays the activity of chemokine pathways in high and low CRD score cancer cells (top) and the interaction networks among CRC samples, macrophages, cancer cells, and T cell populations (bottom). The chemokine pathway activity was found to be significantly different between high and low CRD score cell populations, indicating a role of chemokines in regulating cellular interactions in a CRD-dependent manner. We also observed differential activation of cytokine signals between high and low CRD score cell populations, including IL3 and IFNL (Figure 3J). This suggests that cytokine signaling pathways are differentially regulated in a CRD-dependent manner, potentially influencing the immune response to cancer. In summary, our analysis reveals significant differences in cellular communication between high and low CRD score cell populations, particularly in co-inhibitory factors and cytokine signaling pathways.

Validation of CRD Scores in Bulk RNA-Seq Data

In the present study, we conducted an analysis of CRC samples derived from the GSE221103 and GSE221173 datasets to ascertain the clinical relevance of copper death-related gene signatures (CRDGs). Figure 4A presents heatmaps and boxplots illustrating CRD scores, which reveal significant disparities between groups with high and low CRD scores. The boxplots demonstrate distinct distributions of CRD scores, suggesting a potential correlation with clinical outcomes. Figure 4B depicts scatter plots that elucidate the relationship between CRD scores and patient prognosis across the two datasets. The data imply that higher CRD scores may be associated with a poorer prognosis, although further statistical analysis is required to substantiate this trend.

Figure 4 Validation of CRDGs in Bulk RNA-seq Data. (A) Heatmaps and boxplots of CRD scores in the GSE221103 and GSE221173 datasets, showing differences in CRD scores between high and low score groups. * indicates a p-value less than 0.05; **** represents a p-value less than 0.0001. (B) Scatter plots depicting the relationship between CRD scores and prognosis in the two datasets. (C) Immune scores and tumor purity of CRC samples with high and low CRD scores from the TCGA database. The left panel shows immune scores, while the right panel shows tumor purity. **** represents a p-value less than 0.0001. (D) Heatmap of the correlation between CRD scores and immune cells. (E) Expression profiles of immune response-related factors (antigen presentation-related factors, co-inhibitory factors, and co-stimulatory factors) in CRC samples.** represents a p-value less than 0.01; **** indicates a p-value less than 0.0001; ns indicates a p-value more than 0.05. F. Enrichment of pathways in CRC samples from the TCGA database, comparing high CRD score (top) and low CRD score (bottom) groups.

To investigate the immunological characteristics of CRC samples, we leveraged the TCGA database to evaluate immune scores and tumor purity in samples exhibiting high and low CRD scores. The left panel of Figure 4C indicates that the immune score was significantly elevated in the group with low CRD scores compared to the group with high CRD scores. In contrast, the right panel suggests that the group with high CRD scores displayed greater tumor purity. Figure 4D presents a heatmap that delineates the correlation between CRD scores and various immune cell types. The data imply a positive correlation between CRD scores and the majority of immune cells, indicating a potential link between CRDGs and immune cell infiltration.

Figure 4E illustrates the expression profiles of immune response-related factors in CRC samples, encompassing antigen presentation-related factors, co-inhibitory factors, and co-stimulatory factors. Notable expression differences were observed between the high and low CRD score groups, further corroborating the association between CRDGs and immune responses. Lastly, Figure 4F depicts the enrichment of pathways in CRC samples from the TCGA database, contrasting groups with high and low CRD scores. Although pathways such as glycolysis and gluconeogenesis exhibited marginally higher enrichment scores in the high CRD score group (NES = 1.16 for glycolysis and 1.32 for gluconeogenesis) compared to the low CRD score group (NES = 0.97 for glycolysis and 1.02 for gluconeogenesis, respectively), these differences did not achieve statistical significance (Adjusted P-value > 0.05). This suggests that while CRD scores may influence the activity of these metabolic pathways, the differences were not statistically significant.

Expression Landscape of CRD Scores Across Pan-Cancer

In this study, we analyzed the expression landscape of CRD scores across 22 cancer types from the TCGA database, encompassing over 7,500 transcriptomic samples. Figure 5A illustrates the distribution of CRD scores in each cancer type. The scatter plots show that the distribution of CRD scores is generally similar across most cancer types, with some variations observed in certain types such as LAML and KIRC. Figure 5B presents a bar plot depicting the distribution of high and low CRD score groups based on the median CRD score across the 22 cancer types. The proportion of patients classified into high and low CRD score groups varies slightly among different cancer types, indicating a differential impact of CRD scores on these cancers. Figure 5C shows the differences in CRD scores between tumor samples and adjacent normal samples across various cancer types. The data points represent the log2 fold change, with a larger size indicating a higher number of samples. The results indicate that CRD scores are consistently higher in tumor samples compared to normal tissues across multiple cancer types, suggesting a potential role of CRD in tumorigenesis. Figure 5D presents the results of a univariate Cox proportional hazards model, showing the association between CRD status and overall survival across different cancer types. The size of the dots represents the number of samples, and the color indicates the hazard ratio. The analysis reveals that high CRD scores are associated with poorer prognosis in adrenocortical carcinoma (ACC) and prostate adenocarcinoma (PRAD) samples, confirming the potential of CRD scores as a prognostic predictor in various cancers. Overall, these findings demonstrate the diverse expression landscape of CRD scores across pan-cancer and their potential as a prognostic biomarker in specific cancer types.

Figure 5 Expression Landscape of CRD Scores Across Pan-Cancer. Panel (A) Scatter plot showing the distribution of CRD scores across 22 cancer types from the TCGA database, with each point representing a tumor sample. Panel (B) Bar plot depicting the distribution of high and low CRD score groups based on the median CRD score across the 22 cancer types. The x-axis represents different cancer types, while the y-axis shows the proportion of samples in each CRD score group. Panel (C) Box plots showing differences in CRD scores between tumor samples and adjacent normal samples across various cancer types in the TCGA database. The y-axis represents CRD scores, while the x-axis shows different cancer types. Panel (D) Forest plot from the univariate Cox proportional hazards model showing the association between CRD status (high vs low) and overall survival across different cancer types. The hazard ratio (HR) and 95% confidence intervals are shown for each cancer type, with the y-axis representing different cancer types and the x-axis showing the HR values.

Construction of the CRDGs Model for Predicting Immunotherapy Response

We evaluated the ability of CRD scores to predict patient prognosis in two independent CRC cohorts (GSE39582 and GSE72970) using Kaplan-Meier survival curves. The results showed that CRD scores could partially predict patient outcomes, with higher CRD scores associated with shorter relapse-free survival (RFS). This indicates that CRD scores have independent prognostic value in CRC patients (Figure 6A). CRD scores were also identified as a valuable tool for assessing responses to bevacizumab and Adoptive Cell Transfer (ACT) therapies. In samples from the GSE19860 and GSE72970 datasets, significant differences in treatment responses were observed between high and low CRD score groups. Patients with higher CRD scores were more likely to have a poor response to these treatments, while those with lower CRD scores were more likely to benefit (Figure 6B).

Figure 6 The Role of CRDGs in Predicting Immunotherapy Responses. Panel (A) Kaplan-Meier survival curves showing the prediction of relapse-free survival (RFS) based on CRD scores for CRC patients from the GSE39582 and GSE72970 datasets. High CRD scores are associated with better outcomes. Panel (B) Boxplots illustrating differences in CRD scores between responders (R) and non-responders (NR) to bevacizumab and Adoptive Cell Transfer (ACT) therapies in CRC samples from the GSE19860 and GSE72970 datasets. * indicates a p-value less than 0.05. Panel (C) Results of the LASSO algorithm for selecting prognostic CRDGs associated with patient outcomes, showing the parameters used for CRDG selection. The displays plot the coefficients and log (λ) values used in the LASSO regression. Panel (D) Comparison of CRDG scores, gene expression profile (GEP) scores, and PD-L1 expression between non-responders (NR) and responders (RD) across four representative cohorts (GSE67501, GSE78220, GSE188249, and PMC7499153) including RCC1, melanoma, ovarian cancer, and RCC2. * indicates a p-value less than 0.05; ** represents a p-value less than 0.01; **** represents a p-value less than 0.0001; ns represents a p-value more than 0.05. Panel (E) ROC curves for CRDGs predicting immunotherapy response in the four representative validation datasets. The curves show sensitivity vs 1-specificity for each cohort. Panel (F) Bar chart summarizing AUC values for CRDGs and comparisons with other immune response prediction features across the four cohorts. The chart highlights the predictive power of CRDGs relative to other features.

To further explore the potential of CRDGs as indicators for immunotherapy, we collected CRDGs score data from various cancer types and applied the LASSO algorithm to select prognostic CRDGs in four representative cohorts: renal cell carcinoma (RCC1), melanoma, ovarian cancer, and renal cell carcinoma (RCC2). The results showed that CRDGs scores outperformed the traditional T-cell inflammation gene expression profile (GEP) scores and PD-L1 expression in predicting immunotherapy responses. In four representative cohorts (GSE67501, GSE78220, GSE188249, and PMC7499153), CRDGs scores were compared with GEP scores and PD-L1 expression levels between non-responders (NR) and responders (RD) (Figure 6C). The findings indicated that CRDGs scores demonstrated superior predictive performance in RCC1, melanoma, and ovarian cancer cohorts, although the performance in RCC2 was relatively weaker. Specifically, CRDGs scores showed higher discriminative ability than GEP scores and PD-L1 expression, especially in melanoma and ovarian cancer cohorts (Figure 6D).

ROC curve analysis further validated the predictive performance of CRDGs scores. The prediction AUC values for CRDGs scores were 0.679 for RCC1, 0.938 for melanoma, 0.83 for ovarian cancer, and 0.532 for RCC2. These AUC values were higher than those of GEP scores and PD-L1 expression, indicating that CRDGs scores have significant advantages in predicting immunotherapy responses (Figure 6E and F). In summary, CRD scores and CRDGs demonstrated good prognostic and treatment response prediction capabilities across multiple cancer types, particularly in melanoma and ovarian cancer. As an emerging indicator for immunotherapy response prediction, CRDGs scores show high clinical application potential and could be further explored in more cancer types in the future.

Discussion

Disruption of circadian rhythms influences the biological characteristics of colorectal cancer (CRC) through multiple pathways, encompassing alterations in gene expression, disruption of physiological functions, modification of the tumor microenvironment, and optimization of chemotherapy efficacy. Research has demonstrated that several clock genes, including Per1, Per2, Cry1, and BMAL1, exhibit circadian rhythmicity in colonic epithelial cells and play crucial roles in orchestrating gastrointestinal functions such as cell proliferation and migration.¹⁶ Furthermore, physiological processes like absorption rhythms and mucosal enzyme activity show a pattern of minimal activity at night, indicating that gut motility and other functions are regulated by circadian rhythms.^17,18 Circadian rhythm genes may impact the development of colorectal cancer by modulating the immune microenvironment. In metastatic colorectal cancer (mCRC), understanding the relationship between clock genes and the immune microenvironment could provide new insights into treatment, particularly for immunotherapy.¹⁹ Chronotherapy, which optimizes drug administration timing according to the patient’s biological clock, has been shown to identify optimal dosing times for patients receiving irinotecan-based treatment for mCRC, reducing toxicity and enhancing efficacy.²⁰ Experimental data also suggest that the network of clock genes impacts the progression of colorectal cancer, providing potential prognostic biomarkers and revealing new molecular targets for therapy; for instance, overexpression of the CRY2 gene is associated with tumor progression and poor prognosis.^21,22 This study is pioneering in elucidating the role of CRDGs in CRC, highlighting their potential as both prognostic biomarkers and therapeutic targets. Our findings reveal distinct expression patterns of CRDGs between primary CRC and liver metastases, with notably higher levels observed in metastatic lesions. These observations imply that CRDGs could critically influence tumor metastasis, possibly by modulating metabolic and immune pathways beneficial to tumor growth and immune evasion.

Firstly, the prognostic model based on CRD scores shows strong predictive capability, particularly with higher CRD scores correlating with poorer responses to immunotherapy. This implies that CRD status could be an important factor in determining the efficacy of immunotherapy in CRC patients. This finding is consistent with previous research indicating that circadian rhythm disruptions can affect immune function and cancer progression. Moreover, gene set variation analysis (GSVA) of high and low CRD score groups identified significant pathways, such as angiogenesis, epithelial-mesenchymal transition (EMT), DNA repair, and the p53 pathway, which have been implicated in CRC development. For instance, Yang et al reviewed the molecular mechanisms, signaling pathways, microenvironment, and therapies related to angiogenesis in CRC.²³ EMT plays a crucial role in tumor progression, metastasis, and resistance to therapy in CRC.²⁴ Fawaz N Al-Shaheri et al found an association between DNA repair gene polymorphisms and CRC risk and treatment outcomes.²⁵ Additionally, WDR43, by binding with RPL11, enhances MDM2-mediated ubiquitination of p53, reducing p53 protein stability and inducing CRC cell proliferation and chemoresistance.²⁶

Secondly, our findings indicate that non-pCR samples exhibit higher CRD scores in cancer cells, suggesting that circadian rhythm might be significantly associated with CRC responses to immunotherapy. Gene set enrichment analysis reveals significant upregulation of oxidative phosphorylation pathways in the pCR group, indicating that enhanced metabolic activity may facilitate pathological complete response. Despite downregulation of T-helper 17 cell differentiation in the pCR group, other T-cell related pathways did not show significant changes, suggesting that T-cell-mediated immune responses might be suppressed in this context. These results imply that circadian rhythm, through modulation of metabolism and immune responses, may play a key role in influencing CRC patients’ sensitivity to immune checkpoint inhibitor therapy. Furthermore, malignant cells in high and low-risk groups participate in several CRC-associated pathways, including the Toll-like receptor signaling pathway, oxidative phosphorylation, and canonical glycolysis. Toll-like receptor signaling is crucial for activating innate and adaptive immune responses and plays a key role in inflammation-induced CRC.²⁷ Lin et al found that PHB2-mediated oxidative phosphorylation promotes CRC cell proliferation and tumorigenesis.²⁸ Wang et al identified the significant role of glycolysis in chemoresistance of TRPC5-induced human CRC cells.²⁹

Thirdly, in cancer treatment, the interaction between cancer cells and T cells is a critical factor. To better understand the impact of circadian rhythm-related genes (CRDGs) on this process, we used the CellPhoneDB software package for cell communication analysis. The results indicate that, despite no significant differences in the overall expression levels of chemokines and co-stimulatory signals between high and low CRD score cancer cell groups, co-inhibitory factors such as FAM3C-PDCD1 are significantly elevated in the high score group. This suggests that CRD may regulate cancer cell immune evasion by enhancing inhibitory signals. In the communication analysis between T cells and macrophages, we observed that interactions between T cell subtypes and macrophages are primarily mediated by CCL20-CXCR3. Further analysis revealed similar co-inhibitory and co-stimulatory signals of M1 and M2 macrophages towards T cells in both high and low CRD score cancer cell groups. This phenomenon suggests that CRD scores do not significantly alter the fundamental communication patterns between macrophages and T cells. However, network analysis identified significant differences in communication patterns between cancer cells/macrophages and T cells under different CRD statuses, indicating that cancer cells with varying CRD states might adopt different immune evasion strategies. Additionally, we assessed the impact of CRD on cytokine signaling patterns at the single-cell level, finding significant differences in pathway activities between high and low score cell groups, particularly with frequent interactions between T cells and macrophages as receptors and ligands. We also observed differential activation of cytokine signaling pathways, such as CCL20-CXCR3 and FAM3C-PDCD1, between high and low CRD score cell groups. Wang et al found that CCL20 derived from CRC cells recruits regulatory T cells through FOXO1/CEBPB/NF-κB signaling to promote chemoresistance.²⁰ Jin et al found that CRC cells co-expressing CXCR3 and CXCR4 exhibited more metastases in lung tissues following intravenous injection in mouse models.³⁰ Additionally, NUTM2A-AS1 promotes CRC malignancy by upregulating FAM3C expression.³¹ These findings suggest that CRD may influence cytokine signaling and thereby modulate the immune microenvironment, impacting CRC responses to immunotherapy.

Fourthly, we explored the dynamic changes in CRD scores during tumor progression. Both neuroblastoma and osteosarcoma datasets showed a trend of decreasing CRD scores over time, suggesting that CRD may undergo dynamic changes during tumor progression. High and low score groups enriched in various pathways, such as glycolysis, PI3K/AKT/mTOR signaling pathways, and pentose phosphate, exhibited significant differences. Shen et al found that m6A-dependent glycolysis promotes CRC progression.³² Targeting EGFR/VEGFR related to PI3K/AKT/mTOR, Wnt/β-catenin, or MAPK inhibitors could provide new therapeutic targets for CRC patients.³³ Additionally, ATP13A2 activates the pentose phosphate pathway and promotes CRC growth through the TFEB-PGD axis.³⁴ Furthermore, the assessment of immune scores and tumor purity in CRC samples indicated that tumors with high CRD scores might have stronger immune evasion capabilities and a larger proportion of tumor cells in the microenvironment. In the comparative analysis of metabolic pathways between high and low score groups, although gluconeogenesis and glycolysis pathways showed slightly higher enrichment scores in the high score group, these differences did not reach statistical significance. This suggests that while CRD scores may have some impact on the activity of these metabolic pathways, the differences are not statistically significant. These results indicate that CRD may influence tumor biology through modulation of immune responses and metabolic activity, but the specific mechanisms require further investigation.

Finally, this study underscores the potential of targeting CRDGs in the development of new therapies for CRC. By modulating CRD-related pathways, there is potential to enhance the efficacy of existing treatments and improve patient outcomes. The consistent results across multiple datasets, particularly the pan-cancer analysis, suggest that the role of CRDGs in cancer may extend beyond CRC, laying the groundwork for future research into their potential as universal cancer biomarkers.

Conclusion

This study comprehensively investigated the roles of CRDGs and CRD in CRC progression, metastasis, and immunotherapy response. Our findings revealed that high CRDG expression is associated with liver metastasis and critical pathways for tumor proliferation and immune evasion. CRD scores effectively distinguished responders from non-responders to immune checkpoint inhibitors, with higher scores correlating with poorer outcomes. Additionally, CRD scores showed superior predictive performance compared to traditional biomarkers in assessing immunotherapy response across multiple cancer types, particularly in melanoma and ovarian cancer. These results highlight the potential of CRDGs and CRD scores as biomarkers for prognosis and immunotherapy response prediction. Future work should focus on validating these findings in larger cohorts and exploring the therapeutic potential of targeting CRDGs and CRD pathways in CRC treatment.