Identification and validation of a prognostic signature of cuproptosis-related genes for esophageal squamous cell carcinoma

Data download and arrangement

We used the TCGABiolinks R package [14] to download the ESCC dataset (TCGA-ESCC) from The Cancer Genome Atlas database (TCGA, https://portal.gdc.cancer.gov/) and analyzed it as the test set. We removed samples lacking clinical information data and retained only primary tumor (-01A) tissue types, resulting in a count of sequencing data comprising 80 ESCC samples with clinical information. We normalized sequencing data from the TCGA cohort from raw counts to fragments per kilobase million (FPKM). Corresponding clinical information of patients was obtained from the UCSC Xena database (http://genome.ucsc.edu) [15], and mutation data were taken from MAF files downloaded from the TCGA database. Furthermore, we used a unified standardized TCGA-ESCC dataset, which integrated the TCGA and GTEx databases and the sequencing data in the form of normalized read counts as Transcripts per Million (TPM). The TCGA-ESCC dataset consisted of 745 samples (653 normal samples from GTEx and 92 tumor tissues from TCGA), and was used for subsequent differential expression analysis of CRGs between normal and ESCC groups.

Furthermore, we downloaded two microarray gene-expression profiles, GSE2034 [16] and GSE38129 [17], from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) [18] using the GEOquery package [19] for use as the verification cohort. These cohorts were obtained from Homo sapiens and were converted into corresponding genes using the annotation information available from the gene platforms, namely GPL571 (HG-U133A_2) Affymetrix Human Genome U133A 2.0 Array. The GSE20347 cohort included 17 para-carcinoma tissue samples and 17 carcinoma tissue samples, while the GSE38129 cohort included 30 tumor tissue samples and 30 tumor-adjacent tissues. Additionally, we used the gene expression dataset of ESCC patients obtained from our own sequencing for correlation analysis.

The individual ESCC dataset comprised 12 samples, including expression profile data from six esophageal squamous carcinoma tissue samples (grouping: ESCC) and six adjacent normal esophageal tissues (grouping: Normal). These samples were obtained from surgical patients at Fujian Cancer Hospital between July 1st and July 17th, 2020, and none of the patients had received treatment before surgery. The research scheme (including specimen collection) was reviewed and approved by the Biomedical Ethics Committee of Fujian Cancer Hospital (batch number: K2021-027-01). The clinical data for sequencing of patients with esophageal squamous cell carcinoma in our institution was shown in Supplementary Table 1. The GSE20347, GSE38129, and ESCC datasets were used as validation sets for subsequent validation.

We then searched for twelve CRGs from published literature [4] by using “Cuproptosis-related genes” as the search keyword on PubMed. The 12 CRGs were ATPase copper transporting beta (ATP7B), cyclin dependent kinase inhibitor 2A (CDKN2A), dihydrolipoamide S-acetyltransferase (DLAT), dihydrolipoamide dehydrogenase (DLD), FDX1, glutaminase (GLS), lipoic acid synthetase (LIAS), LIPT1, metal regulatory transcription factor 1 (MTF1), PDHA1, pyruvate dehydrogenase E1 (PDHB), and solute carrier family 31 member 1 (SLC31A1).

Differential expression analysis and consensus clustering analysis for cuproptosis-related genes

First, we compared the expression profiles (HTSeq-TPM) of cuproptosis-related genes (CRGs) between para-carcinoma tissues and carcinoma tissues. Next, we preprocessed somatic mutation data obtained from the TCGA-ESCC cohort using the “VarScan” software and visualized it using the “maftools” R package [20].

To verify the rationality of clustering, we employed the consensus clustering algorithm, which involves multiple iterations over subsamples of the dataset to provide an indication of cluster stability and parameter decisions. We used subsampling to induce sampling variability [21] and classified ESCC into various clusters to explore the function of CRGs. We applied the “ConsensusClusterPlus” R package [22] with 50 iterations, 80% resampling rate Pearson correlation, clusterAlg = “km”, and distance = “Euclidean” (available at http://www.bioconductor.org/). We investigated gene expression arrays in the ESCC groups to compare the group expression difference of CRGs among different disease subtypes from the TCGA-ESCC cohort using the Wilcoxon rank sum test. CRGs with a P-value < 0.05 were considered as threshold values for identification.

Gene set variation analysis

Gene Set Variation Analysis (GSVA) [23] is a non-parametric unsupervised analysis method primarily used to assess the results of gene set enrichment in microarrays and transcriptomes. It evaluates whether different pathways are enriched between samples by converting the gene expression matrix between samples into the expression matrix of gene sets between samples. To further analyze the differences in biological behavior between ESCC disease subtypes, we conducted GSVA analysis. We used the gene set “h.all.v7.4.symbols.gmt” from the Molecular Signature Database (MSigDB) [24] as the reference set for GSVA. A threshold value of P < 0.05 was considered significant.

Single-sample gene-set enrichment analysis and cuproptosis score

The single-sample gene-set enrichment analysis (ssGSEA) algorithm quantifies the relative enrichment of each gene in a given dataset sample. We used the R package GSVA to calculate the cuproptosis score (CPs) using the ssGSEA algorithm based on the expression matrix of 12 CRGs for each sample in the TCGA-ESCC, GSE20347, GSE38129, and ESCC datasets. We then analyzed the expression differences in CPs between the normal subtype and ESCC subtype in the GSE20347 and GSE38129 cohorts, as well as the expression differences in grouping of different subtypes of ESCC disease in the TCGA-ESCC dataset using the Wilcoxon rank sum test. A P-value < 0.05 was set as the cutoff for significance.

Construction of a prognostic risk model based on CRGs

To construct a prognostic signature and screen prognostic CRGs in the TCGA-ESCC dataset, we utilized the glmnet package [25] to perform Least absolute shrinkage and selection operator (LASSO) Cox regression [26]. The LASSO regression is based on linear regression and adds a penalty term (absolute value of lambda × slope) to reduce overfitting and improve the generalization ability of the model. We performed the analysis using CRGs in the TCGA-ESCC dataset with a random seed number of “2021” and tenfold cross-validation. The results were visualized, and we illustrated the grouping of each sample in the LASSO regression prognostic model according to risk score and survival outcome, as well as the molecular expression of prognostic CRGs in each subgroup by risk factor plots.

The ESCA patients were stratified into low-risk and high-risk subgroups based on the median risk score in the TCGA cohort. Kaplan-Meier curves were generated using the “survival” package to compare the overall survival (OS) differences between the low- and high-risk groups. Furthermore, receiver operating characteristic (ROC) curves [27] and corresponding area under the curve (AUC) values were calculated using the “timeROC” package [28] to evaluate the performance of the prognostic risk model for predicting patient outcomes.

Drug sensitivity analysis

We collected ESCC cell lines and their drug response data from The Genomics of Drug Sensitivity in Cancer (GDSC, https://www.cancerrxgene.org/) database [29]. To compare the difference in immune drug sensitivity between high- and low-risk patients, we used the “oncoPredict” package [30] to perform drug sensitivity analysis in the TCGA-ESCC cohort. Additionally, the “DESeq2” R package [31] was utilized to screen for differentially expressed genes (DEGs) between the high- and low-risk groups in the TCGA-ESCC cohort, with screening criteria of an adjusted P-value < 0.05 and | log2 Fold Change (FC)| ≥1.0. Upregulated genes were screened using the criteria of logFC > 1.0 and an adjusted P-value < 0.05. Downregulated genes were screened using the criteria of logFC < −1.0 and an adjusted P-value < 0.05. We visualized the results of the expression analysis using the “ggplot2” R package.

Genomic and functional analysis

Gene Ontology (GO) functional analysis [32] is a widely used approach for large-scale functional enrichment studies, encompassing biological process (BP), molecular function (MF), and cellular component (CC) categories. In our study, we employed the “ClusterProfiler” package [33] to conduct GO enrichment analysis of the differentially expressed genes (DEGs) between low- and high-risk groups in TCGA-ESCC cohort. We considered a statistically significant result when the false discovery rate q-value was less than 0.10 and the P-value was less than 0.05. The Benjamini-Hochberg method (BH) was utilized to correct the P-values.

Gene Set Enrichment Analysis (GSEA) [34] is a computational method that determines whether a defined set of genes exhibits statistically significant concordant differences between two biological states. In this study, GSEA was performed using the “clusterProfiler” R package to elucidate significant functional and pathway differences between the high and low expression groups of DEGs in the TCGA-ESCC cohort. Each analysis procedure was repeated 1000 times. The c2.cp.v7.2.symbols set was downloaded from the MSigDB database and used as the reference gene set. A functional or pathway term with an adjusted P-value < 0.05 and a false discovery rate (FDR) < 0.25 was considered to have statistically significant enrichment.

Immune infiltration analysis

CIBERSORTx [35] is a suite of machine learning tools for assessing immune cell types and their abundance from mixed cell populations. In this study, we uploaded the gene expression data to the CIBERSORTx online analytical platform (https://cibersortx.stanford.edu/) to assign values for immune cell infiltration using the LM22 gene signature. We used boxplots to visualize differences in immune cell infiltration between low- and high-risk groups in the TCGA-ESCC cohort. The Spearman algorithm and “ggplot2” package were used to visualize correlations between low- and high-risk groups. The “pheatmap” package in R software was used to visualize correlations of infiltrating immune cells and CRGs among the low- and high-risk groups in the TCGA-ESCC cohort using a heatmap.

Construction and validation of the prediction nomogram

To define the clinical value of the CRGs prognostic model constructed by LASSO regression for ESCC, we performed univariate Cox regression analysis on the expression of prognostic CRGs in the TCGA-ESCC dataset. We included all factors from the univariate Cox regression analysis in the multivariate Cox regression analysis to construct a multivariate Cox regression model. Based on the results of the multivariate Cox regression analysis, we established a nomogram using the rms R package. A nomogram [36] is a graph that represents the functional relationship among multiple independent variables through a cluster of disjoint lines in a rectangular coordinate system. It is based on the multivariate Cox regression analysis by setting a certain scale to score and characterize the individual variables within the multivariate Cox regression model, and finally calculating the total score to predict the probability of the event occurrence situation. Furthermore, we used decision curve analysis (DCA) to evaluate the accuracy and discrimination of the nomogram. DCA [37] is a simple method to evaluate clinical predictive models. We also used the R package “ggDCA” [38] to draw the DCA curve to assess the effect of the COX regression model.

Statistical analysis

All statistical analyses and plots were performed using R software (Version 4.1.2). The difference between two groups was compared using Student’s t-test or Wilcoxon rank-sum test. The Kruskal-Wallis test was used to compare differences among three or more groups. For categorical variables, Fisher’s exact test and chi-square test were employed to compare differences between two groups. Spearman’s correlation coefficient was used to assess correlation. All P-values were two-sided, and P < 0.05 was considered statistically significant.

Differential expression of cuproptosis-related genes in ESCC dataset

Firstly, we collected RNA-seq data of 12 CRGs (ATP7B, CDKN2A, DLAT, DLD, FDX1, GLS, LIAS, LIPT1, MTF1, PDHA1, PDHB, SLC31A1) from the TCGA-ESCC dataset, GSE20347 dataset, GSE38129 dataset, and ESCC dataset. We then analyzed the expression differences of these 12 CRGs between ESCC tumor tissues and normal tissues in the four ESCC datasets using the Mann-Whitney U test (Wilcoxon rank-sum test). The results of the expression difference correlation analysis were presented through subgroup comparison plots (Figure 1). Our findings revealed that in the TCGA-ESCC dataset, nine out of the 12 CRGs (ATP7B, DLAT, FDX1, GLS, LIAS, MTF1, PDHA1, PDHB, SLC31A1) showed highly statistically significant differences (P < 0.001) between the ESCC and normal groups, except for CDKN2A, DLD, and LIPT1 genes (P > 0.05) (Figure 1A). In the GSE20347 dataset, the expression of CRGs DLAT, GLS, LIAS, MTF1, and PDHB was highly statistically significant (P < 0.001) between the ESCC and normal groups. Additionally, the expression of FDX1 gene was also highly statistically significant (P < 0.01) between different subgroups (Figure 1B).

We performed differential expression analysis in the GSE38129 dataset and found that seven out of the 12 CRGs showed statistically significant differences in expression between the ESCC and normal groups, with GLS, LIAS, MTF1, and PDHB genes being highly statistically significant (P < 0.001) (Figure 1C). In addition, in the ESCC dataset, the expression of GLS, LIAS, MTF1, PDHB, and SLC31A1 CRGs was statistically significant (P < 0.05) between the ESCC and normal groups (Figure 1D). Taken together, these results suggested that the 12 CRGs may play an important role in ESCC based on the differential expression analysis of these genes in the four ESCC datasets.

Genetic variants of cuproptosis-related genes in ESCC

To investigate somatic mutations in 12 CRGs (ATP7B, CDKN2A, DLAT, DLD, FDX1, GLS, LIAS, LIPT1, MTF1, PDHA1, PDHB, SLC31A1) in ESCC from the TCGA-ESCC cohort, we summarized the results of genetic mutation analysis. Four regulators (CDKN2A, ATP7B, DLD, LIAS) were found to have somatic mutations. The mutations mainly included missense and splice mutations, with missense mutations being the most common. Additionally, all genetic mutations of the 4 CRGs in ESCC were shown to be single nucleotide polymorphisms (SNPs), with C>T being the top-class SNV (Figure 2A).

Figure 2B showed that 7 ESCC samples had somatic mutations, and 5 (71.43%) of the 7 ESCC samples had mutations in the 4 regulators (CDKN2A, ATP7B, DLD, LIAS). CDKN2A mutation samples accounted for 29% of all ESCC genetic mutations, involving missense and splice mutations. The remaining 3 regulator (ATP7B, DLD, LIAS) mutation samples accounted for 14% of all ESCC genetic mutations and only consisted of missense mutations.

Furthermore, we used the RCircos package to map the chromosomal location of the 12 CRGs (Figure 2C). The 12 CRGs were mainly located on chromosomes 1, 2, 4, 9, 11, and X. CRGs were most commonly distributed on chromosomes 2, 9, and 11. Chromosomes 2, 9, and 1 each had 2 CRGs, indicating that CRGs located on the same chromosome might be closely related at the genomic level.

Development and validation of ESCC subtypes

Initially, we investigated the differential expression of CRGs in ESCC patients from the TCGA-ESCC cohort by using the “ConsensusClusterPlus” R package for consensus clustering analysis based on 12 CRGs (ATP7B, CDKN2A, DLAT, DLD, FDX1, GLS, LIAS, LIPT1, MTF1, PDHA1, PDHB, SLC31A1). This led to the separation of tumor samples into 2 clusters (Figure 3A). Cluster 1 consisted of 59 tumor samples, while cluster 2 contained 21 tumor samples. Next, we performed principal component analysis (PCA) on the expression matrix from these 2 ESCC subtypes. The results showed significant differences in the cuproptosis transcription profiles between cluster 1 and cluster 2 (Figure 3B). We then used a heatmap to identify the gene expression pattern of each subtype (Figure 3C), which revealed a few differences in expression level of the 12 CRGs between the two ESCC subtype from the TCGA-ESCC cohort.

Furthermore, we compared the prognosis signature of these 2 ESCC subtypes using Kaplan-Meier analysis, and the survival analysis showed that cluster 1 had the worst survival situation (Figure 3D). We then used the Wilcoxon rank sum test to compare the level of expression of the 12 CRGs (ATP7B, CDKN2A, DLAT, DLD, FDX1, GLS, LIAS, LIPT1, MTF1, PDHA1, PDHB, SLC31A1) between cluster 1 and 2 from the TCGA-ESCC cohort, and the results are shown in Figure 3E. As a result in Figure 3E, the expression of two regulators (CDKN2A, SLC31A1) showed significant differences between cluster 1 and 2 (P < 0.05).

Gene set variation analysis and construction of cuproptosis score prognosis and diagnostic signature

To explore the biological differences between cluster 1 and 2, we conducted GSVA enrichment analysis. Compared to cluster 1, cluster 2 exhibited significantly higher GSVA enrichment scores in two pathways, namely MYC targets V1 and unfolded protein response. The difference in these pathways between the two clusters was further analyzed using Mann-Whitney U test and visualized in Figure 4A. Additionally, we used the “pheatmap” R package to create a heatmap to present the specific difference analysis results of these two pathways between cluster 1 and 2, as shown in Figure 4B. The heatmap results confirmed that there were significant differences in these two pathways between cluster 1 and 2 in the TCGA-ESCC cohort.

Spearman correlation analysis was used to investigate the expression levels of 12 CRGs (ATP7B, CDKN2A, DLAT, DLD, FDX1, GLS, LIAS, LIPT1, MTF1, PDHA1, PDHB, and SLC31A1) from the TCGA-ESCC cohort, and the results are presented in Figure 4C. As shown in Figure 4C, only the expression level of GLS was negatively correlated with that of other CRGs, while the expression levels of other CRGs were positively correlated. Additionally, we displayed the correlation scatter plots for the two pairs of genes with the highest and lowest correlation (Figure 4D, 4E). Among the 12 CRGs in the TCGA-ESCC cohort, DLAT and FDX were significantly positively correlated (r = 0.536, P < 0.001, Figure 4D), while GLA and PDHB were significantly negatively correlated (r = −0.388, P < 0.001, Figure 4E).

Furthermore, to investigate the extent of copper-related cell death in ESCC patients, we utilized the single-sample gene-set enrichment analysis (ssGSEA) algorithm to calculate the Cuproptosis score (CPs) of ESCC patients from the TCGA-ESCC dataset. Subsequently, we used the Wilcoxson rank sum test to examine the difference of CPs between the two subtypes (cluster 1 and 2) in the TCGA-ESCC cohort. We observed a statistically significant difference of CPs between cluster 1 and 2 (P < 0.001) (Figure 4F). Then, based on the CRs, we divided the TCGA-ESCC dataset samples into two subtypes (high score group and low score group). Next, we employed KM analysis to investigate the prognostic significance between the two groups using the clinical information of samples from the TCGA-ESCC dataset. As demonstrated in Figure 4G, the CPs might not be able to accurately predict the survival outcome of ESCC patients (P > 0.05).

In addition, based on the transcriptome expression matrix of 12 CRGs from three verification datasets (GSE20347, GSE38129, and ESCC), we calculated the correlation expression values among the 12 CRGs using Pearson correlation analysis, and presented the results in heatmaps (Figure 5A, 5D, 5G). According to these heatmaps, we found that the GLS values were negatively correlated with the majority of the other CRGs values in these three verification datasets, while the other CRGs were positively correlated with one another.

Moreover, based on the expression level of the 12 CRGs from these three verification datasets, we utilized the ssGSEA algorithm to evaluate the CPs of each sample. Then, we used the Wilcoxon rank sum test to quantify the difference of CPs between the normal group and the ESCC group in the three verification datasets. The correlation analysis revealed a statistically significant difference of CPs between the normal group and the ESCC group in GSE20347 and GSE38129 (Figure 5B, 5E). Conversely, we found that the difference of CPs was not statistically significant between different groups in the ESCC dataset (Figure 5H).

To explore the clinical utility of CPs evaluation in these three datasets, we utilized a ROC curve to demonstrate its ability to discriminate ESCC diagnosis. As indicated by the ROC curves, the AUC values in the GSE20347 cohort, GSE38129 cohort, and ESCC cohort were 0.862 (Figure 5C), 0.771 (Figure 5F), and 0.750 (Figure 5I), respectively. It can be inferred that the CPs has moderate accuracy in diagnosing ESCC in these three verification datasets.

Construction and prognostic risk model based on 11 cuproptosis-related genes

To evaluate the predictive value of 12 CRGs (ATP7B, CDKN2A, DLAT, DLD, FDX1, GLS, LIAS, LIPT1, MTF1, PDHA1, PDHB, SLC31A1) for clinical outcomes in the TCGA-ESCC dataset, we used the transcriptome expression matrix of 79 patients with ESCC who had a survival time greater than zero to construct a prognosis signature through the LASSO algorithm (Figure 6A, 6B). Supplementary Table 2 shows the signature coefficients for each gene under different penalty coefficients in the LASSO regression. The LASSO regression analysis, which is based on linear regression, removes overfitting and improves generalization ability by increasing the penalty term (lambda times the absolute value of the slope). Based on the LASSO regression analysis, we used 11 CRGs, including ATP7B, CDKN2A, DLAT, DLD, FDX1, LIAS, LIPT1, MTF1, PDHA1, PDHB, and SLC31A1, according to the minimum-minimum criterion and optimal-minimum criterion. We also calculated the penalty coefficients of the CRGs using LASSO analysis and established a risk index by multiplying the gene expression by the corresponding coefficients. We further visualized the sample groupings in the constructed prognostic model for CRGs using risk factor plots (Figure 6C). Kaplan-Meier survival analysis revealed that overall survival was significantly lower in the high-risk group of ESCA patients than in the low-risk group (P < 0.05, Figure 6D). The AUC values demonstrated that the risk scores significantly predicted the prognosis of patients with ESCA. As shown in Figure 6E, the AUC values for 1-, 2-, and 3-year were 0.847, 0.958, and 1.000, respectively.

Evaluation of the sensitivity of antitumor drugs using the signature

The TCGA-ESCC cohort was stratified into low- and high-risk groups using the median risk score. To investigate the differences in drug sensitivity between these groups, we compared the IC50 values of multiple antitumor drugs in both low- and high-risk groups (Figure 7A–7G). Our drug sensitivity analysis revealed significant differences between the low- and high-risk groups in seven drugs, including BMS.536924 (Figure 7A), BMS.754807 (Figure 7B), CGP.60474 (Figure 7C), NVP.TAE684 (Figure 7D), PF.02341066 (Figure 7E), PLX4720 (Figure 7F), and Sunitinib (Figure 7G). Notably, Sunitinib (Figure 7G, P < 0.05) was more effective in patients in the low-risk group, whereas the other six drugs were more effective in patients in the high-risk group.

Furthermore, we compared the cuproptosis score (CPs) between the low- and high-risk groups using the Mann-Whitney U test and found a significant difference between the two groups (P < 0.01, Figure 7H). We then performed a correlation analysis between the risk score and cuproptosis score in patients from the TCGA-ESCC dataset, and the results were visualized by a correlation scatterplot (Figure 7I). Our findings showed that the risk score was negatively correlated with the cuproptosis score in patients from the TCGA-ESCC dataset (r = −0.387, P < 0.001, Figure 7I).

To identify the differentially expressed genes (DEGs) associated with risk score in these patients from the TCGA-ESCC dataset, we compared the transcription matrix expression levels in the low-risk and high-risk groups using the DESeq2 package. After analyzing the transcriptomic changes, a total of 18,057 DEGs were identified in the volcano plot (Figure 7J) (Supplementary Table 3). We further screened 71 critical genes (|log (FC)| >1 and FDR <0.05) using R software. Among these, we found 36 upregulated (downregulated in the low-risk group, log (FC) >0) and 35 downregulated (upregulated in the low-risk group, log(FC) <0) transcription factors in the high-risk group.

Gene ontology (GO) and gene set enrichment analysis (GSEA)

After being converted into gene IDs, the 71 differentially expressed genes (DEGs) were analyzed using Gene Ontology (GO) (Supplementary Table 4). The GO annotations of the DEGs consisted of three parts: CC (cellular component), BP (biological process), and MF (molecular function), which were used to analyze the functional enrichment of the DEGs. The results were visualized in a bubble plot (Figure 8A). The main biological processes involved adult locomotory behavior and neuron migration, while the most abundant cellular component terms were the integral component of postsynaptic membrane, intrinsic component of postsynaptic membrane, integral component of synaptic membrane, and integral component of postsynaptic specialization membrane. The most abundant molecular function terms were extracellular matrix structural constituent. Additionally, we also visualized the results of GO functional enrichment analysis with a circos plot (Figure 8B).

Subsequently, based on the functional enrichment analysis, we calculated the z-score of each gene of the DEGs by providing logFC values obtained from the difference analysis of TCGA-ESCC dataset. The results of GO functional enrichment analysis, which were standardized by logFC values, were presented in a bubble plot (Figure 8C), mainly enriched in cellular component (CC) pathway. The relationship between these pathways and genes was shown in Figure 8D.

We further performed Gene Set Enrichment Analysis (GSEA) on the expression of the DEGs in TCGA-ESCC dataset to uncover signaling pathways that are differently active in ESCC. The GSEA showed that the DEGs of TCGA-ESCC dataset were significantly enriched in four pathways (Figure 8E–8I, Supplementary Table 5, P < 0.05 and q-value < 0.25), including the proteasome degradation pathway (Figure 8F), Biocarta classic pathway (Figure 8G), complement activation pathway (Figure 8H), and integrin 3 pathway (Figure 8I).

Characterization of immune infiltration between the low- and high-risk groups

To explore differences in immune cell infiltration between the low- and high-risk groups, we used CIBERSORT to evaluate the abundance of 22 immune cell types in the TCGA-ESCC cohort. The distribution of overall immune cell abundances among the samples from the TCGA-ESCC cohort is shown in Figure 9A, and the results indicate that Macrophages M0 and Macrophages M1, as well as T cells CD4 memory activated, are the main infiltrating immune cells.

Furthermore, we examined the relationships among the infiltration abundances of the 22 distinct types of immune cell (Monocytes, B cells memory, Plasma cells, B cells naive, T cells follicular helper, Eosinophils, T cells CD4 memory activated, T cells CD8, NK cells activated, T cells regulatory Tregs, Macrophages M1, T cells gamma delta, Mast cells resting, Dendritic cells resting, Macrophages M2, T cells CD4 memory resting, Macrophages M0, T cells CD4 naive, Mast cells activated, NK cells resting, Dendritic cells activated, and Neutrophils) based on the samples with abundances greater than 0 in the low- and high-risk groups (Figure 9B, 9C). Negative correlations were found between most of the 22 types of immune cells in the low-risk group (Figure 9B), while negative correlations were observed between most of the 21 types of immune cells, except B cells memory, in the high-risk group (Figure 9C).

Additionally, we investigated the correlation between the expression of 11 CRGs (ATP7B, CDKN2A, DLAT, DLD, FDX1, LIAS, LIPT, MTF1, PDHA1, PDHB, and SLC31A1) and the infiltration abundance of 22 types of immune cells in the low-risk (Figure 9D) and high-risk groups (Figure 9E), respectively. The screening criterion for correlation analysis was a P value < 0.05. Our results showed that CRG expression levels were statistically positively correlated with the infiltration abundance of Mast cells resting in the low-risk group, while the infiltration abundances of Mast cells activated and T cells CD4 naive were positively correlated with PDHB expression levels in the low-risk group (Figure 9D). Conversely, most immune cell infiltration abundances were positively correlated with the expression levels of 12 CRGs in the high-risk group (Figure 9E).

Validation of the prognostic value of the signature

To further validate the prognostic risk model of CRGs, we conducted a statistical analysis of clinical data from ESCC patients obtained from the TCGA-ESCC dataset (Supplementary Table 6). Firstly, we performed univariate and multivariate Cox regression analyses to demonstrate the correlation between patients’ prognosis and the expression levels of 11 CRGs (ATP7B, CDKN2A, DLAT, DLD, FDX1, LIAS, LIPT1, MTF1, PDHA1, PDHB, SLC31A1). In this study, we first conducted univariate Cox regression analysis based on the expression levels of 11 CRGs, followed by multivariate Cox regression analysis including all variables in the univariate Cox analysis (Table 1). Our results revealed that 8 CRGs, namely CDKN2A, DLAT, FDX1, LIAS, LIPT1, PDHA1, PDHB, and SLC31A1, were significantly correlated with ESCC survival (P < 0.05). The results of the univariate and multivariate Cox regression analyses were then presented in a forest plot (Figure 10A).

Subsequently, based on the statistically significant prognostic factors identified in the multivariate analysis, a nomogram prognostic evaluation model was constructed. A nomogram characterizes multiple variables in the multivariate regression model, and the total score is calculated to predict the probability of events. As shown in Figure 10B, unlike other CRGs, the expression levels of three CRGs, including ATP7B, DLAT, and DLD, exhibited limited usefulness in the Cox regression model.

Afterwards, a Decision Curve Analysis (DCA) was performed to evaluate the clinical value of the LASSO-Cox regression prognosis model for 1-year (Figure 10C), 2-year (Figure 10D), and 3-year (Figure 10E) time points. In the DCA curve, the x-axis represents the threshold probability, while the y-axis displays the net benefit. The results were determined by observing the x-value range of the model line that is higher than all positive line and all negative line. The larger the range of X-value, the better the clinical judgment utility. The DCA analysis revealed that the risk score model constructed by 11 CRGs (ATP7B, CDKN2A, DLAT, DLD, FDX1, LIAS, LIPT1, MTF1, PDHA1, PDHB, SLC31A1) had the best clinical judgment utility for 3-year followed by 2-year and 1-year.

We verified the diagnostic performance of our risk model using the TCGA-ESCC dataset. Firstly, we calculated the risk score for each sample in the low- and high-risk groups based on the LASSO-Cox regression prognosis model. We then used Mann-Whitney U tests to compare the risk scores between the two groups, and found a significant difference (P < 0.001, Figure 11A). Subsequently, we generated an ROC curve for the CRG prediction model to evaluate the ability of the risk score system to predict outcomes for low- and high-risk patients (Figure 11B). We observed that the risk score levels exhibited high diagnostic accuracy for grouping patients into low- and high-risk categories (AUC = 0.998, Figure 11B).

Based on the results of multivariate Cox regression analysis, we identified and included 11 CRGs in the construction of the LASSO-Cox regression prognosis model. The formula was as follows:

RiskScores = ATP7B × 1.114 + CDKN2A × −1.180 + DLAT × −2.618 + DLD × −1.091 + FDX1 × 1.302 + LIAS × −2.163 + LTPT1 × −2.366 + MTF1 × 0.326 + PDHA2 × 2.349 + PDHB × 3.816 + SLC31A1 × 1.499.

To validate the diagnostic accuracy of the LASSO-Cox regression prognosis model, we utilized Mann-Whitney U tests to calculate the difference in risk scores between the ESCC and normal groups in the TCGA-ESCC cohort (Figure 11C, 11D). Our analysis demonstrated a highly significant difference in risk score between the two groups (P < 0.001, Figure 11C), indicating the model’s effectiveness in distinguishing between the two groups. Additionally, we generated an ROC curve to evaluate the diagnostic capacity of the risk score system, with the resulting analysis shown in Figure 11D. The area under the curve (AUC) for the TCGA-ESCC prognostic model was 0.842.

Furthermore, based on the expression levels of 11 CRGs across three cohorts (GSE20347, GSE38129, and ESCC), we calculated the risk score for each sample within the ESCC and normal groups using the formula. We then used Mann-Whitney U tests to compare the risk score between ESCC and normal groups in each of the three datasets. Our results showed no statistically significant difference in risk score between the two groups in any of the datasets, including GSE2034 (Supplementary Figure 1A), GSE38129 (Supplementary Figure 1C), and ESCC (Supplementary Figure 1E) (P > 0.05). Finally, ROC curve analysis was conducted, and the resulting AUC values for the prognosis model score in GSE2034 (Supplementary Figure 1B), GSE38129 (Supplementary Figure 1D), and ESCC (Supplementary Figure 1F) were 0.623, 0.577, and 0.694, respectively.