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

Abstract

Background: Colorectal cancer (CRC) is one of the most prevalent malignancies worldwide. With the increasing incidence of CRC, there is an urgent need for effective strategies for early diagnosis and treatment. Circadian rhythm, a natural biological clock, regulates various physiological processes, and its disruption has been implicated in the onset and progression of cancer. However, the specific roles of circadian rhythm-related genes (CRDGs) in CRC remain unclear.

Methods: In this study, we analyzed the expression patterns of CRDGs in CRC using single-cell RNA sequencing (scRNA-seq) data and bulk RNA sequencing data from the GSE178318 dataset. We constructed a CRC prognostic model based on CRD scores. Additionally, we explored the potential mechanisms of CRDGs in tumor progression through weighted gene co-expression network analysis (WGCNA) and gene set enrichment analysis (GSEA), and assessed their impact on the response to immune checkpoint inhibitors.

Results: The analysis revealed that CRDGs were significantly upregulated in liver metastasis samples compared to primary CRC samples and were closely associated with several metabolic and immune-related pathways. The prognostic model based on CRD scores indicated that higher CRD scores were associated with poorer outcomes in immunotherapy. These findings were further validated in multiple datasets, underscoring the potential of CRDGs as prognostic indicators in CRC.

Conclusion: This study systematically reveals, for the first time, the expression characteristics of CRDGs in CRC and their relationship with tumor progression and response to immunotherapy. CRDGs may serve as effective prognostic biomarkers and therapeutic targets, offering new strategies for the personalized treatment of CRC.

Keywords: circadian rhythm-related genes; colorectal cancer; immune checkpoint inhibitors; prognostic model; single-cell RNA sequencing.

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).

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.

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.

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.

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.

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.