Predicting prognosis and immunotherapy response of ferroptosis-related genes in colorectal cancer

Abstract

See More

To identify ferroptosis markers, which were associated with prognosis and immune response in colorectal cancer (CRC), we integrated spatial transcriptome, single-cell transcriptome, and bulk transcriptome data for analysis, and this multi-omics integration approach is an innovation that distinguishes this study from previous ones. We assessed the prognostic value of ferroptosis-related genes by Kaplan–Meier survival analysis, subject operating characteristic (ROC) curves, and Lasso analysis. A total of 135 ferroptosis-associated genes were identified, with the high-risk subgroup having significantly worse overall survival. The prognostic model showed strong predictive power. Immunohistochemistry (IHC) assay showed that Sorting Nexin 5 (SNX5) and Ubiquitin-specific protease 7 (USP7) were significantly up-regulated in CRC tissues compared with adjacent CRC tissues. Immune infiltration analysis revealed higher infiltration levels of CD8 + T cells, regulatory T cells, and resting NK cells in the high-risk group, as well as higher dysfunction scores, suggesting possible immunotherapy tolerance. Functional enrichment analysis revealed that the highly expressed genes were mainly enriched in extracellular matrix, integrin binding, and cell adhesion-related pathways. This study provides a comprehensive multi-omics framework for identifying prognostically relevant ferroptosis gene markers and immune response predictors in CRC, providing an important research foundation for precision oncology.

Introduction

Colorectal cancer (CRC) is currently the leading cause of death in men, and ranks second as the cause of death in women. The incidence of colon cancer among individuals younger than 55 years is increasing by 1–2% per year [1]. Tumor immune escape and tumor microenvironment (TME) formation are important factors in CRC development and progression. Although the emergence of immunotherapy has brought hope to patients with CRC, the majority of patients do not benefit from it. Thus, it is crucial to prioritize the assessment of tumor spatial characteristics and identify cell types that exhibit enhanced interactions with tumor cells. Identifying markers to distinguish between cell types in different regions can enable personalized treatment and lead to better patient outcomes. Recent advances in spatial transcriptomics and single-cell sequencing have enabled the differentiation of immune environments within tumor and mesenchymal regions. This approach, when combined with multi-omics, facilitates the identification of markers to predict the immune environment status of various patients with CRC.

Ferroptosis is an iron-dependent, regulated cell death characterized by excessive accumulation of iron, lipid peroxides, and reactive oxygen species (ROS). However, its role in CRC development remains a topic of active research. Key genes associated with ferroptosis have proven effective in predicting CRC prognosis [2, 3], with low ferroptosis scores correlating with better progression-free survival and chemotherapy response in patients with CRC [4]. Several studies have reported that ferroptosis affects the immune microenvironment of CRC. For example, cetuximab, an approved treatment for metastatic CRC that targets EGFR, promotes ferroptosis by regulating the p38/Nrf2/HO-1 axis in KRAS-mutant CRC cell lines [5]. In addition, interferon (IFN)-ɣ promotes tumor cell lipid peroxidation and ferroptosis by regulating the expression of solute carrier (family-3 member 2 [SLC3A2 and family-7 member 11 [SLC7A11]), two subunits of the glutamate-cysteine antiporter system xc- in CRC [6, 7]. Further, IFN-ɣ derived from CD8⁺T cells induces immunogenic tumor ferroptosis by stimulating acyl-COA synthetase long-chain family member 4 (ACSL4) [8] or downregulating SLC3A2 and SLC7A11 [9]. Ferroptosis participates in CRC pathogenesis via various key pathways involved in immune regulation [10]. Therefore, identifying relevant ferroptosis-related genes could aid in predicting responses to CRC immunotherapy.

Currently, most relevant studies mainly rely on bulk RNA sequencing data, which makes it difficult to distinguish between tumor- and mesenchymal region-specific gene expression features, thus limiting accurate prediction of prognosis and immunotherapy response. Although single-cell sequencing technology can reveal cell type composition, it disrupts the spatial structure of tissues; whereas spatial transcriptome technology retains tissue structural information but with limited resolution. To compensate for these shortcomings, this study integrated spatial transcriptome, single-cell transcriptome and bulk transcriptome data to systematically characterize ferroptosis genes associated with CRC prognosis and immune response. The multi-omics and multi-dimensional analysis strategy enabled us to precisely resolve the expression patterns of ferroptosis genes under different spatial and cellular heterogeneity in colorectal cancer, providing new insights into the prediction of patient outcome and immunotherapy efficacy, and effectively filling an important knowledge gap in the field. Our findings provide a comprehensive theoretical foundation for future precision oncology research and clinical applications.

Methods

Data resources

Spatial transcriptome data from a human patient with colorectal cancer (CRC)were downloaded from the official website of 10X Genomics (https://www.10xgenomics.com/cn/resources/datasets/human-colorectal-cancer-whole-transcriptome-analysis-1-standard-1−2–0). All data were analyzed using 10X Visium spatial transcriptome technology with an Illumina NovaSeq 6000 sequencing instrument. The sample captured 3,138 spots with a median gene count of 3,538 per spot and a median UMI per spot of 8,906 spots.

Gene expression sequencing results in transcripts per million (TPM) format for 647 primary tumor tissues from patients with CRC, as well as patient clinical data, including sex, age, stage, grading, and survival prognosis, were obtained from the official site of TCGA GDC (https://portal.gdc.cancer.gov/). Additionally, for all the patients included in this study, we downloaded the “Masked Copy Number Segment"of the copy number variation (CNV) as somatic cell copy number data, and the"Masked Somatic Mutation” of the single nucleotide variation (SNV) as somatic mutation data.

Gene expression data of 441 patients with CRC and clinical survival data of patients from three datasets, GSE17536 [11], GSE17537 [11], and GSE87211 [12], were downloaded from the GEO database, and the microarray platforms were based on GPL570 ([HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array) and GPL13497 (Agilent-026652 Whole Human Genome Microarray 4 × 44 K v2 (Probe Name version)). The data were subsequently homogenized, and batch effects were removed using R-based limma [13] and sva [14] packages.

The ferroptosis-related gene set was obtained from the FerrDB database, the official website of which is http://www.zhounan.org/ferrdb[15]. Enrichment analysis was performed using Metascape [16] and the profile for enrichment analysis was based on a gene ontology (GO) [17] gene set.

ESTIMATE analysis

Spatial transcriptome analysis was conducted using the Seurat R package [18] to identify spots with different transcriptome features using dimensionality reduction clustering. Subsequently, ESTIMATE analysis was performed using the ESTIMATE R package [19] to identify tumor and mesenchymal regions by calculating the tumor purity and immune and mesenchymal scores. The criteria for tumor regions were an average tumor purity ≥ 8, an immune score ≤ 350, and a mesenchymal score ≤ 100, and vice versa for mesenchymal regions.

MCPCOUNTER analysis

Using MCPCOUNTER analysis [20], the infiltration levels of different immune cells (including T cells, B cells, NK cells, monocyte lineages, neutrophils, myeloid-derived cells, fibroblasts, and endothelial cells) in each spot of the spatial transcriptome were identified, and the spatial immune infiltration distribution information of CRC was subsequently obtained.

CRC single-cell data processing

First, we obtained CRC single-cell transcriptome data from GSE146771 [21] using the TISCH database [22] and differentiated the single-cell data into malignant, immune, and mesenchymal cell datasets using annotation information. Subsequently, clustering and cellular annotations were performed on Dimplot after Uniform Manifold Approximation and Projection (UMAP) downscaling of single-cell data, while heatmaps identified the highly expressed characteristic genes of the three cell types.

Identification of ferroptosis genes characterizing CRC at multiple levels

In terms of feature selection and downscaling techniques, we performed differential expression analysis on different data types (spatial transcriptome, single-cell, and large-scale RNA sequencing) for feature selection, with a special focus on genes associated with ferroptosis, using stringent criteria (p-value < 0.05) to refine the gene set. First, differential expression analysis was performed on the tumor and mesenchymal regions of the spatial transcriptome using Seurat's FindAllMarker function. Subsequently, the intersection of the differentially expressed genes and ferroptosis gene sets was identified in the spatial transcriptome to obtain the spatial level of differentially expressed ferroptosis genes in CRC. Subsequently, differential expression analysis of malignant, immune, and mesenchymal cells in the single-cell transcriptome was performed using Seurat FindAllMarker function. The intersection of differentially expressed genes and ferroptosis gene sets obtained from tumor cells in the single-cell transcriptome was determined to obtain CRC single-cell-level differentially expressed ferroptosis genes. Finally, differential gene expression analysis was performed on the CRC and control groups in TCGA dataset using the limma R package [13]. The intersection of the differentially expressed genes and the ferroptosis gene set of the tumor samples in the bulk transcriptome was determined to obtain the CRC bulk-level differentially expressed ferroptosis genes. In the three differential expression analyses, a p value < 0.05 was used to define differentially expressed genes. Ferroptosis-related genes exhibiting differential expression in all three sets of differential expression analyses were combined to obtain multilevel CRC-characteristic ferroptosis genes.

Construction of CRC prognostic models

Based on the multilevel CRC-characteristic ferroptosis gene set, we first screened the genes using one-way Cox analysis to identify the CRC prognosis-related ferroptosis genes. Subsequently, we constructed a prognostic model of CRC related to ferroptosis using least absolute shrinkage and selection operator (LASSO) analysis and Cox stepwise regression analysis using the TCGA colorectal cancer cohort as the training cohort and the GEO colorectal cancer cohort as the test group. By summing the product of each modeled gene multiplied by its coefficients, we obtained the risk score for each sample and categorized the patients into low- and high-risk CRC subgroups based on the median value of the risk score. The predictive validity of the model was tested by performing K-M survival analysis and time-dependent ROC curve analysis on both cohorts.

Construction of a clinical phenotype nomogram

To better validate the model’s independent prognostic capability, we combined the risk score and common clinical phenotypes of colorectal cancer patients and used univariate and multivariate Cox analyses to remove covariates and further demonstrate the predictive validity of the prognostic model. Subsequently, we constructed a clinical nomogram and calibration curve (calibration) by combining the risk score and clinical phenotypes with prognostic values in a multifactorial Cox regression, which allowed us to derive nomogram scores and predict the prognosis of patients with CRC by weighting the scores of clinical phenotypes and risk score.

Functional annotation enrichment analysis

First, we analyzed the differential expression of genes in the low- and high-risk subgroups of CRC using the limma package [13], defining differentially expressed genes as those with p-values < 0.05, and |log2(Fold Change [FC])|> 1. Using the Metascape online tool [14], we performed GO enrichment analysis on the differentially expressed genes between the two subgroups to investigate the differential functional characteristics between the different disease subtypes using GO analysis of biological process (BP), cellular component (CC), and molecular function (MF) gene subsets. The enriched functionally related pathways were determined according to -log10(p) values in descending order (p-values in ascending order).

Immune cell infiltration and immune microenvironment characterization

CIBERSORT (https://cibersort.stanford.edu/) is a valuable tool for deconvoluting the composition of complex tissues, particularly in the context of immune cell subtypes, based on gene expression data [23]. This method is based on a known reference set that provides a set of gene expression features for 22 immune cell subtypes to calculate the immune cell infiltration. The correlation between immune cells has an impact on immune pathways and functions of the immune system. Therefore, after obtaining the immune infiltration information, we compared the scores of different immune cells between the two subgroups and obtained immune cells with different infiltration levels between the two groups using the Wilcoxon test. Subsequently, we calculated the correlation coefficients between immune cells using Spearman’s analysis and constructed a correlation heat map to visualize the correlations between immune cells.

Immunotherapy correlation analysis

Tumor mutation burden (TMB) and microsatellite instability (MSI) play significant roles in affecting the interaction between immune cells and tumor cells during immunotherapy, and a large number of studies have illustrated their importance during immunotherapy; therefore, we compared the TMB and MSI scores in the two subgroups. T cell dysfunction and exclusion (TIDE) [24] is an algorithm used to predict the functional status of T cells in tumors. We calculated the dysfunction and exclusion scores for each sample and compared their distributions within the two subgroups.

Construction of risk gene-related protein–protein interaction (PPI) networks and various regulatory networks

In this study, based on CRC ferroptosis-related prognostic model risk genes, we used the Search Tool for Retrieving Interacting Genes (STRING) database [25] to construct protein–protein interaction (PPI) networks between genes. Based on risk genes, the miRTarBase [26] database (https://mirtarbase.cuhk.edu.cn/php/download.php) was employed for forecasting information about miRNA-mRNA interactions of risk genes; the chEA database [27] (https://maayanlab.cloud/chea3/) to predict information about transcription factors (TF) of risk genes interacting with genes; using the Comparative Toxicogenomics Database (CTG) [28] database (http://ctdbase.org/) to predict risk gene-drug interaction information to explore drugs that target risk genes. Subsequently, Cytoscape was used to visualize the risk gene-related PPIs, mRNA-miRNAs, gene-TFs, and risk gene-drug molecule interaction networks.

Immunochemistry (IHC)

We determined USP7 and SNX5 expression in five human colon cancer tissues and five adjacent colon cancer tissues using IHC. After dewaxing, the slices were put into EDTA (pH9.0) solution in microwave for antigen repair. Tissue sections were covered evenly with 3% BSA and closed at room temperature for 30 min. Then, rabbit anti-SNX5 (1:200) (17,918–1-AP, Proteintech, China) or anti-USP7 (1:200) (BSM-54298R, BIOSS, China) was added dropwise to the sections and incubated overnight at 4 °C. After that, the sections were shaken dry and covered with HRP-labelled goat anti-rabbit secondary antibody in a circle and incubated at room temperature for 50 min. The sections were covered with freshly prepared DAB chromogenic solution in a circle, controlled under the microscope, and the color was developed. The color development time was controlled under the microscope, the positive color was brownish yellow, and the color development was terminated by rinsing the section with tap water. Hematoxylin re-staining for about 3 min, tap water wash, hematoxylin return blue solution return blue, rinse with running water. After dehydration and sealed, the slides were placed under a white light microscope for interpretation of the results.

Statistical analysis

R software (version 4.1.1) was used for data processing and analysis. Differences between non-normally distributed variables were analyzed using the Mann–Whitney U test (Wilcoxon rank-sum test). Correlation coefficients between different non-normally distributed data were calculated using the Spearman’s correlation analysis. Differences in patient survival between the two groups were analyzed using the log-rank test for K-M survival analysis. IHC scoring was performed by three independent pathologists and quantitatively analyzed based on the intensity of cell staining and the proportion of positive cells. During the analysis, a semi-quantitative scoring system was used to assess the level of U2AF1 expression. The scoring criteria included the proportion of positive cells and staining intensity, which were combined for a composite score. Student-t test was used to compare the expression of SNX5 or USP7 in two groups. P < 0.05 was regarded as statistically significant.

Results

The overall analytical flow of the study

The overall analytical flowchart of this study is presented in Fig. 1.

Fig. 1

The overall analytical flowchart of this study

Identification of tumor and mesenchymal regions in the CRC spatial transcriptome

Using principal component analysis (PCA) downscaling and clustering, we categorized the spots of the spatial transcriptome samples into 16 clusters (0–15). Spots within each cluster were spatially proximate (Fig. 2A, B). Subsequently, using ESTIMATE analysis, we identified the different clusters as tumor and mesenchymal ecotopes (Fig. 2C). The tumor purity, immune, and mesenchyme scores in ESTIMATE analysis are shown in Fig. 2D–F, respectively, to identify tumor and mesenchymal regions with the criteria of average tumor purity ≥ 8, immune score ≤ 350, and mesenchymal score ≤ 100, which is in good agreement with the pathology results.

Fig. 2

Identification of tumor and mesenchymal regions of spatial transcriptome samples. A Pathological sectioning of spatial transcriptome samples; B Clustering of spatial transcriptome samples; C Delineation of tumor and mesenchymal regions of spatial transcriptome samples; D Tumor purity scoring by ESTIMATE analysis; E Immunity score scoring by ESTIMATE analysis; F Mesenchymal scoring by ESTIMATE analysis

Distribution of immune infiltration based on spatial transcriptome

The infiltration levels of T cells, B cells, NK cells, monocyte lineages, neutrophils, myeloid-derived cells, fibroblasts, and endothelial cells were identified by MCPCOUNTER analysis. T cells had a scattered distribution (Fig. 3A), B cells were predominantly located in the mesenchymal region (Fig. 3B), NK cells had a lower level of infiltration (Fig. 3C), monocyte lineages and neutrophils showed scattered distributions (Fig. 3D, E), cells of myeloid origin had a low level of infiltration (Fig. 3F), and fibroblasts and endothelial cells were mainly distributed in the mesenchymal area (Fig. 3G, H). After statistical comparison, we found that B cells, monocytes, myeloid cells, endothelial cells, and fibroblasts showed high levels of infiltration into the mesenchymal region, whereas NK cells and neutrophils were highly distributed in the tumor area (Fig. 3I).

Fig. 3

Immune infiltration analysis of spatial transcriptome samples. A: T-cell infiltration in the spatial transcriptome samples; B: B-cell infiltration in the spatial transcriptome samples; C: NK-cell infiltration in the spatial transcriptome samples; D: monocyte lineage cell infiltration in the spatial transcriptome samples; E neutrophil infiltration in the spatial transcriptome samples; F myeloid-derived cell infiltration in the spatial transcriptome samples; G fibroblast infiltration in the spatial transcriptome samples; H endothelial cell infiltration in the spatial transcriptome samples. I Comparison of the eight cellular immune infiltration scores in the tumor and mesenchymal areas of the samples

Annotation of CRC single-cell data

First, we performed clustering (Supplemental Figure S1A) and cellular annotation (Supplemental Figure S1B) of the CRC single-cell transcriptome data of GSE146771 using UMAP plots. Single-cell data were grouped into 20 clusters (clusters 0–19) and annotated as immune, malignant, and mesenchymal cells. The three clusters were analyzed for differential expression. Differentially expressed genes were identified as those with a p-value < 0.05, and the top 10 characterized genes of the three cell lines are shown in Supplemental Figure S1C.

Identification of multilevel CRC-characteristic ferroptosis genes

Through the integrated multi-omics analysis, we screened a total of 135 multilevel characteristic ferroptosis-related genes, which showed significant differential expression in the spatial transcriptome, single-cell transcriptome and bulk transcriptome (Supplementary Figure S2, Table 1). These key genes laid the foundation for subsequent prognostic model construction and immune-related analysis.

Table 1 Part of 135 multilevel CRC genes obtained

SNX5 and USP7 was up-regulated in CRC tissues

To verify the expression of Sorting Nexin 5 (SNX5) and Ubiquitin-specific protease 7 (USP7), as two of the above ferroptosis genes, we performed immunochemistry (IHC), and the result showed SNX5 and USP7 were significantly high-expressed in CRC tissues compared to adjacent CRC tissues (Fig. 4A and B) (P = 0.03; P = 0.00, respectively).

Fig. 4

Expression of SNX5 and USP7 in CRC tissues via IHC. A High expression of SNX5 was observed in CRC tissues (200x). B USP7 was up-regulated in CRC tissues compared to adjacent tissues(200x)

Construction and validation of ferroptosis-related prognostic model and clinical nomogram for CRC

To forecast patient survival outcomes according to the expression and differences in multilevel CRC-characteristic ferroptosis genes, we screened prognosis-related genes using one-way Cox analysis and built a prognostic model using the lasso-Cox stepwise regression method. The corresponding coefficients of the best and simplest LASSO regression models are shown in Fig. 5A. Figure 5B shows the relationship between the LASSO regression entry characteristics and the absolute values of the coefficients. The models were constructed based on the training dataset (TCGA CRC cohort), and validated using the test dataset (GEO CRC cohort). The risk scores (Risk score) of each sample in the TCGA and GEO cohorts and information on survival outcomes are shown in Fig. 5C and E. Patients were categorized into high- and low-risk subgroups based on the median risk score. We found that the high-risk subgroup had significantly worse prognosis than the low-risk subgroup in both cohorts (Fig. 5D, F).

Fig. 5

Construction and validation of lasso-Cox stepwise regression prognostic models. A Obtaining the best and simplest models for LASSO regression; B The relationship between LASSO regression entry characteristics and absolute values of coefficients; C Presentation of model scores and survival outcomes for each sample of the TCGA cohort; D Survival differences between the high- and low-risk subgroups in the TCGA cohort; E Presentation of model scores and survival outcomes for each sample of the GEO cohort; and F Survival differences between the high-risk and low-risk subgroups in the GEO cohort

The ROC curves of both cohorts had high area under the curve (AUC) values, with 1-, 3-, and 5-year values of 0.702, 0.728, and 0.724 in the TCGA cohort (Supplemental_Fig_S3A) and 0.689, 0.661, and 0.660 in the GEO cohort, respectively (Supplemental Figure S3B), suggesting a high prognostic predictive ability. We subsequently validated the independent prognostic power of the risk score using univariate and multivariate Cox analyses after removing other factor covariates (Supplemental Figure S3C, D). Subsequently, we established a clinical nomogram using risk scores, age, stage, and T, M, and N staging to more precisely predict the prognosis of patients (SupplementalFig. S3E) and established 1-, 3-, and 5-year calibration curves (SupplementalFig. S3F), indicating an excellent prognostic predictive ability.

Functional enrichment analysis of the two subgroups

To investigate the functional pathway differences between the two CRC subgroups, we performed a GO enrichment analysis of the upregulated and downregulated differentially expressed genes. The upregulated differentially expressed genes in the GO enrichment analysis were significantly enriched in the extracellular matrix, integrin binding, and positive regulation of cell adhesion pathways (Supplemental Figure S4A). The downregulated differentially expressed genes in the GO enrichment analysis were significantly enriched in ribosome biogenesis, nucleobase-containing small molecule metabolic, and peptide metabolic processes (Supplemental Figure S4B).

Immune cell infiltration and its correlations and differences

Gene expression data were analyzed for the types of infiltrating immune cells using CIBERSORT, and a panorama of 22 infiltrating immune cell types in each sample was obtained by removing cells with scores of zero in all samples (Fig. 6A). After comparing the immune infiltration scores of the two subgroups (Fig. 6B), the infiltration of immune cells, such as plasma cells, resting memory CD4 T cells, and activated memory CD4 T cells, was significantly higher in the low group. In the high group, the infiltration of immune cells, such as CD8 T cells, regulatory T cells (Tregs), and resting NK cells, was significantly higher. We then analyzed and visualized the correlation between immune cells (Fig. 6C) and found that follicular helper T cells had a high positive correlation with CD8 T cells and M1 Macrophages, while resting memory CD4 T cells had a strong negative correlation with CD8 T cells, Tregs, and memory CD4 T cells.

Fig. 6

Panorama of immune cell infiltration in the two subgroups and analysis of differences and correlations. A Panorama of infiltration of 22 immune cells in the two groups; B Differences in infiltration of each immune cell type between the two groups (* p < 0.05; ** p < 0.01; *** p < 0.001); C Heat map showing the correlation between the degree of immune cell infiltration

Characterization of immunotherapy-related indicators in the two subtypes

We analyzed the differences in dysfunction and exclusion scores of mutation counts, MSI, and TIDE in relation to the subgroups (Supplemental_Fig_S5) and found that the mutation count was higher in the high-risk subgroup, which predicted the possible existence of better immunotherapy efficacy (Supplemental Figure S5A). The MSI in the high-risk subgroup had a higher score, which also predicted better immunotherapy efficacy (SupplementalFig_S5B). In the TIDE analysis, the high-risk subgroup had a higher dysfunction score, which predicted a higher level of T-cell dysfunction (SupplementalFig. S5C), while the exclusion score was not significantly different between the two groups (Supplemental Figure S5D).

Construction of risk model gene-related regulatory networks

A risk model gene-related PPI network was constructed based on the STRING database (Fig. 7A). A risk gene-related mRNA-miRNA interaction network was constructed based on the miRTarBase database (Fig. 7B). A TF network of risk genes was constructed based on the ChEA database (Fig. 7C). Finally, risk gene-drug interaction information was predicted using the Comparative Toxicogenomics Database to explore drugs targeting risk genes (Fig. 7D).

Fig. 7

Construction of risk gene-related PPI, mRNA-miRNA, gene-TF, gene-drug networks. A A risk gene-based PPI network, where each node represents a different protein, and darker color represents more interactions; B Construction of a risk gene-based mRNA-miRNA network, where orange represents mRNA, darker color represents more interactions, and blue color represents miRNA; C Construction of a risk gene-based TF-gene interaction networks, where red represents genes and blue represents TF; D Construction of a gene-drug networks based on risk genes, where orange represents genes, darker colors represent more interactions with drugs, and blue represents drugs. TF, Transcription factor

Discussion

The search for molecular markers capable of predicting prognosis and immune response in colon cancer is essential for tailoring precision therapy and maximizing patient benefits. Current research has primarily focused on the identification of markers using whole-tissue sequencing. However, the heterogeneity of tumors and the role of mesenchymal tissues in tumors have largely been neglected. Therefore, in this study, we integrated data from the spatial transcriptome, single-cell transcriptome, and bulk transcriptome of colon cancer to delineate colon cancer characteristics. We performed differential analyses across various regions to identify markers to refine prognosis and predict the immune response in colon cancer.

Tumor heterogeneity is a hallmark of CRC [29] in which different tumors or stromal cells exhibit distinct molecular characteristics. This variability may account for the lack of specificity or sensitivity of current CRC markers such as carcinoembryonic antigen (CEA). For example, CRC contains various subsets of fibroblasts, resulting in varied effects on CRC development [30]. In the present study, we identified distinct cell clusters in CRC and observed their distribution using a spatial transcriptome analysis. This observation suggests that combining different dimensions of analysis helps accurately stratify patients and guide treatment decisions.

A recent study concluded that individuals who consume iron-rich red meat daily are at a higher risk of developing colon cancer [31]. Furthermore, studies have shown that ferroptosis plays a critical role in the regulation of immune responses in CRC [32,33,34]. A triple combination of the ferroptosis inducer withaferin A, CXCR2 inhibitor SB225002, and α-programmed death protein 1 (α-PD-1) has demonstrated efficacy in reducing CRC metastatic growth in the liver [35]. In agreement, a recent study reported that increased expression of ferroptosis-related genes GPX4, FTH1, and FTL, coupled with lower ACSL4 expression, correlated with worse 5-year overall survival (OS) among male CRC patients with KRAS mutations [36]. In addition, ferroptosis may hold the potential for predicting both CRC prognosis and immunotherapeutic responses [37]. In light of these reports, ferroptosis-related genes appear to hold promise as therapeutic targets for predicting survival outcomes and immunotherapy responses in CRC.

In the present study, we found that low ferroptosis scores were associated with better PFS and chemotherapy response, and that cetuximab promoted ferroptosis in KRAS-mutant CRC cells. This apparent contradiction suggests that we need a deeper understanding of the mechanism of ferroptosis in different therapeutic settings. The role of ferroptosis may depend on the specific biological environment, tumor subtypes and treatment strategies. In our analysis, the low ferroptosis score represents a specific tumor microenvironment characteristic that may be associated with the sensitivity of some CRC patients to conventional chemotherapy. The anti-tumor mechanism of cetuximab not only involves the inhibition of EGFR signaling pathway, but also plays a role by regulating ROS level, inducing ferroptosis or changing the immune microenvironment, especially in KRAS mutant CRC cells. Therefore, although ferroptosis may be associated with poor prognosis in some cases, it may become a mechanism to enhance the therapeutic effect in specific therapeutic settings, such as the use of cetuximab. In this study, we established a prognostic and immune response prediction model for colorectal cancer based on ferroptosis-related genes by integrating spatial transcriptome, single-cell transcriptome, and large-scale transcriptome data, and verified its robustness in multiple cohorts. We found that ferroptosis-related genes were not only effective in differentiating patients’ prognostic risk, but also closely associated with immune cell infiltration and immunotherapy-related molecular features. In particular, the high expression of SNX5 and USP7 at the tissue level provides new clues for subsequent studies on ferroptosis mechanisms and targeted therapies. In recent years, studies on SNX5 and cancer have gradually increased, showing its potential role in tumor progression [38, 39]. Study has shown that SNX5 was able to affect the proliferation and migration of clear cell renal cell carcinoma cells by inhibiting the internalization of CD44 and its associated epithelial-mesenchymal transition (EMT) process [38]. However, the role of SNX5 in CRC has not been studied yet. USP7 plays an important role in oxidative stress and cell death mechanisms, including ferroptosis, by regulating p53 and its downstream signaling pathways [40]. Although this study confirmed the high expression of SNX5 and USP7 in colorectal cancer tissues by IHC, and found that they are related to the prognosis of patients and the immune microenvironment, the specific regulatory mechanism in the process of ferroptosis is still unclear. At present, there are few reports on whether SNX5 and USP7 affect the progression of colorectal cancer by regulating ferroptosis, so further experimental studies are needed to verify their specific effects. In the future, we plan to explore the regulatory role of SNX5 and USP7 in the process of ferroptosis through in vitro cell experiments, gene knockout/knockdown experiments and ferroptosis inhibitor intervention, so as to further reveal its potential mechanism in CRC. In conclusion, although the specific mechanisms of SNX5 or USP7 in colorectal cancer still need to be further explored, there is evidence that it may be involved in the development of this disease through multiple pathways and may serve as a promising biomarker and therapeutic target. This study emphasizes the unique advantage of multi-omics integration for revealing tumor heterogeneity and microenvironmental features, which lays the foundation for precise typing and individualized treatment of CRC. In the future, the mechanism of key ferroptosis genes will be further explored through functional experiments to promote clinical translational applications.

In addition to providing valuable insights into the role of ferroptosis-related genes in CRC development, this study had a few limitations. First, the availability of datasets for the spatial transcriptomics of CRC is limited and the sample sizes are small. Future large-scale studies of CRC spatial transcriptomes are required. Second, ferroptosis-related genes were screened based on multi-omics data analysis, and the expression levels of some genes were verified by IHC. However, due to the limitation of time and experimental conditions, we have not further explored the specific mechanism of SNX5 and USP7 in the ferroptosis process. In future studies, we plan to further verify the specific regulatory role of SNX5 and USP7 in the ferroptosis process through functional experiments (such as gene knockout, lipid peroxidation level detection, cell viability experiments, etc.) to enhance the scientific and clinical value of the study. Third, although we elucidated the characteristics of immune cell distribution, we did not analyze the interactions between immune cells and tumor cells; Although our study demonstrated good predictive performance on an independent validation set with high reliability and stability of results, we recognize that further work such as increasing the sample size and combining clinical samples and ex vivo and in vivo models have a more complete complement to our study. We plan to consider these additional efforts in future studies to enhance the credibility and comprehensiveness of the study.

Conclusions

Our results identified distinct tumor and mesenchymal regions and demonstrated patterns of immune cell infiltration. Furthermore, we illustrated the link between ferroptosis-related genes and CRC prognosis and immune response through comprehensive bioinformatics analysis at multiple levels.