RNA-seq (HTSeq-Counts/log2(FPKM+1)), Illumina Human Methylation 450K BeadChip (HM450, $\beta$-value), miRNA-seq (log2(TPM)), and clinical data of COAD and normal samples were downloaded from the TCGA database using the UCSC Xena platform (https://xena.ucsc.edu/). The transcriptome dataset included 446 COAD samples and 39 normal samples, the methylome dataset included 297 COAD samples and 18 normal samples, and the miRNAome dataset included 425 COAD samples and seven normal samples. The expression levels of mRNAs, lncRNAs, and miRNAs were transformed into FPKM and TPM values. Gene expression microarray data for COAD were downloaded from the Gene Expression Omnibus (GEO) database (GSE39582) [17]. The matrix of each dataset was normalized using the ‘normalize.quantiles’ function in the R “preprocessCore” (v1.64.0) package (https://bioconductor.org/).

In total, 288 ferroptosis-related genes, including drivers, suppressors, and markers were obtained from the FerrDb database [18]. Among them, genes with relationships between expression and COAD survival were identified using univariate Cox survival analysis. The ssGSEA function in the R package “GSVA” (v1.50.0) (https://bioconductor.org/) was used to calculate the enrichment scores of gene sets positively or negatively associated with COAD survival [19]. The FS, representing the ferroptosis level in each COAD sample, was defined as the difference in enrichment scores between positive and negative components.

Cancer cell line drug sensitivity data (IC50) were downloaded from the Genomics of Drug Sensitivity in Cancer database (GDSC, https://www.cancerrxgene.org/) [20, 21]. The drug sensitivity score of each COAD patient was estimated using the R package “oncoPredict” (v0.2) (https://cran.r-project.org/) based on gene expression data [22]. Drugs with a mean predicted sensitivity value $\le$0.5 in patients were retained. Differences in response sensitivity for each drug between the two groups were tested using Student’s t-test. Correlations between ferroptosis scores (FSs) and estimated drug sensitivity scores were determined using Pearson correlation analysis. A p-value 0.05 was considered statistically significant.

Survival differences between the two groups were calculated using the Kaplan-Meier method and log-rank test in the R package “survival” (v3.5-7) (https://cran.r-project.org/). Samples were classified into two subgroups based on the median FS or risk score values.

DE-mRNAs and DE-lncRNAs were identified using the R package “DESeq2” (v1.42.0) (https://bioconductor.org/) based on the thresholds: $|$log2Fold-change$|$ $>$1 and adjusted p 0.01 between the tumor and normal groups and $|$log2Fold-change$|$ $>$0.5 and adjusted p 0.01 between the high-FS and low-FS COAD subgroups. DE-miRNAs were identified using the R package “limma” (v3.58.1) (https://bioconductor.org/) based on the thresholds: adjusted p 0.01 between the tumor and normal groups and adjusted p 0.01 between the high-FS and low-FS COAD subgroups. DMCs were also identified using the R package “limma” based on the thresholds: $|$methylation difference$|$ $>$0.2 and adjusted p 0.01 between the tumor and normal groups and $|$methylation difference$|$ $>$0.1 and adjusted p 0.01 between the high-FS and low-FS COAD subgroups.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed for the DE-mRNAs using DAVID (https://david.ncifcrf.gov/tools.jsp) [23]. A Benjamini-Hochberg (BH)-adjusted p-value 0.05 was considered statistically significant.

LASSO Cox regression model analysis was used for variable selection with the “glmnet” R package. Multivariate analysis based on the Cox proportional hazards regression model was used to calculate the coefficients of the selected variables and construct the prognostic signatures. Receiver operating characteristic (ROC) curves and area under the ROC curves (AUC) were used to evaluate the effectiveness of the prognostic model with the “survivalROC” R package (v1.0.3.1) (https://cran.r-project.org/).

The “mixOmics” R package was employed to select key variables within the three -omics data signatures that effectively differentiate between the high-FS and low-FS COAD subgroups according to the provided instructions (http://mixomics.org/) [24]. In brief, the continuous quantitative values of variables, such as the expression levels of lncRNAs and miRNAs and methylation levels of CpGs, within the three signatures served as input data for the “mixOmics” package (v6.26.0) (https://bioconductor.org/). Supervised analysis and feature selection with sparse Partial Least Square-Discriminant Analysis (sPLS-DA) were performed to identify specific variables of different molecular types critical for discriminating the high-FS and low-FS groups. All statistical analyses were performed using the R program (v4.3.1) (https://cran.r-project.org/).

Table 1 displays the clinical information of the 446 COAD patients and 39 normal samples. Survival analysis showed that age and pathological stage of the COAD patients were both associated with COAD survival (p 0.05). Univariate Cox regression analyses identified eight ferroptosis-related genes with a positive association and 12 with a negative association between their expression and COAD survival (p 0.05). The FS index of ferroptosis activity was estimated based on the expression of the 20 ferroptosis-related genes using ssGSEA (see Materials and Methods). Results showed that the FSs ranged from -0.87 to 0.65 (Supplementary Fig. 1). The young group (age 60 years) had higher FSs compared to the old group (age $>$60 years), but no significant differences in FSs were found between the male and female samples (Fig. 1A,B). Notably, FSs were lower in patients at the late pathological stage, indicating an association between FSs and tumor progression (Fig. 1C). The COAD patients were then divided into two subgroups according to median FSs. Kaplan-Meier survival analysis indicated that the high-FS patients showed good survival (p = 0.00011) (Fig. 1D), as validated by multivariate Cox regression analysis after correcting for age, sex, and pathological stage (Fig. 1E).

Fig. 1.

Association between ferroptosis scores (FSs) and colon adenocarcinoma (COAD) overall survival. (A–C) FS differences in patients of different age, sex, and pathological stage. (D) Kaplan-Meier survival curve of high- and low-FS patients. (E) Forest plot of association between FSs and COAD survival in multivariate Cox proportional hazards model. **, p 0.01, ***, p 0.001.

To explore the potential association between FSs and COAD treatment precision, the drug sensitivity score of each patient was calculated using the R package “oncoPredict” based on drug sensitivity data from the GDSC database. In total, 11 agents showed sensitivity differences between the high-FS and low-FS COAD subgroups, with high drug sensitivity in the high-FS patients (Fig. 2A). Among these drugs, several have been used in COAD therapy, such as camptothecin, vinblastine, and staurosporine [25, 26, 27]. Based on further correlation analysis, 10 of the anti-cancer drugs exhibited significant associations between sensitivity and FSs (Fig. 2B).

Fig. 2.

Drug sensitivity between high- and low-FS COAD subgroups. (A) Differences in predicted drug sensitivity between high- and low-FS subgroups (*, p 0.05). (B) Correlations between FSs and predicted sensitivity values of anti-cancer drugs in COAD patients.

Next, changes in protein-coding mRNA expression were analyzed between the two COAD subgroups. In total, 1065 DE-mRNAs exhibited differences between the COAD patients and normal samples, as well as between the high-FS and low-FS COAD subgroups (see Materials and Methods). Functional enrichment analysis of the DE-mRNAs was conducted to explore differences in biological functions and signaling pathways between the two subgroups. Notably, GO analysis identified 376 up-regulated DE-mRNAs enriched in the inflammatory response, chemokine-mediated signaling pathway, and immune response (Fig. 3A), and 689 down-regulated DE-mRNAs enriched in the Wnt signaling pathway, lipid metabolism, and cell adhesion (Fig. 3B). Furthermore, KEGG pathway analysis revealed that the up-regulated DE-mRNAs were mainly enriched in immune-related pathways, such as the TNF signaling pathway, chemokine signaling pathway, and cytokine-cytokine receptor interaction (Fig. 3C), while down-regulated DE-mRNAs were enriched in the cell adhesion molecule, chemical carcinogenesis, and cAMP signaling pathways (Fig. 3D). Therefore, these findings highlight the close association between ferroptosis and immune processes in tumor progression.

Fig. 3.

Functional enrichment analyses of differentially expressed mRNAs (DE-mRNAs) between high- and low-FS COAD subgroups. (A,B) GO enrichment analysis of up-regulated and down-regulated DE-mRNAs in COAD patients with high-FS. (C,D) Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of up-regulated and down-regulated DE-mRNAs in COAD patients with high-FS.

Ferroptosis-related mRNA signatures associated with COAD prognosis were determined utilizing expression data of identified DE-mRNAs from the TCGA database. LASSO Cox regression analysis was performed on the 1065 DE-mRNAs, with 12 DE-mRNAs selected to build a risk model (risk-score${}_{(\text{mRNA})}$). Subsequently, the multivariate Cox model was used to calculate the coefficients, as follows:

$\text{Risk}-\text{score}(\mathrm{mRNA})$
$=(0.5975)\times \mathrm{ATP6V1B1}+(-0.4779)\times \mathrm{FZD3}+(0.8540)\times \mathrm{GABRD}$
$+(-0.2871)\times \mathrm{GPR15}+(0.8738)\times \mathrm{LRRN4}+(-0.8904)\times \mathrm{MS4A15}$
$+(-2.1418)\times \mathrm{NRG1}+(0.2937)\times \mathrm{NSMF}+(0.1770)\times \mathrm{PCOLCE2}$
$+(0.2700)\times \mathrm{PHACTR3}+(0.5333)\times \mathrm{PLAC1}+(0.2952)\times \mathrm{RSPO4}$

The risk score for each COAD patient was then calculated, with median values used to categorize patients into high- and low-risk groups. Kaplan-Meier analysis revealed a significantly poorer prognosis in the high-risk group compared to the low-risk group (Fig. 4A). The ROC curves indicated AUC values for the risk score model at 1, 3, and 5 years of 0.833, 0.799, and 0.835, respectively (Fig. 4B). Additionally, an independent gene expression microarray dataset of COAD from the GEO database was analyzed to validate the prognostic value of the 12-DE-mRNA signature. This analysis confirmed that COAD patients with a low risk score demonstrated better survival, with AUCs of 0.582 at 3 years and 0.592 at 4 years (Fig. 4C,D). In addition, we also observed that the patients with high FS presented good survival in this independent COAD cohort (Supplementary Fig. 2).

Fig. 4.

Predictive efficacy of ferroptosis-related mRNA signatures in COAD patients from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases. (A,B) Survival and receiver operating characteristic (ROC) curves of the 12-DE-mRNA signature in COAD patients from the TCGA. (C,D) Survival and ROC curves of the 12-DE-mRNA signature in COAD patients from GSE39582.

To further explore the ferroptosis-related regulatory features associated with COAD prognostic outcomes, 91 DMCs, 324 DE-lncRNAs, and 39 DE-miRNAs were identified displaying differences between COAD patients and normal samples, as well as between the high-FS and low-FS COAD subgroups (see Methods). A subset of 132 high-FS and 160 low-FS COAD patients with transcriptomic, miRNAomic, and methylomic data were selected for subsequent analyses. LASSO Cox regression analysis was then performed on the 91 DMCs, 324 DE-lncRNAs, and 39 DE-miRNAs to identify the most predictive features for COAD prognosis. In total, 10 DMCs, 12 DE-lncRNAs, and 14 DE-miRNAs were selected as variables for the risk score models at different molecular levels (i.e., risk-score${}_{(\text{DNAm})}$, risk-score${}_{(\text{lncRNA})}$, and risk-score${}_{(\text{miRNA})}$). The coefficients of the selected DMCs, DE-lncRNAs, and DE-miRNAs were recalculated in each multivariate Cox model, leading to the construction of three risk prediction models (Supplementary Fig. 3), detailed as follows:

$\text{Risk}-\text{score}(\text{DNAm})$
$=(-1.6307\times \mathrm{cg}00987461)+(0.995\times \mathrm{cg}02415779)$
$+(-1.3444\times \mathrm{cg}03616722)+(1.1853\times \mathrm{cg}06120945)$
$+(1.4154\times \mathrm{cg}14192291)+(0.8314\times \mathrm{cg}15507690)$
$+(-0.8891\times \mathrm{cg}18693345)+(1.1020\times \mathrm{cg}21308365)$
$+(-1.5025\times \mathrm{cg}21494776)+(0.5229\times \mathrm{cg}25588389)$

$\text{Risk}-\text{score}(\text{lncRNA})$
$=(-0.46348\times \mathrm{AC016831}.7)+(0.51557\times \text{ATP2A1\_AS1})$
$+\left(1.26266\times \mathrm{CCDC144NL}\mathrm{\_}\mathrm{AS1}\right)+\left(-1.69841\times \mathrm{RP11}\mathrm{\_}109\mathrm{J}4.1\right)$
$+\left(0.33265\times \mathrm{RP11}\mathrm{\_}1143\mathrm{G}9.5\right)+\left(0.35439\times \mathrm{RP11}\mathrm{\_}353\mathrm{N}14.2\right)$
$+\left(-0.36330\times \mathrm{RP11}\mathrm{\_}386\mathrm{I}14.4\right)+\left(0.41455\times \mathrm{RP11}\mathrm{\_}60\mathrm{A}8.1\right)$
$+\left(0.08611\times \mathrm{RP11}\mathrm{\_}638\mathrm{I}8.1\right)+\left(0.69306\times \mathrm{RRP4}\mathrm{\_}673\mathrm{M}15.1\right)$
$+\left(0.14529\times \mathrm{RP5}\mathrm{\_}1059\mathrm{L}7.1\right)+\left(0.55365\times \mathrm{RP5}\mathrm{\_}884\mathrm{M}6.1\right)$

$\text{Risk}-\text{score}\text{ (miRNA})$
$=(0.49643\times \text{hsa\_let\_7e})+(-0.24373\times \text{hsa\_mir\_144})$
$+(0.24081\times \text{hsa\_mir\_146a})+(0.20736\times \text{hsa\_mir\_147b})$
$+(0.37224\times \text{hsa\_mir\_16\_1})+(0.26853\times \text{hsa\_mir\_217})$
$+(-0.07030\times \text{hsa\_mir\_223}+(0.13606\times \text{hsa\_mir\_3074})$
$+(0.14476\times \text{hsa\_mir\_33b})+(-0.37786\times \text{hsa\_mir\_3615)}$
$+(-0.31860\times \text{hsa\_mir\_3677})+(-0.18818\times \text{hsa\_mir\_375)}$
$+(0.05530\times \text{hsa\_mir\_577})+(0.20584\times \text{hsa\_mir\_625)}$

The risk score of each COAD patient was then calculated at three molecular levels using the three different datasets. According to the median values of the risk scores inferred by each molecular signature (i.e., DNAm, lncRNA, and miRNA), the COAD patients were categorized into high- and low-risk groups (Fig. 5A,E,I). Kaplan-Meier analyses revealed significantly poorer prognosis for patients in the high-risk group than those in the low-risk group for each -omics data type (Fig. 5B,F,J). Supplementary Fig. 4 presents the vital statuses of COAD patients classified into the high- and low-risk groups, as determined by the signatures from the three -omics-level datasets. The risk score models from the three datasets were then used to predict patient survival status at 1, 3, and 5 years. ROC curve analysis showed that the AUC values for the DNAm risk score model were 0.725, 0.718, and 0.774 at 1, 3, and 5 years, respectively (Fig. 5C); the AUC values for the lncRNA risk score model were 0.797, 0.739, and 0.773 at 1, 3, and 5 years, respectively (Fig. 5G); and the AUC values for the miRNA risk score model were 0.717, 0.721, and 0.707 at 1, 3, and 5 years, respectively (Fig. 5K). In addition, heatmaps presented the methylation and expression level changes of the 10 DMCs, 12 lncRNAs, and 14 miRNAs within corresponding prediction models alongside their respective risk scores (Fig. 5D,H,L). Notably, a significant proportion of molecular features, including 8 DMCs, 11 lncRNAs and 8 miRNAs, demonstrated statistically significant correlations between their methylation and expression levels and the associated risk cores (p 0.05; Supplementary Table 1). These findings suggest that ferroptosis-related multi-omics features are effective in predicting COAD prognostic outcomes.

Fig. 5.

Construction of ferroptosis-related multi-omics signatures (10 differentially methylated CpG sites (DMCs), 12 DE-lncRNAs, and 14 DE-miRNAs) in COAD patients. (A–D) Risk score plot, survival curve, ROC curve, and heatmap plot of 10-DMC signatures. (E–H) Risk score plot, survival curve, ROC curve, and heatmap plot of 12-DE-lncRNA signatures. (I–L) Risk score plot, survival curve, ROC curve, and heatmap plot of 14 DE-miRNA signatures. In this panel, figures (C,G,K) show predictive performance of three molecular signatures for COAD survival.

To explore the potential of ferroptosis-related multi-omics features as biomarkers for COAD classification, the “mixOmics” R package was used to integrate three omics datasets. This integration encompassed 292 samples with corresponding RNA-seq, HM450, and miRNA-seq data, and included the previously identified 10 DMCs, 12 DE-lncRNAs, and 14 DE-miRNAs. Results showed that the high-FS and low-FS COAD subgroups were distinguishable based on the first and second components of the models for the DNAm, lncRNA, and miRNA datasets, with AUC values of 0.77, 0.76, and 0.80, respectively (Fig. 6A). The loading weights of each selected variable on each component and each dataset are presented in Fig. 6B. In addition, the Circos plot showed significant correlations between and within variables at different -omics levels (r = 0.2) (Fig. 6C). These results indicate that the multi-omics signatures can classify COAD patients into two subtypes based on their ferroptosis activity.

Fig. 6.

Discriminant analysis of COAD patients. (A) ROC curves of three molecular signatures (viz., DNAm, lncRNA, and miRNA) discriminating between high- and low-FS patients. Plot suggests that variables (i.e., markers) selected by sparse Partial Least Square-Discriminant Analysis (sPLS-DA) on components 1 and 2 in each of the signatures can discriminate high-FS samples from low-FS samples. (B) Loading weight plots of variables selected by sPLS-DA on components 1 and 2 in each signature. Variables are ranked according to their loading weight (most important at bottom to least important at top). Colors indicate group for which a particular variable has a higher value than the other group. (C) Circos plot of positive or negative correlations between selected variables from three signatures (presented by brown and black links, respectively). Plot represents correlations greater than 0.2 between variables of different molecular types. Outer lines show levels of each variable in the two groups (i.e., high-FS vs. low-FS).