Enrichment score of nucleotide metabolism-related genes in scRNA-seq

Twelve samples were obtained for the study of colorectal cancer heterogeneity and to assess differences between tumor and normal samples using single-cell sequencing data. Following the identification of 2000 highly variable genes, we employed principal component analysis (PCA) for dimensionality reduction, focusing on the top 20 principal components (PCs). Subsequently, 30 clusters were generated, and known marker genes were referenced to annotate cellular subpopulations. Visualization using t-distributed Stochastic Neighbor Embedding (t-SNE) illustrated distinct samples, tissue types, clusters, and annotated cellular subpopulations (Fig. 1A-D). The relative expression of marker genes, calculated as Top5, in each cellular subpopulation was depicted through a heatmap (Fig. 1E). Furthermore, Fig. 1F displays the percentage distribution of different cell types across the 12 samples. Violin plots were utilized to showcase the expression patterns of common marker genes in each cell type (Fig. 1G). Subsequently, we evaluated the level of nucleotide metabolism in scRNA-seq data based on the expression of 882 nucleotide metabolism-related genes. Using five common algorithms (AddModuleScore, UCell, GSVA, AUCell, and singscore), we scored gene sets to assess nucleotide metabolism. Figs. 1H and 1I demonstrate that nucleotide metabolism scores (NMS) were relatively higher in epithelial and myeloid cells. Additionally, we compared NMS across various cell types in tumor and normal samples, revealing that NMS in cells such as myeloid cells, fibroblasts, T/NK cells, and epithelial cells was relatively higher in tumors (Fig. 1J). Furthermore, a differential analysis was conducted for epithelial cells in the tumor and normal tissues, identifying a total of 102 NMRGs with expression differences in both settings (Supplementary Table 3).

Fig. 1

Classification of Cell Subpopulations and Gene Expression Scores Related to Nucleotide Metabolism in Colorectal Cancer. (A-D) t-SNE plots depicting diverse samples, tissue origins, cell clusters, and cell subpopulations, color-coded for clarity. (E) Heatmap illustrating the relative expression of marker genes within eight distinct cell subpopulations. Genes with high expression are represented in red, while those with low expression are displayed in blue. (F) Histogram displaying the distribution of cell types across different samples. (G) Expression patterns of commonly used marker genes for cellular annotation within these cell subpopulations. (H) Bubble plots demonstrating the enrichment scores of nucleotide metabolism-related genes per cell type in colorectal cancer. (I) t-SNE plots illustrate the enrichment scores of nucleotide metabolism-related genes for each cell type, with darker shades of green indicating higher scores. (J) The discrepancy in enrichment scores of nucleotide metabolism-related genes for each cell type between cancer and normal tissues. ns, Not significant; * p< 0.05; ** p< 0.01; *** p< 0.001; **** P< 0.0001

Analysis of cellular interactions in scRNA-seq

Cell trajectory analysis provides valuable insights into cellular differentiation relationships, developmental trajectories, and changes in tumor immune cell dynamics at single-cell resolution. In our study, we utilized the “monocle” R package to determine cell trajectories and pseudotime distributions of epithelial cells in tumor tissues, identifying a total of five cell states in epithelial cells during development, with cluster 5 corresponding to the end state of cell development (Fig. 2A). The heatmap in Fig. 2B illustrates the expression patterns of the top 40 nucleotide metabolism-related genes with the highest differential expression at various stages of epithelial cell development, highlighting, for instance, the predominantly high expression of ACOX1 at the end of epithelial cell development. We classified epithelial cells into nucleotide metabolism (NM) high and NM low groups based on the median value of Nucleotide Metabolism Score (NMS) in all epithelial cells in tumor tissues. Bubble plots were then utilized to visualize the results of signaling pathway activity analysis (Fig. 2C), revealing, for instance, enhanced Notch signaling pathway activity in epithelial cells with high nucleotide metabolism scores, particularly with fibroblasts. From a molecular pathology perspective, the accumulation of DNA mutations, particularly in molecules within the Notch signaling pathway, plays a crucial role in malignancy development (Meurette and Mehlen 2018). The circled graph in Fig. 2D demonstrates the strength of ligand-receptor signaling between different cell types, with epithelial cells exhibiting higher nucleotide metabolism scores showing stronger cellular communication with fibroblasts. Further analysis revealed intensified ligand-receptor pair relationships between high-scoring epithelial cells and fibroblasts, exemplified by the TGFA-EGFR interaction (Fig. 2E). The abnormal activation of the TGF-α/EGFR autocrine loop, observed in many malignant tumors, underscores its close association with tumorigenesis and progression (Tang et al. 2016). Finally, we inferred the existence of ligand-receptor relationships and corresponding transcription factors (TFs) between epithelial cells and fibroblasts with higher nucleotide metabolism scores (Fig. 2F), shedding light on potential regulatory mechanisms underlying cellular interactions in the tumor microenvironment

Fig. 2

Cell Developmental Trajectory Analysis and Cell Communication Analysis. (A) Cell trajectory and pseudo-time analysis for malignant cells. (B) Heatmap illustrating the expression patterns of 40 genes related to nucleotide metabolism that exhibit differential expression during cell development. Low expression is represented in blue, while high expression is depicted in red. (C) Bubble diagram showcasing the activity analysis of signaling pathways across various cell types. (D) Circle diagram visualizing the strength of ligand-receptor interactions between different cell types. (E) Identification of highly ranked ligand-receptor pairs and their associated transcription factors between epithelial cells and fibroblasts. (F) Assessment of ligand-receptor strength between diverse cell types

Characterization of nucleotide metabolism in spatial transcriptome sequencing

We employed SCTransform's approach to correct for spatial sequencing depth and conducted a series of normalization processes, resulting in the identification of 14 distinct cellular subpopulations in space following dimensionality reduction clustering (Fig. 3A). Notably, subpopulations 1, 3, 4, and 11 were predominantly situated in the tumor core of colorectal cancer, as depicted in the original representation of the spatial transcriptome. Bubble plots showcasing the expression patterns of the top 20 nucleotide metabolism-related genes with the largest expression differences are illustrated in Fig. 3B. The metabolic activity of different cell subpopulations was further analyzed using the “scMetabolism” R package. Subpopulations 1, 3, 4, and 11, located in the core region of the tumor, exhibited close associations with purine and pyrimidine metabolic activities (Fig. 3C). This metabolic activity enrichment was primarily observed in the core region of the tumor, as depicted in Figs. 3D and 3E. Subsequently, employing Python's Scanpy and stlearn packages, we conducted cell developmental trajectory analysis on spatially resolved cell subpopulations. Following normalization and clustering of spatial transcriptome data, a total of 11 distinct cell subpopulations were identified (Fig. 3F). Intriguingly, cluster 1, situated in the core region of the tumor, exhibited differentiation towards cluster 8 in the peripheral region of the tumor, as observed in the trajectory analysis (Fig. 3G). Furthermore, we employed the RCTD method to back-convolute annotated cell types from single-cell data to spatial data, inferring the predominant cell types at each spatial location. Epithelial cells with high nucleotide metabolism scores were primarily located in the core region of the tumor, whereas those with low scores were predominantly concentrated in the peripheral region of the tumor (Fig. 3H). Finally, according to MISTy's results, epithelial cells with higher nucleotide metabolism scores showed congruence of clustering and higher correlation of spatial interactions with fibroblasts in the internal space (Fig. 3I, J).

Fig. 3

Characterization of nucleotide metabolism in the spatial transcriptome of CRC. (A) Spatial representation illustrating the identification of 14 clusters through stRNA-seq. (B) Bubble plot displaying the expression levels of genes related to nucleotide metabolism within distinct clusters. Red signifies high expression, while blue indicates low expression. (C) Bubble chart presenting the metabolic intensity across various clusters. (D) Spatial depiction of pyrimidine metabolism intensity. (E) Spatial visualization of purine metabolic intensity. (F) Spatial representation of the 11 cellular clusters identified using Python. (G) Spatial map showcasing the developmental trajectory of clusters 1 through 8. (H) An algorithm is used to identify the predominant distribution of different cell types within the CRC spatial map using RCTD. (I) Extrapolation of spatial clustering of different cell types based on MISTy. (J) Projection of spatial correlations among different cell types based on MISTy

Construction of a prognostic model related to nucleotide metabolism

To leverage the potential of nucleotide metabolism-related gene signatures for clinical decision support, we utilized 102 expression-differentiated NMRGs to develop prognostic models for colorectal cancer. These models were constructed using both high-throughput sequencing data and microarray data. Our methodology involved utilizing a training set consisting of 584 CRC samples with available survival data from the TCGA dataset to construct prognostic risk models. Additionally, we employed 232 and 579 CRC patient samples with survival information from the GSE17538 and GSE39582 cohorts for external validation. Initially, we conducted a Univariate Cox analysis to identify five NMRGs significantly influencing OS in CRC patients (Fig. 4A). To address the risk of overfitting and refine the gene selection for OS prediction, LASSO regression analysis was performed, resulting in the selection of four candidate genes from the initial five (Fig. 4B and C). Subsequently, stepwise multifactorial Cox analysis identified ACOX1, ALDOB, CPT2, and TKT as independent prognostic factors. The risk score was then calculated by summing the expression levels of individual genes, each weighted by their corresponding regression coefficients (Fig. 4D). Patients were stratified into low-risk and high-risk groups using the median score as the cutoff point. Survival analyses conducted for both the TCGA-trained group and the GEO-validated group consistently demonstrated that patients in the high-risk category exhibited poorer OS compared to those in the low-risk group (Fig. 4E-H). Additionally, high-risk patients displayed worse progression-free survival (PFS) (Supplementary Figure 1A). Furthermore, receiver operating characteristic curves illustrated the strong predictive capability of the risk score for OS in the TCGA cohort, as depicted in Fig. 4I. The risk plots provided detailed survival outcomes for individual patients across the TCGA cohort, as well as the GSE17538 and GSE39582 cohorts (Fig. 4J-L). These findings underscore the potential clinical utility of our prognostic models based on nucleotide metabolism-related gene signatures, offering valuable insights for personalized treatment strategies in CRC patients.

Fig. 4

Calculation of Risk Scores Associated with Nucleotide Metabolism and Development of Prognostic Models. (A) Forest plot presenting the five prognostic genes identified through univariate Cox analysis. (B) Profiles of LASSO coefficients. (C) Ten-fold cross-validation for selecting tuning parameters in the LASSO model. (D) Results of multivariate Cox analysis for model genes and their corresponding coefficients. (E-H) Kaplan-Meier survival curves for overall survival (OS) of patients categorized into low-risk and high-risk groups in the TCGA cohort, the complete GEO cohort, the GSE17538 cohort, and the GSE39582 cohort. (I) Area under the curve (AUC) values for risk scores at 1, 3, and 5 years in the TCGA cohort. (J-L) Distribution of scores among low-risk and high-risk groups in the TCGA cohort, the GSE17538 cohort, and the GSE39582 cohort, along with patient survival data

Validation of clinical features and construction of nomograms

Considering the robust correlation observed between our nucleotide metabolism-based risk model and adverse prognosis, we sought to assess the potential of our 4-NMRG signature as an independent prognostic predictor in colorectal cancer patients. In the TCGA cohort, we conducted univariate Cox analysis, revealing that the risk score could serve as an independent prognostic indicator, surpassing other common clinical characteristics (such as age, grade, stage, and histologic type) (Fig. 5A). This trend persisted even after multifactorial analysis, further establishing the risk score as the most reliable independent predictor within the cohort (Fig. 5B). Consistent with these findings, in the GEO external validation cohort, the risk score demonstrated its potential as an independent prognostic indicator for patients (Supplementary Figure 1B, C). To enhance the clinical utility of our risk model and aid clinicians in making informed decisions, we developed a nomogram for predicting 1-, 3-, and 5-year survival rates for CRC patients based on correlations between clinicopathologic features and risk scores (Fig. 5C). The risk score displayed a greater impact on OS prediction, underscoring the superior prognostic potential of our 4-NMRG-based risk model. Calibration curves affirmed the accuracy of the nomogram’s predictions (Fig. 5D), while the area under the curve (AUC) at 3 years significantly outperformed other clinicopathologic features (Fig. 5E). Decisions Curve Analysis curves at 3 years (Fig. 5F) and C-index values (Fig. 5G) consistently demonstrated that our constructed nomograms and risk scores provided the highest net benefit, surpassing traditional models, thus wielding more substantial influence on clinical decision-making. We further illustrated associations between risk groupings, clinical characteristics, and the expression of our four modeled genes using heatmaps. Chi-square tests revealed significant associations between risk groupings and patient Stage and histologic type. Interestingly, all four model genes exhibited higher expression in the low-risk group (Fig. 5H), and the high-risk group displayed more advanced Stage staging (Fig. 5I). These analyses reinforce the reliability of the risk score and nomogram as a clinical predictive scoring system.

Fig. 5

Independent Prognostic Analysis of Risk Scores and Clinicopathologic Factors in the TCGA Cohort. (A, B) Univariate and multivariate Cox regression analyses of clinicopathologic variables and risk scores for overall survival (OS) in the TCGA training cohort. (C) Integrated nomograms combining age, grade, and stage for the prediction of OS at 1, 3, and 5 years in colorectal cancer patients. (D) Calibration curves for the nomograms. (E) Area under the curve (AUC) values for risk scores and clinical characteristics at 3 years in the TCGA cohort. (F) Decision curve analysis (DCA) curves for risk scores, nomogram scores, and other clinical characteristics. (G) Assessment of predictive performance using C-Index for different clinical characteristics, nomogram scores, and risk scores. (H) Heatmap displaying the expression profiles of the four model genes and clinical characteristics associated with subgroups, as determined by the chi-square test. (I) Distribution of clinical stages within various score subgroups. ns, Not significant; * p< 0.05; ** p< 0.01; *** p< 0.001; **** P< 0.0001

Mutational landscape and microsatellite instability

The efficacy of immunotherapy often varies greatly among patients with different tumors, and in addition to factors related to tumor type, pathological stage, and immune infiltration, genetic mutations may also affect the efficacy of immunotherapy. The waterfall plot in Fig. 6A illustrates the somatic mutation spectrum of CRC patients, where the most frequent form of mutation in these genes is missense mutation. Figure 6B includes statistical plots of various mutation classifications, summary plots of base alterations, box plots of the various mutation classifications in the samples, statistical plots of the 10 genes with the highest number of mutations, plots of the number of mutation counts included in each sample, and the sample percentage profile. We examined the distribution of the most commonly mutated genes in CRC in risk score subgroups (Fig. 6C). Patients in the high-risk score subgroup exhibited higher tumor mutation load (TMB) relative to patients in the low-risk score subgroup (Fig. 6D, E). Next, we categorized patients into four groups based on median TMB values and median risk scores (high-TMB+ high-risk score, high-TMB+ low-risk score, low-TMB+ high-risk score, and low-TMB+ low-risk score), and the results showed that patients with high-risk scores and low mutations had the relatively worst OS (Fig. 6F). When the DNA mismatch repair function is aberrant, replication errors occurring in the microsatellites result in microsatellite instability (MSI). Microsatellite instability is categorized into Microsatellite High Instability (MSI-H), Microsatellite Low Instability (MSI-L), and Microsatellite Stable (MSS) based on degree. Compared to the MSI-low group, the MSI-high group in CRC had a greater risk score (Fig. 6G, H).

Fig. 6

Mutational Landscape and Microsatellite Instability in CRC Samples. (A) Overview of the mutation landscape in 542 CRC samples. (B) Detailed breakdown of mutation types, with missense mutations being the most common. Single-nucleotide polymorphisms (SNPs) constituted the majority of mutations, with C>T mutations occurring most frequently. Horizontal histograms present the top 10 mutated genes in CRC. (C) Mutation status and tumor mutation load (TMB) of the 20 genes with the highest mutation frequency across different risk subgroups. (D) Comparison of TMB among different subgroups. (E) Correlation analysis between risk scores and TMB. (F) Survival disparities among four subgroups: H-TMB+ high-risk score, H-TMB+ low-risk score, L-TMB+ high-risk score, and L-TMB+ low-risk score. (G) Differences in risk scores of CRC patients in three subgroups based on microsatellite instability: microsatellite high instability (MSI-H), microsatellite low instability (MSI-L), and microsatellite stable (MSS). (H) Percentage of MSI classifications for patients in high-risk and low-risk groups

Prediction of immune infiltration and biological mechanisms

The tumor microenvironment (TME) is a critical determinant of patient clinical outcomes and therapeutic response. Tumor-infiltrating lymphocytes (TILs), comprising various cell types such as effector, regulatory, and inflammatory cells, engage in complex interactions mediated by cytokines and soluble factors. Additionally, tumor cells themselves release immunosuppressive cytokines, influencing immune cell recruitment within the microenvironment. Consequently, the composition of cells and their interactions with cytokines in the TME collectively shape the anti-tumor immune response (Xie et al. 2022).

In this study, we explored the immune landscape of both high and low-risk score groups using various algorithms (Fig. 7A). To delve deeper into the relationship between risk scores and immune-related functions, we evaluated the enrichment scores of different immune cell subpopulations and functions using the ssGSEA approach. Our findings revealed that the high-risk score group exhibited heightened infiltration scores of immune cells and elevated immune pathway scores (Fig. 7E). We employed the “estimate” method to estimate tumor purity by calculating stromal and immune cell ratios across different risk groups (Fig. 7B). Considering the significant impact of immune checkpoint molecules on tumor immunotherapy, we analyzed the expression of immune checkpoint genes within distinct risk score subgroups, revealing higher expression levels across almost all immune checkpoints in the high-risk score group (Fig. 7C). A heatmap depicted immune checkpoint genes, immune scores, immune cell infiltration, and tumor microenvironment scores across various risk score groups (Fig. 7D).

Fig. 7

Analysis of the immune microenvironment and immune-related functions in different risk-scoring subgroups of the TCGA cohort. (A) Evaluation of variations in immune infiltration across risk score subgroups employing seven different algorithms. (B) Assessment of differences in immune scores and stromal scores calculated via ESTIMATE for distinct risk score subgroups. (C) Examination of variations in immune checkpoint expression within different risk score subgroups. (D) Heatmap displaying distinctions in tumor microenvironment (TME) score, immune checkpoint expression, and immune cell infiltration among diverse risk subgroups. (E) Radar chart depicting variations in immune cell infiltration and immune-related pathways assessed via ssGSEA among patients in different risk groups. (F) Correlation analysis between cancer RNA stemness score (RNAss) and risk score. ns, Not significant; * p< 0.05; ** p< 0.01; *** p< 0.001; **** P< 0.0001

Furthermore, we investigated the correlation between the RNA stemness score (RNAss) and the risk score, uncovering a notable negative correlation (Fig. 7F). This suggests that CRC cells with lower risk scores exhibit more prominent stem cell characteristics and lower levels of cell differentiation. These findings imply that patients in the high-risk score group may have a less favorable prognosis, accompanied by heightened immune activity, likely indicative of an immunosuppressive tumor microenvironment in colorectal cancer. This could potentially result in a reduced response rate to immunotherapy.

Moreover, our risk score signature demonstrated a strong positive correlation with various tumorigenic pathways, including epithelial-mesenchymal transition, angiogenesis, and NF-KB signaling pathways (Fig. 8A). Notably, we observed significant distinctions in nucleotide metabolism-related pathways between the risk groups (Fig. 8B). Differentially expressed genes (DEGs) between the two nucleotide metabolism-related risk subgroups were enriched in hormone metabolism and metabolism-related diseases (Fig. 8C). In terms of the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) terms derived from Gene Set Enrichment Analysis (GSEA), we noted distinct enrichment patterns between the low-risk and high-risk groups. Specifically, the low-risk group exhibited enrichment in nucleotide metabolism and nitrogen metabolism, while the high-risk group showed significant enrichment in JAK-STAT and tumor necrosis factor (TNF) signaling pathways (Fig. 8D).

Fig. 8

Biological characteristics of different risk score groups in the TCGA cohort. (A) MsigDB-based GSVA analysis describing the biological properties of the two nucleotide metabolism-related score groups. (B) Metascape-based enrichment analysis of differentially expressed genes between the two risk-scoring groups. (C) t-SNE plots of both GO and Reactome terms describing the differences in nucleotide metabolic pathway activities in the two risk-scoring groups. (D) . GSEA of GO and KEGG terms for the risk signature. ns, Not significant; * p< 0.05; ** p< 0.01; *** p< 0.001; **** P< 0.0001

Prediction of the effects of immunotherapy and chemotherapy

Tumor immunotherapy, particularly immune checkpoint inhibitors (ICBs), has transformed cancer treatment by activating T-cells, reversing CD8 T-cell depletion, and stimulating immune cells to recognize and eliminate tumor cells. However, the effectiveness of ICB is limited to a subset of tumor patients, with many not experiencing long-term benefits (Xiong et al. 2023). To deepen our understanding of how risk scores influence immunotherapy, we utilized TIDE and IPS scores to evaluate patients with tumors and regional lymph nodes, assessing their potential for an immunocompetent response. This approach aimed to more effectively identify suitable candidates for immunotherapy, recognizing that the efficacy of immunotherapy can vary depending on the level of immune infiltration, often influenced by tumor progression. Expanding on these findings, we explored the feasibility of a prognostic model for predicting the response to immune checkpoint blockade (ICB) in colorectal cancer patients.

The violin plot in Fig. 9A illustrates the relationship between IPS and risk groups, with higher IPS values indicating an improved likelihood of responsiveness to PD-1 and CTLA-4 inhibitors. Notably, individuals in the low-risk group exhibited superior immune responses to immune checkpoint inhibitors, particularly CTLA-4 inhibitors. Since the immune microenvironment modulates ICB responses, we conducted an in-depth analysis of the correlation between risk scores and ICB response characteristics. Our findings revealed a significant positive correlation between the risk score and proteasome and APM_signal while displaying a significant negative correlation with other ICB response attributes. Furthermore, the risk score exhibited substantial and meaningful associations with key stages of the tumor immune cycle, including cancer cell antigen release (step 1), cancer antigen presentation (step 2), initiation and activation (step 3), and immune cell infiltration into the tumor (step 4) (Fig. 9B).

Fig. 9

Prediction of the effects of immunotherapy and chemotherapy. (A) Comparison of the relative distributions of immunization scores (IPS) in the high-risk scoring group and the low-risk scoring group. (B) The relationship between risk scores, ICB response traits, and the various stages of the tumor-immunity cycle. (C) Heatmap of modeled gene-immunity gene correlations. (D) Differences in TIDE between CRC patients in the high-risk scoring group and those in the low-risk scoring group. (E) Correlation of risk scores with IC50 values for cisplatin, imatinib, and doxorubicin. ns, Not significant; * p< 0.05; ** p< 0.01; *** p< 0.001; **** P< 0.0001

To further investigate variations in immune responses among different subgroups, we conducted correlation analyses involving 12 model genes and classical immune-related genes (Fig. 9C). Higher tumor TIDE prediction scores are associated with reduced responsiveness to immune checkpoint blockade (ICB) and diminished patient survival. Our findings revealed that individuals in the high-risk score group demonstrated elevated dysfunction and exclusion scores, along with relatively higher TIDE scores, as depicted in Fig. 9D. Lastly, we explored the relationship between risk scores and the IC50 values of three clinically utilized chemotherapeutic agents. Our findings showed a significant negative correlation between risk scores and the IC50 values of cisplatin, imatinib, and doxorubicin (Fig. 9E). Taken together, these results suggest that individuals in the low-risk group may be more likely to benefit from both immunotherapy and chemotherapy.

ACOX1+ and CPT2+ tumor cells may serve as prognostic influencers and targets for immunotherapy

A total of four nucleotide metabolism-related genes (ACOX1, ALDOB, CPT2, and TKT) were incorporated into our risk model. The spatial maps illustrate the expression patterns of these four genes (Fig. 10A, F, Supplementary Figure 2A, C). In our scRNA-seq analysis, we categorized epithelial cells into two groups: those expressing the four genes (expression-positive) and those not expressing them (expression-negative). Utilizing the ssGSEA algorithm, we estimated the abundance of these cells in the TCGA dataset based on marker genes specific to the positively expressing epithelial cells. For survival analysis, we determined the optimal cutoff value and subsequently divided the CRC patients from the TCGA dataset into two groups.

Fig. 10

ACOX1 and CPT2 are protective genes in CRC patients. (A) Spatial map demonstrating the expression of ACOX1 in colorectal cancer. (B) Kaplan-Meier survival curves of OS for patients in the ACOX1+ epithelial cells high and low expression groups. The proportion of ACOX1+ epithelial cells in patients producing different immune responses. (C) Spatial maps of different cell types were obtained by the algorithm of reverse convolution. Included here are ACOX1 expression-positive and expression-negative epithelial cells, endothelial cells, myeloid cells, mast cells, fibroblasts, and T/NK cells. (D, E) Heatmaps and network diagrams to predict the strength of communication between different cell types based on the stlearn method. (F) Spatial map demonstrating the expression of CPT2 in colorectal cancer. (G) Kaplan-Meier survival curves for OS in patients in the CPT2+ epithelial cells high and low expression groups. The proportion of CPT2+ epithelial cells in patients who produced different immune responses. (H) Spatial maps of different cell types were obtained by an algorithm of reverse convolution. Included here are CPT2 expressing positive and expressing negative epithelial cells, fibroblasts, T/NK cells, endothelial cells, myeloid cells, and B cells. (I, J) Heatmaps and network diagrams of the strength of communication between different cell types were extrapolated according to the method of stlearn

Patients with higher proportions of ACOX1+ and CPT2+ epithelial cells exhibited relatively more favorable prognosis and a higher likelihood of responding effectively to immunotherapy (Fig. 10B, G). Conversely, patients with higher proportions of TKT+ epithelial cells had a relatively better prognosis. However, the difference in cell proportions between patients who did or did not produce an effective immune response, as indicated by TIDE analysis, was not statistically significant (Supplementary Figure 2B). Although patients with a higher proportion of ALDOB+ epithelial cells also demonstrated a relatively better overall survival, the log-rank test results showed no significant difference between the two groups (Supplementary Figure 2D). These findings suggest that ACOX1 and CPT2 expression-positive epithelial cells might serve as protective factors for colorectal cancer patients.

The Reverse Cell Type Delineation (RCTD) method was employed to back-convolute well-annotated cell types from the single-cell data into spatial data (Fig. 10C, H). Interestingly, the extrapolation of cell communication relationships in space revealed that ACOX1-negative expressing epithelial cells exhibited stronger cellular interactions with fibroblasts compared to ACOX1-expressing positive epithelial cells (Fig. 10D, E). Similarly, CPT2-negative expressing epithelial cells also displayed enhanced cellular communication with fibroblasts (Fig. 10I, J).

We evaluated the expression of ACOX1 and CPT2 in 80 pairs of colorectal cancer tissues and adjacent normal tissues using IHC methods. Our results showed that ACOX1 and CPT2 were significantly downregulated in colorectal cancer tissues compared with normal tissues, as shown in Fig. 11A. Finally, we compared the expression levels of ACOX1 and CPT2 in normal intestinal epithelial cells and four CRC cell lines by PCR assay and found that the expression levels of ACOX1 and CPT2 genes were significantly down-regulated in tumor cells (Fig. 11B). These results strongly support the potential of ACOX1 and CPT2 as biomarkers for CRC diagnosis and prognosis. Subsequently, the expression levels of ACOX1 and CPT2 were evaluated after 5 days of transfection using qRT-PCR to validate the effect of siRNA knockdown of ACOX1 and CPT2 in RKO and HCT116 cell lines (Fig. 11C). Based on the knockdown efficiency, we chose the No. 1 and No. 2 siRNA knockdown cell lines for ACOX1-related functional experiments, whereas the No. 2 and No. 3 siRNA knockdown cell lines for CPT2-related functional experiments.

Fig. 11

Expression of ACOX1 and CPT2 in CRC. (A) Immunohistochemical staining results showed the protein expression levels of ACOX1 and CPT2 in colorectal cancer tissues. (B) Compared with human intestinal epithelial NCM cell lines, ACOX1 and CPT2 were expressed at lower levels in CRC cell lines. (C) Relative expression of ACOX1 and CPT2 in CRC cells transfected with si-RNA or negative control (NC) was detected by RT-qPCR. ns, Not significant; * p< 0.05; ** p< 0.01; *** p< 0.001; **** P< 0.0001

Subsequently, CCK-8 cell experiments showed that the knockdown-induced reduction of ACOX1 and CPT2 significantly enhanced the proliferation of RKO and HCT116 cell lines (Fig. 12B). Tumor cells transfected with si-ACOX1 and si-CPT2 also exhibited enhanced migration and invasion in transwell assays (Fig. 12A). Together, these findings suggest that ACOX1 and CPT2 are oncogenes for colorectal cancer.

Fig. 12

Functional experiments of ACOX1 and CPT2 in CRC. (A) Transwell assay showed that down-regulation of ACOX1 and CPT2 expression promoted the migration and invasion ability of CRC cells. (B) CCK8 assay showed that the proliferation ability of CRC cells with reduced expression of ACOX1 and CPT2 was significantly enhanced compared with the NC group. ns, Not significant; * p< 0.05; ** p< 0.01; *** p< 0.001; **** P< 0.0001