A total of 471 ferroptosis-related genes (FRGs) were downloaded from the FerrDb database (http://www.zhounan.org/ferrdb). Original counts of RNA-seq transcriptome data and clinical data derived from TCGA-COAD tissues (n=473) were obtained from The Cancer Genome Atlas (TCGA) database (https://tcga-data.nci.nih.gov/tcga/). The GSE17536 (n=177), GSE17537 (n=55), and GSE39084 (n=70) from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) were downloaded. The 8 rectal cancer samples in GSE39804 were deleted.

The downloaded TCGA and GEO data were merged while maintaining common genes. The “sva” package of the R software was used to remove batch effects between different datasets. The expression matrix of 390 FRGs was obtained by the intersection of the expression matrix of merged cohort and the 471 FRGs. We employed consensus clustering using k-means algorithms to identify clusters. The number and stability of clusters were established by implementing the “ConsensuClusterPlus” package (10). Moreover, the categorization was repeated 1000 times to ensure its accuracy and stability (10). Lastly, the “survival” package in R software was used to evaluate the overall survival (OS) among different clusters through the Kaplan-Meier curve. Principal component analysis (PCA) could determine whether the clusters can be intuitively separated.

Differentially expressed genes (DEGs) were identified between different clusters using the “limma” package with an adjusted P value of 0.05 and an absolute value of |log₂FC|=0.585. Using the Gene Ontology (GO) database, information of DEGs in cell components (CC), biological processes (BP), and molecular functions (MF) was identified (11). Using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database, we were able to obtain information on possible signaling pathways, genes, and diseases associated with DEGs (12). DEGs were enriched and visualized by GO and KEGG analyzes using the packages “clusterProfiler”, “org.Hs.eg.db”, “enrichplot”, and “ggplot2”. Gene Set Variation Analysis (GSVA) enrichment analysis was used to analyze the differences in signaling pathways between different clusters (13). The”GSVA”package was used to identify the 20 most enriched signaling pathways that are differentially expressed in different clusters.

Univariate Cox analysis was used to identify 16 DEGs. Patients in the merged cohort were randomly divided into a training cohort and a testing cohort. Then, the Least Absolute Shrinkage and Selection Operator (LASSO) regression analysis was used to eliminate genes that might overfit the model and minimize variables (14). Multivariate Cox analysis was used to identify the prognostic genes. The risk model was constructed by three prognostic genes (SLC2A2 CDKN2A and FABP4). Risk score = Coef_gene1 * Exp_gene1+Coef_gene2 * Exp_gene2+…+Coef_genen * Exp_genen. The Coef_gene was the risk coefficient of each gene calculated by the LASSO Cox model, and the Exp_gene was each gene expression level. Patients in the training cohort were divided into high- and low-risk groups based on the median risk score. Survival analysis was performed using the “survival” and “survminer” packages to compare the OS of patients in the high-and low-risk groups. Time-dependent receiver operating characteristic (ROC) curve analysis was performed with the “timeROC” package to evaluate the predictive value of the prognostic signature. Next, PCA and t-distributed stochastic neighbor embedding (t-SNE) analysis were performed using the “Rtsne” and “ggplot2” packages to determine whether the high- and low-risk groups can be intuitively separated. The testing cohort and the merged cohort were used to validate the prognostic value of the three genes.

To further validate the prognostic value of the three FRGs in CC. The Tumor-immune System Interactions Database (TISIDB) database (http://cis.hku.hk/TISIDB/) was used to analyze the differences in expression of prognostic genes at different stages of CC. The Kaplan-Meier curves of the high and low expression levels of the three prognostic genes in CC were obtained from the Gene Expression Profiling Interactive Analysis (GEPIA) database (http://gepia.cancer-pku.cn/), respectively.

To explore the differences in immune cells between high- and low-risk groups, we analyzed the correlation between tumor infiltrating immune cells (TIICs) and prognostic genes, risk score. Cell-type Identification by Estimating Relative Subsets of RNA Transcript (CIBERSORT) was an algorithm that revealed the composition ratio of 22 TIICs in TIME and was used to analyze the correlation between the risk score and the TIICs by the “limma” and “tidyverse” packages (15). The 22 TIICs of gene expression files could be obtained from the official CIBERSORT website (http://cibersort.stanford.edu/).

In order to assess the value of risk scores for the immunotherapy effect. The tumor immune dysfunction and exclusion (TIDE) analysis (http://tide.dfci.harvard.edu) was an algorithm that calculates tumor response to immune checkpoint inhibitors. The TIDE score integrated both mechanisms of T cell dysfunction and T cell rejection in tumor immune escape (16). The higher the TIDE score, the worse the immunotherapy effect and the prognosis of patients. Immunophenoscore (IPS) is used to predict the therapeutic effect of immune checkpoint inhibitors in cancer patients and is calculated based on the expression of four major types of genes that determine immunogenicity (17). The higher the IPS value, the better the patient had received immunotherapy with CTLA-4 monoclonal antibody or PD-1 monoclonal antibody (17). Obtain IPS of CTLA-4 monoclonal antibody and PD-1 monoclonal antibody in CC patients from the cancer immunome atlas (TCIA) (https://tcia.at/home), and compare immunotherapy effect between high- and low-risk groups.

The CC cell lines (HCT116, DLD-1, and LOVO) and the normal colon epithelial cell line (NCM460) were obtained from the American Type Culture Collection. These cell lines were cultured in RPMI-1640, McCoy’s 5A or F12-K medium (Gibco) containing 10% FBS (Hyclone) and 1% penicillin–streptomycin (Gibco) and incubated at 37°C with 5% CO2.

In this study, a retrospective cohort study was conducted to collect postoperative CC specimens and clinical data from CC patients who underwent surgery in the Department of General Surgery of Tianjin Medical University General Hospital from January 2021 to April 2021. Inclusion criteria for clinical patients: a) primary adenocarcinoma of the colon confirmed by postoperative pathology; b) negative surgical margins, complete resection of the primary lesion; c) complete clinical and pathological data, and follow-up can be completed; Exclusion criteria: a) patients with non-primary colon cancer or patients with other primary malignancies; b) patients who underwent emergency surgery; c) patients who were lost to follow-up postoperatively; d) patients whose postoperative specimens were not collected. All patients were followed up for 2 years, and follow-up ended when the patient died. Patient survival was determined from the time of surgery until the end of follow-up. According to statistics, from January 2021 to April 2021, a total of 69 patients with CC were treated surgically, except for 18 patients who did not have postoperative specimens, 3 patients who underwent emergency surgery, 3 patients who were lost to follow-up, and 2 patients with CC complicated with rectal cancer. In this study, a total of 43 patients were collected. Of the 43 CC patients, 21 patients had tumors located in the right colon, and the remaining 22 patients had tumors located in the left colon. All of them received radical tumor resection surgery, and 27 patients completed chemotherapy with capecitabine plus oxaliplatin (XELOX) after surgery. There were 3 patients with cancer-related death and 3 patients with liver metastases after CC surgery.

Total RNA was extracted using TRIzol reagent (Vazyme) following the manufacturer’s protocols. The concentration of RNA was measured using a NanoDrop-2000 spectrophotometer (Thermo Fisher Scientific, Inc.) and reverse transcribed using a FastQuant RT Super mix kit (Tiangen Biotech, Co, Ltd.). Subsequently, qPCR was performed using SYBR-Green qPCR Master Mix (Bimake Biotechnology) under the following thermal cycle conditions: initial denaturation at 95°C for 15 min, followed by 40 cycles at 95°C for 10 s, 60°C for 20 s, and 72°C for 20 s. The data were quantified using 2−ΔΔCT method. The GAPDH gene was used to normalize target gene expression. The PCR primer used were as follows:

GAPDH forward primer: 5′‐TGGCACCGTCAAGGCTGAGAA‐3′;

GAPDH reverse primer: 5′‐TGGTGAAGACGCCAGTGGACTC‐3′;

SLC2A3 forward primer: 5′‐GCATCGTTGTTGGAATTCTGGT‐3′;

SLC2A3 reverse primer: 5′‐TGTAGGATAGCAGGAAGGATGG‐3′;

CDKN2A forward primer: 5′‐GGGTTTTCGTGGTTCACATCC‐3′;

CDKN2A reverse primer: 5′‐CTAGACGCTGGCTCCTCAGTA‐3′;

FABP4 forward primer: 5′‐ACTGGGCCAGGAATTTGACG‐3′;

FABP4 reverse primer: 5′‐CTCGTGGAAGTGACGCCTT‐3′.

In this study, all statistical analyses were performed using R software (version 4.1.2) and perl-5.34.0. The t tests, nonparametric tests, and chi-square tests were used to test for differences between variables, where appropriate. The Kaplan–Meier method was used to draw survival curves for the high- and low-risk groups in the training, testing and merged cohorts, respectively, and the log-rank test was used to determine the statistical significance of differences (18). A p value <0.05 indicated a statistically significant difference.

When integrating the TCGA and GEO datasets, a total of 24 repeated samples from the 13 CC patients in the TCGA cohort were deleted. The consensus clustering analysis was used to identify clusters from 743 CC patients of the merged cohort to further evaluate the prognostic genes in CC. When increasing the cluster variable (k) from 2 to 10, the intragroup correlations were found to be the highest when k=2 ( Figures 1A–C ). There were 398 patients in cluster-A and 345 patients in cluster-B. The two clusters can be clearly distinguished by PCA analysis ( Figure 1D ). The Kaplan-Meier curve was drawn between the two clusters, and the difference was statistically significant (p<0.001) ( Figure 1E ). The clinicopathological characteristics, such as age, sex, grade, pathological stage, T stage, N stage, and M stage, were compared in the heatmap between the 2 clusters, and no significant difference was found in the clinicopathological characteristics ( Figure S1 ). In addition, GSVA enrichment analysis revealed that cluster-B was enriched in the peroxisome pathway compared with cluster-A ( Figure 1F ). These results indicated that the two clusters identified based on FRGs accurately reflected the prognostic difference.

To explore the potential biological behavior of ferroptosis in the two clusters, DEGs between cluster-A and cluster-B were analyzed. There were 35 DEGs found between the two clusters, including 4 protective genes and 31 risk genes. GO and KEGG enrichment analyses were implemented to gain insights into signal pathways of clusters. The GO enrichment analysis indicated that risk genes were predominantly enriched in multicellular organismal homeostasis (biological process), membrane raft (cellular component) and heme binding (molecular function) ( Figure 2A ). The GO enrichment analysis indicated that protective genes were predominantly enriched in response to oxidative stress (biological process), NADPH oxidase complex (cellular component) and superoxide−generating NADPH oxidase activity (molecular function) ( Figure 2B ). In addition, the KEGG pathway analysis showed that the risk genes were enriched in the ferroptosis, IL-17 signaling pathway, NOD-like receptor signaling pathway and HIF-1 signaling pathway ( Figure 2C ). The KEGG pathway analysis showed that the protective genes were enriched in the thyroid hormone synthesis ( Figure 2D ). The DEGs were mainly enriched in ferroptosis-related pathways, and ferroptosis might play a dual effect on tumor cells.

We constructed a risk signature to better diagnose the prognosis of CC. With filtering the 23 CC samples with missing survival status, missing survival time and 0 survival time in the TCGA database. A total of 720 CC samples were used for the construction and validation of the risk model. Using the univariate Cox analysis, 16 DEGs were identified as prognostic signatures ( Figure 3A ). Patients in the merged cohort were randomly divided into training cohort (n=360) and testing cohort (n=360). Based on the expression levels of the 16 DEGs and clinical data, three genes significantly associated with the prognosis of CC were identified by applying the LASSO regression analysis and the multivariate Cox analysis, including SLC2A3, CDKN2A, and FABP4 ( Figures 3B, C ). The risk score was calculated according to the following formula: Risk score = SLC2A3 Exp*0.136407939746576 + CDKN2A Exp *0.175273043748396 + FABP4 Exp *0.165558961729006. According to the risk model, the heatmap, risk score distribution plot, and survival status scatter plot were used to evaluate the model ( Figures 3D–F ; Table S1 ). The CC patients in the training cohort were divided into high- and low-risk groups based on the median risk score. The number of patients with death status increased with the risk score. The risk score of cluster-A was significantly higher than the risk score of cluster-B (p<2.22e-16, Figure 3G ). Then, the Kaplan-Meier plot showed that low-risk patients had a significantly better OS than high-risk patients (p<0.001) ( Figure 3H ). In this study, the area under curve (AUC) of training cohort predicting the sensitivity and specificity of 5-year OS was 0.718 ( Figure 3I ). Finally, after dimensionality reduction by the PCA and the t-SNE analyses, it was found that the CC samples could clearly visualize the distribution between the high- and low-risk groups ( Figures 3J, K ). We identified three FRGs and constructed a prognostic signature that could well predict the prognosis of CC patients.

In order to validate the prognostic value of the risk model, the same methods were performed in the merged cohort and testing cohort ( Figures 4A–F ; Table S1 ). Then, the Kaplan–Meier plot revealed that the high-risk group had a significantly worse OS compared with the low-risk group in the both groups (p_merged<0.001, p_testing<0.001) ( Figures 4G, H ). The AUCs for 5-year OS were 0.706 and 0.673, respectively, indicating that the risk model had a strong prognostic value for CC patients in the merged cohort and testing cohort ( Figures 4I, J ). Lastly, the high- and low-risk groups can be clearly distinguished by PCA and t-SNE analyses ( Figures 4K–N ). Both the merged cohort and the testing cohort had well validated the value of the prognostic signature.

To identify the potential value of SLC2A3, CDKN2A, and FABP4 on the prognosis of CC, we generated gene expression differential diagrams of CC stages I-IV and Kaplan-Meier survival curves from the TISIDB database and the GEPIA database, respectively. The expression levels of SLC2A3, CDKN2A and FABP4 genes increased with the growth of CC stage (rho_{SLC2A3 =} 0.108, p_{SLC2A3 =} 0.0229, rho_CDKN2A=0.152, p_CDKN2A=0.00126, rho_{FABP4 =} 0.122, p_{FABP4 =} 0.0102) ( Figures 5A–C ).High expression of CDKN2A and FABP4 were associated with poor prognosis in CC, but there was no statistical difference in survival curves of SLC2A3 (p_{SLC2A3 =} 0.077, p_CDKN2A=0.011, p_{FABP4 =} 0.021) ( Figures 5D–F ). The analysis results of the TISIDB and the GEPIA databases further proved that the three FRGs were all risk genes in CC.

The status of the tumor TIME determines the immunotherapy effect. The results suggest that there is a close correlation between prognostic genes and TIICs was investigated. The correlation between prognostic genes and TIICs was essentially consistent with the risk score. The risk score was mainly positively correlated with neutrophils and macrophages infiltration and negatively correlated with infiltration of CD4⁺T cells, CD8⁺T cells, regulatory T cells, B cells, plasma cells and dendritic cells (P <0.001) ( Figure 6 ). In this study, there are significant correlations between SLC2A3, CDKN2A, FABP4, risk score, and TIICs.

To further validate the value of the risk model for immunotherapy, we predicted response to immunotherapy in high- and low-risk groups. TIDE analysis calculated response to immunotherapy from CC in the TCGA cohort. The TIDE score, dysfunction score, and exclusion score of the high-risk group were significantly higher than those of the low-risk group (p_TIDE<0.005, p_Dysfunction<0.005, p_Exclusion<0.001) ( Figures 7A–C ). The IPS score was used to analyze potential response to immunotherapy. The results showed that the low-risk group had higher IPS score after treatment with CTLA-4 monoclonal antibody or PD-1 monoclonal antibody (p_{mAb-CTLA-4 =} 3e-08, p_{mAb-PD-1 =} 4.1e-10) ( Figures 7D, E ). Therefore, we predicted that the immunotherapy effect might be better in low-risk groups than in high-risk groups.

In order to further validate the prognostic value of the risk model, we performed RT-qPCR on 43 CC samples and obtained the relative expression level of prognostic genes in each sample ( Table S2 ). The risk score of 43 CC samples was calculated. There were no statistically significant differences in clinicopathological characteristics between high- and low-risk groups ( Table S3 ). The expression levels of prognostic genes were higher in the high-risk group than in the low-risk group (p_{SLC2A3 =} 0.0002, p_CDKN2A<0.0001, p_{FABP4 =} 0.0001) ( Figures 8A–C ). Combined with the 2-years follow-up data of 43 CC patients, the OS and DFS were calculated between the high- and low-risk groups. There was no statistical difference in OS, but there was a statistical difference in DFS (p_OS=0.0831, p_DFS=0.0108) ( Figures 8D, E ). Validation of the risk model with clinical samples further demonstrated the prognostic value of the risk model for CC patients, and the high-risk group patients might have the worse OS and DFS.

In this study, the expression differences between normal tissues and cancer tissues were also compared. The diagram for difference analysis of expression of three prognostic genes in normal and cancer tissues was downloaded from the GEPIA database. The result showed that there was no significant difference in the expression of SLC2A3 between normal and cancer tissues; the expression of CDKN2A in cancer tissues was higher than that in normal tissues (p<0.01); the expression of FABP4 in cancer tissues was lower than that in normal tissues (p<0.01) ( Figures 9A–C ). The expression of three prognostic genes was verified by RT-qPCR in CC cell lines. The expression results of the three prognostic genes between CC cell lines (HCT116, DLD-1, and LOVO) and normal cell line (NCM460) were consistent with the expression differences between tissues in the GEPIA database (p_{CDKN2A-HCT116}<0.0001, p_{CDKN2A- LOVO}<0.0001, p_{CDKN2A-DLD-1 =} 0.0013, p_FABP4-HCT116<0.0001, p_FABP4-LOVO<0.0001, p_FABP4-DLD-1<0.0001) ( Figures 9D–F ). Therefore, CDKN2A and FABP4 not only had the prognostic value for CC patients, but also were expected to become the therapeutic targets of CC in the future.