A risk prognostic model for patients with esophageal squamous cell carcinoma basing on cuproptosis and ferroptosis

The data of RNA sequencing (RNA-seq) and its corresponding clinicopathologic information of ESCC were downloaded from public databases, namely the Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO, https://cancergenome.nih.gov/↗, https://www.ncbi.nlm.nih.gov/geo/↗, GSE53624↗). The FRGs were obtained from the Molecular Signatures Database (http://www.gsea-msigdb.org/gsea/msigdb/↗), while the CRGs were obtained from previous publications by Tsvetkov et al. (2022). As the data were available online with usage allowance, additional ethical approval was unnecessary.

To obtain the ferroptosis and cuproptosis scores of all samples, we conducted a GSVA analysis using the “GSVA” package in the TCGA (Tomczak et al. 2015) and GTEx database (Nangraj et al. 2020). We performed WGCNA and identified 4 optimal cuproptosis and ferroptosis-related genes (OCFRGs) based on the scores obtained above. To determine the optimum soft thresholding power, we performed the function “pickSoftThreshold” in R and found that it was 4. To analyze modules with different power, we calculated the topological overlap matrix (TOM) and diss TOM (1-TOM) to generate the gene dendrogram (Lin et al. 2020). Finally, we conducted a Pearson correlation coefficient analysis for the co-expression modules of cuproptosis-score and ferroptosis-score.

The detailed work process of our study is presented in Fig. 1. After screening out the CFRGs, we performed five machine learning methods, including Decision tree, extreme gradient boosting (XGBoost), gradient boosting decision tree (GBDT), random survival forest, and least absolute shrinkage and selection operator (LASSO), to estimate their prognostic value. Then, we performed univariate Cox regression and constructed a Lasso Cox model using the top-ranked 30 selected genes. The median risk score was considered the optimal cutoff value to divide all patients into two groups (high risk and low risk).

We evaluated the predictive performance of CFRGs using survival status distribution, principal component analysis (PCA), Kaplan–Meier (KM) analysis, and 1-, 3-, and 5-year receiver operating characteristic (ROC) curves, using the “survival”, "survminer," and "timeROC" R packages, respectively (Kong et al. 2020; Wang et al. 2020; Zhu et al. 2020). Finally, we used the GEO cohorts’ datasets as the external validation cohorts.

We conducted cox regression analysis via “ComplexHeatmap” to investigate the correlation between the CFRGs prognostic model and several clinicopathological features, including sex, gender, and TNM stage (Moreno-Ayala et al. 2021). To present the results intuitively, we used “ggalluvial” to estimate the relationship among survival status, risk score, and significant prognostic factors (Ma et al. 2022).

We calculated several prognosis-related tumor-infiltrating immune cells, including CD4 + T cells, CD8 + T cells, B cells, neutrophils, macrophages, and mast cells, in ESCC patients using algorithms such as “TIMER,” “CIBERSORT”, “CIBERSORT-ASB”, “QUANTISEQ”, “MCPCOUNTER”, “XCELL”, and “EPIC” (Lv et al. 2022). To further investigate the difference in immune infiltration level between two risk groups, we conducted single sample gene set enrichment analysis (ssGSEA) using the “ssGSEA” R package (Lv et al. 2022). We then used the “cibersort” algorithm to evaluate the Pearson correlation coefficient of the expression level of four OCFRGs with the immune cells, and the most significant immune cell and its R-value were calculated.

We conducted GO enrichment analyses between the two risk groups to further investigate the potential biological functions and pathways of our prognostic signature of CFRGs. Additionally, we collected several cuproptosis and ferroptosis-related pathways from previous publications, including Hippo, ErbB, Wnt, MAPK, NF-kB, PI3K/AKT, TGF-β, Notch, Ras, AMPK, JAK-STAT, TNF, PD-1/PD-L1, mTOR, and HIF-1 (Gao et al. 2014; Mizushima et al. 2002; Tian et al. 2006; Du et al. 2009; Naganuma et al. 2012; Qiu et al. 2014; Zhang et al. 2020).

Using GSVA enrichment analysis, we selected a hub gene set (https://www.gsea-msigdb.org/gsea/downloads.jsp↗) to conduct correlation analysis with the risk score, aiming to further investigate the correlation between the risk score and enrichment score. Additionally, we used the “GSVA” R package to estimate the active signaling pathway between the two risk groups. Finally, we calculated the differential expression levels of four CFRGs between ESCC tissue and its paired paracancerous samples.

To evaluate the differences in mutation status between high and low risk subgroups, we used the MAFTOOL software to obtain their respective top mutational genes, types, and frequencies of mutations (Mayakonda et al. 2018). Based on the expression data, we obtained immune scores, stromal scores, and estimate scores via the “ESTIMATE” algorithm to assess the relationship between risk score and tumor purity (He et al. 2020). Additionally, we obtained the simple nucleotide variation data of ESCC patients from the TCGA database.

For the in vitro experiments, we used three ESCC cell lines, including KYSE-150, KYSE-510, and EC9706, as well as the human normal esophageal epithelial cell line Het-1A. These cells were maintained in RPMI-1640 (Gibco, China) or DMEM (Gibco, China) supplemented with 10% fetal bovine serum (FBS, Excell Bio, Shanghai, China) in a humidified incubator at 37 °C with 5% CO₂. KYSE-150, KYSE-510, and Het-1A cells were kindly donated by Dr. Xiangyan Liu, and EC9706 cells were purchased from Jennio Biotech, Guangzhou.

We collected four frozen esophageal squamous cell carcinoma and its adjacent normal tissues from Shandong Provincial Hospital of Shandong University. Prior to surgery, all patients underwent gastroscopy, biopsy, and upper abdomen contrast-enhanced computed tomography (CT) to confirm the diagnosis of ESCC, and none of the patients received any anti-cancer treatment before surgery. Tissue specimens were preserved in RNA preservation solution and stored at − 80 °C until used. The study was conducted in accordance with the Declaration of Helsinki, and the protocol was approved by the Shandong Provincial Hospital (Shandong University) of China.

Total RNA from the cell lines, four ESCC tissue, and their paired paracancerous samples were isolated from tissues using the Accurate Biology kit (China). For cDNA, reverse transcription was performed using the TransScript uni all-in-one first-strand cDNA synthesis supermix for qPCR (TransGen Biotech, China), followed by quantitative reverse transcription-polymerase chain reaction (qRT-PCR) on cDNA using 2× Universal Blue SYBR Green qPCR Master Mix (Servicebio, China). The expression value of the target gene was normalized to that of the internal control gene GAPDH (CT GAPDH). All primers used in this research are listed in supplementary file 1.

In addition, we obtained four tissue microarrays (TMA) from Shanghai Outdo Biotech Company (Shanghai, China), with each TMA containing 30 esophageal squamous cell tumor tissues, paracancerous tissues, and distal normal tissues. We completed immunohistochemical staining with the help of Servicebio (Wuhan, China), using primary antibodies for immunohistochemistry including anti-MIDN (BB08169494), anti-COMTD1 (BB12065083), anti-RAP2B (BB12068902) from Bioss (China) and anti-C15orf65 (R98588↗) sourced from Novus Biologicals (USA). The TMA images were acquired using the TissueFaxs software (version 7.1.119) from Tissue-Gnostics^® (Vienna, Austria).

The gene expression score for each sample was calculated by multiplying the percentage of stained cells and the staining intensity, expressed in arbitrary units (A.U.). Staining percentage was scored on a scale of 0–4, depending on the percentage of stained cells: 0 (0–5%), 1 (5–30%), 2 (30–50%), 3 (50–70%), and 4 (> 70%). Staining intensity was scored on a scale of 0 to 5, where 0 represents no staining, 1 represents moderate staining, and 5 represents high-intensity staining. Part of our scoring criteria is shown in Fig. 2A–H, and Fig. 2A, B, D, and E were obtained from the Human Protein Atlas.

All experiments were repeated at least three times, and all the analyses were performed using R software (version 4.1.0). To filter out CFRGs, Pearson correlation coefficient analysis was utilized, and the Kaplan–Meier method was performed for univariate survival analysis (Parashar et al. 2019). The Wilcoxon test was used to evaluate expression levels, select differentially expressed genes (DEGs), estimate the proportion of tumor cells in the tumor tissue in high- and low-risk subgroups, and determine the difference in the proportion of somatic mutations. P values were exact, two-tailed, and values less than 0.05 were considered statistically significant.

The soft threshold power was calculated using the pickSoftThreshold function in the WGCNA R package to select a soft threshold of 10 (Fig. 3A). Modules were determined using a dynamic cutting method, and each candidate module underwent cluster analysis using the “mergeCloseModules” function. Finally, a total of seven modules (represented by different colors) were identified using the Dynamic Tree Cut method. Next, we performed Pearson correlation analysis to investigate the correlation among co-expression modules, tumor size, cuproptosis score, and ferroptosis score. The “brown” module had the strongest correlation [R (Tumor) = 0.4, R (Cuproptosis) = 0.5, R (Ferroptosis) = 0.81, P < 0.001)] (Fig. 3D).

To develop a predictive signature based on the cuproptosis score and ferroptosis score of ESCC patients, we performed univariate Cox analysis on genes in module “MEbrown” to identify survival-related CFRGs. The weights of every survival related CFRG were evaluated using five machine learning algorithms, including random forest, decision tree, LASSO, XGBoost, and GBDT. We then selected the top 30 CFRGs to screen genes and build a prognostic model using LASSO COX regression analysis. Four genes were chosen for the establishment of the CFRG prognostic model based on the optimal setting for the tuning parameter λ, among which C15orf65, COMTD1, and RAP2B were considered to be causative or risk factors, and MIDN was indicated as a protective factor. The risk score was calculated using the following formula

Independent validation datasets were obtained from the GEO database (http://www.ncbi.nlm.nih.gov/geo/↗), and patients were stratified into a high-risk group and a low-risk group based on the median of their risk scores in the training cohort and analyzed (Fig. 4A). The risk score distribution, gene expression, and OS status of patients in the training and validation cohorts are shown in Fig. 4B. PCA analysis indicated that different risk groups were distributed in two directions (Fig. 4D). After performing KM analysis, we found that the OS was significantly different (P < 0.001) between high- and low-risk subgroups in both TCGA and GEO cohorts (Fig. 4C). ROC curve analysis was performed to assess the sensitivity and specificity of the risk model, and the area under the ROC curve (AUC) was 0.799 for 1 year, 0.722 for 3 years, and 0.715 for 5 years in the training set (Fig. 4E) and 0.724 for 1 year, 0.608 for 3 years, and 0.675 for 5 years in the validation set (Fig. 4E).

To examine the different clinical characteristics between two groups, sex, gender, and TNM stage were included as covariates, and the correlations did not reach significance except for tumor staging (Fig. 5A). Using the Sankey diagram, we could unambiguously assess that patients in the high-risk group experience a larger share of death, and the most associated clinical factor with death is tumor staging (Fig. 5B). To facilitate a better comprehension of the subject matter, we have prepared a comprehensive heatmap (Fig. 5C) that delineates the intricate interplay between tumor purity, immune infiltration, and clinical features.

To further investigate the relationship between infiltrating immune cells and gene expression levels, we used the “CIBERSORT” algorithm and found that MIDN expression was significantly positively correlated with B cells naive, and C15orf65 expression was negatively correlated with plasma cells, while both genes were positively correlated with antigen-presenting cells. Regarding the two other genes, COMTD1 expression was linearly related to macrophages M1, and RAP2B showed a positive correlation with activated mast cells. The detailed correlation between gene expression and immune cell infiltration based on CIBERSORT is listed in Fig. 6.

Gene ontology (GO) was used to enrich the correlated biological process, cellular component, and molecular function. The results showed that the enrichment of DEGs mainly presented in T cell activation, calcium-mediated signaling, and antigen processing and presentation pathways. The most related cellular components were MHC class II protein complex, MHC protein complex, and the external side of the plasma membrane. Furthermore, for molecular function, immune receptor activity, MHC class II protein complex binding, and MHC protein complex binding showed evidence of strong positive correlation (Fig. 7C). GSVA enrichment focused primarily on several important pathways, and the relationship between these pathways and the risk score was displayed in Fig. 7A. We found that the Wnt, MAPK, PI3K/AKT, TGF-β, AMPK, JAK-STAT, PD-1/PD–L1, mTOR, TNF, HIF-1, and ErbB pathways presented a positive correlation with the risk score (Fig. 7A), and Wnt, PI3K/AKT, TGF-β, AMPK, JAK-STAT, PD-1/PD-L1, mTOR, TNF, HIF-1, and ErbB showed strong correlations. Regarding the different activated pathways between high- and low-risk groups, we performed GSVA to enrich the activated pathways between the two subgroups and found that cytokine receptor interaction, intestinal immune network for IGA production, and complement and coagulation cascades were significantly enriched (Fig. 7B). The Wilcoxon rank-sum test was used to evaluate expression levels, and we found that MIDN, COMTD1, and RAP2B were highly expressed in tumor tissues and lowly expressed in normal tissues. However, the difference in COMTD1 expression did not reach significance. For C15orf65, we observed the opposite tendency (Fig. 7D–G). We observed a similar tendency in the qRT-PCR and IHC experiments. Files for qRT-PCR and IHC are available in supplementary files 2 and 3.

With respect to the tumor immune microenvironment, we found that in the high-risk subgroup, activated dendritic cells (aDCs), B cells, dendritic cells (DCs), immature dendritic cells (iDCs), macrophages, mast cells, neutrophils, NK cells, plasmacytoid dendritic cells (pDCs), follicular helper T cells, Th1 cells, Th2 cells, tumor-infiltrating lymphocytes (TILs), and regulatory cells (Tregs) had significantly higher infiltrating abundance (Fig. 8A). In terms of immune-related functions between the high- and low-risk groups, several important functions were significantly different, such as T cell co-inhibition, T cell co-stimulation, and checkpoint functions (Fig. 8B).

To further explore the function of our established model, we analyzed the mutation status of the high- and low-risk groups in the TCGA cohort. We found that TP53, TNN, KMT2D, and NFE2L2 were the most frequently mutated genes, and regardless of the mutated genes, missense mutation was the most common mutation type. The detailed relationship among mutation type, number of samples, and TMB in the high- and low-risk groups is displayed in Fig. 8D, E.

Tumor purity is composed of the immune score and estimate score. To analyze the relationship between the risk score and tumor purity, we used the “ESTIMATE” algorithm to obtain the immune score, stromal score, and estimate score and assessed the proportion of tumor cells in the tumor tissue in the high- and low-risk subgroups using the Wilcoxon rank-sum test. Immune score (R = 0.39), estimate score (R = 0.33), and stromal score (R = 0.22) were all positively correlated with the risk score and showed significant differences.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Data are available in a public, open access repository. Data are available upon reasonable request. All data relevant to the study are included in the article or uploaded as supplementary information.