Identification and validation of necroptosis-related prognostic gene signature and tumor immune microenvironment infiltration characterization in esophageal carcinoma
BMC Gastroenterology volume 22, Article number: 344 (2022)
Esophageal carcinoma (ESCA) is a common malignancy with a poor prognosis. Previous research has suggested that necroptosis is involved in anti-tumor immunity and promotes oncogenesis and cancer metastasis, which in turn affects tumor prognosis. However, the role of necroptosis in ESCA is unclear. This study aimed to investigate the relationships between necroptosis-related genes (NRGs) and ESCA.
Methods and results
The clinical data and gene expression profiles of ESCA patients were extracted from The Cancer Genome Atlas (TCGA), and 159 NRGs were screened from the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. We then identified 52 differentially expressed NRGs associated with ESCA and used them for further analysis. Gene ontology (GO) and KEGG functional enrichment analyses showed that these NRGs were mostly associated with the regulation of necroptosis, Influenza A, apoptosis, NOD-like receptor, and NF-Kappa B signaling pathway. Next, univariate and multivariate Cox regression and LASSO analysis were used to identify the correlation between NRGs and the prognosis of ESCA. We constructed a prognostic model to predict the prognosis of ESCA based on SLC25A5, PPIA, and TNFRSF10B; the model classified patients into high- and low-risk subgroups based on the patient’s risk score. Furthermore, the receiver operating characteristic (ROC) curve was plotted, and the model was affirmed to perform moderately well for prognostic predictions. In addition, Gene Expression Omnibus (GEO) datasets were selected to validate the applicability and prognostic value of our predictive model. Based on different clinical variables, we compared the risk scores between the subgroups of different clinical features. We also analyzed the predictive value of this model for drug sensitivity. Moreover, Immunohistochemical (IHC) validation experiments explored that these three NRGs were expressed significantly higher in ESCA tissues than in adjacent non-tumor tissues. In addition, a significant correlation was observed between the three NRGs and immune-cell infiltration and immune checkpoints in ESCA.
In summary, we successfully constructed and validated a novel necroptosis-related signature containing three genes (SLC25A5, PPIA, and TNFRSF10B) for predicting prognosis in patients with ESCA; these three genes might also play a crucial role in the progression and immune microenvironment of ESCA.
An estimated 604,100 new cases of esophageal cancers (ESCA) were diagnosed, and 544,076 related deaths were observed worldwide in 2020; furthermore, ESCA ranks seventh most commonly diagnosed cancer and sixth in cancer-related mortality . Despite great curative progress in surgical treatment, chemoradiotherapy and systemic treatment, which has largely improved the survival rate of tumor patients, the prognosis of ESCA remains poor . This to a large extent is because the condition is mostly diagnosed at an advanced stage and symptoms are non-specific in the early stages, with the 5-year survival rate being lower than 15% . Thus, there is a need to identify more predictive and diagnostic tumor biomarkers to achieve a good outcome in ESCA patients. Therefore, finding more marker-related genes predicting the prognosis of patients could be an important strategy for the treatment of ESCA.
Currently, most anticancer drugs mainly trigger cancer cell death by induction of apoptosis . However, apoptosis resistance makes cancer cells drug-resistant and limits the efficacy of therapies. With the in-depth study of the mechanism of cell death, more cell death forms are being continuously discovered. Necroptosis, a newly recognized mechanism of programmed necrosis, has similar characteristics to both necrosis and apoptosis . Morphological features of necroptosis include cell and organelle swelling, plasma membrane rupture, and organelle breakdown, leading to the leakage of intracellular contents and subsequently inducing an inflammatory response . Necroptosis is triggered through the activation of cellular receptors, mainly including tumor necrosis receptor (TNFR1), death receptors (Fas/FasL), Toll-like receptors (TLRs), and cytosolic nucleic acid sensors . Necroptosis generally relies on the kinase activity of receptor interacting proteins (RIPs) . Furthermore, the RIPK1/RIPK3/MLKL pathway appears to be a key pathway mediating necroptosis.
There exists a complex relationship between necroptosis and cancer. Necroptosis affects tumorigenesis, infiltration, and metastasis, thereby affecting cancer prognosis. Accumulating evidence suggested that necroptosis not only mediates physiological regulation but also results in inflammatory pathologies and various cancers such as gastric cancer, colorectal cancer, cholangiocarcinoma, pancreatic cancer, and hepatocellular carcinoma [9,10,11,12,13]. ZENG et al.  confirmed that RIPK1 expression level was significantly upregulated in colorectal adenocarcinoma samples, while that of RIPK3 and p-MLKL were downregulated, suggesting that these cancer cells escape necroptosis for their survival. In cholangiocarcinoma, matrine induces necroptosis by enhancing RIP3 expression, producing reactive oxygen species (ROS), and activating the downstream RIP3/MLKL/ROS signaling pathway . Furthermore, apoptosis-resistant tumors have been shown to be regressed by inducing necroptosis . Necroptosis, especially necroptosis-related genes (NRGs), plays an important role in cancer. However, to date, the role of NRGs in ESCA remains unclear. This study, therefore, investigates the relationship between NRGs and ESCA to contribute to the diagnosis and prognosis of the disease.
In this study, we used the TCGA database to explore the relevance of NRGs in predicting the prognosis of ESCA. We thus constructed a necroptosis risk-scoring prognosis signature based on the screened prognostic NRGs. Subsequently, the applicability and prognostic value of the predictive model were validated using a GEO ESCA database. We also investigated the expression levels of the candidate NRGs in 20 ESCA tissues and matched adjacent normal tissues by immunohistochemistry (IHC). In addition, we examined the relationship between NRGs and the immune microenvironment in ESCA. This study may provide additional information regarding diagnosis and prognostic biomarkers for ESCA.
Expression levels of NRGs in ESCA
We acquired 6000 differentially expressed genes (DEGs) between 162 ESCA and 1456 normal tissues obtained from the UCSC Xena database and drew a heat map based on the expression level of each gene (Fig. 1A). Of these 6000 DEGs, 4820 were upregulated and 1180 were downregulated, as revealed by a volcano plot (Fig. 1B). The Gene ontology (GO) functional enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were summarized to clarify the biological significance of DEGs (Fig. 1C, D). The cut-off criteria for DEG were |Log2-fold change |> 1 and adjusted P-values < 0.05. We retrieved 159 NRGs from KEGG (Additional file 1: Table S1), which were further analyzed. Of these 159, 52 NRGs were identified among the 6000 DEGs (Additional file 2: Table S2, Fig. 2A, B). More specifically, of the 52 NRGs that were differentially expressed, 45 were upregulated and 7 were downregulated in ESCA samples compared with their expression in normal samples (Fig. 2C–E).
Genetic variation and expression of NRGs
To evaluate the role of NRGs in tumorigenesis, we summarized the mutation frequency and variant classification of the 52 NRGs screened in ESCA samples. As shown in Fig. 2F, G, genetic mutations of NRGs were observed in 19.02% (35 / 184) of the ESCA samples. The most common of all mutations were missense mutations (Fig. 2F). SNPs were the most common type of variation, and C > T ranked as the highest SNV class. Among the 52 NRGs, GLUD2, PYGM, JAK3, and STAT1 were the genes with the highest mutation frequency (Fig. 2F, G).
Functional enrichment analysis of NRGs
The GO functional enrichment and KEGG pathway analyses were performed to clarify the biological significance of NRGs. The enrichment analysis results of 52 NRGs are summarized in Figs. 3 and 4. GO terms are grouped into three major categories: biological processes (BP), cellular component (CC), and molecular function (MF). We found that the most significantly enriched BP terms were regulation of response to cytokine stimulus, regulation of cytokine-mediated signaling pathway, l-kappaB kinase/NF-kappaB signaling, cysteine-type endopeptidase activity, apoptotic process, and necrotic cell death; enriched CC terms were membrane region, membrane raft, membrane microdomain, cytosolic part, and ESCRT III complex; and enriched MF terms were cytokine receptor binding, ubiquitin-like protein ligase binding, tumor necrosis factor receptor superfamily binding, tumor necrosis factor receptor binding, cytokine receptor binding, and cysteine-type endopeptidase regulator activity involved in the apoptotic process (Fig. 3A). In terms of KEGG pathway analysis, these NRGs were mainly associated with pathways related to necroptosis, influenza A, NOD-like receptor signaling pathways, apoptosis, and Measles (Fig. 4A). Moreover, a network diagram of NRGs was plotted to represent the relationship between items and molecules (Figs. 3B, 4B). We then drew a circle of enriched GO and KEGG terms and combined them with Z-Score to predict the function of 52 NRGs in these pathways (Figs. 3C–E, 4C–E).
Establishment of a necroptosis-related prognostic gene model for ESCA
To establish a prognostic risk signature, we performed a univariate Cox regression analysis of the 52 screened NRGs. As a result, three NRGs, namely SLC25A5, PPIA, and TNFRSF10B, were identified as genes with significant prognostic value in ESCA (Fig. 5A–C, Additional file 3: Fig. S1). The results of the Kaplan–Meier survival (KM) analysis (Fig. 5) suggested that high expression of SLC25A5 was correlated with worse prognosis in ESCA patients (overall survival (OS), p = 0.002; progression-free survival (PFS), p = 0.014; disease-specific survival (DSS), p = 0.007) (Fig. 5D–F). Moreover, PPIA high expression was also associated with poor prognosis in ESCA patients (OS, p = 0.045; DSS, p = 0.048) (Fig. 5H, I). However, high TNFRSF10B expression resulted in a better prognosis in ESCA (OS, p = 0.042) (Fig. 5G). Furthermore, correlation analyses showed that PPIA was positively correlated with SLC25A5 (Fig. 5J). The effect of interaction between PPIA and SLC25A5 may be worthy of further study.
After incorporating the result of the LASSO regression analysis, the corresponding three genes were selected for the signature; the model fitted the data well, and the penalty coefficient was three (Fig. 6A, B). Subsequently, we performed multivariate Cox regression analysis of the three NRGs. The results showed that these three NRGs could act as prognostic predictors when coupled with the beta value of the multivariate Cox regression. The risk score was calculated as, risk score = (0.3481) * SLC25A5 + (− 0.0405) * PPIA + (− 0.0948) * TNFRSF10B. Based on the risk score, ESCA patients were classified into a high-risk group and a low-risk group. Compared to the low-risk group, the high-risk group had higher mortality and worse prognosis (Fig. 6C). As shown in Fig. 6D, ESCA patients belonging to the high-risk group had a higher probability of death than those belonging to the low-risk group (median time = 1.5 years vs. 3.8 years, p = 0.0000198). The AUCs of the risk assessment model for the three NRGs were 0.568, 0.743, and 0.877 at 1-year, 3-year, and 5-year, respectively (Fig. 6E). The model exhibited good accuracy of prediction.
Building and validation of the predictive nomogram
Based on the three NRGs and the clinical factors (TNM stage, age, and gender), we constructed a nomogram to predict the survival probability of ESCA patients. Univariate analysis showed that N1-3, M1, and the expression of SLC25A5 and TNFRSF10B were associated with the prognosis of ESCA patients (Additional file 4: Table S3). Multivariate Cox analysis indicated that N1-3 and M1 were independent predictors of prognosis (Additional file 4: Table S3). Based on the Cox regression algorithm, we produced a nomogram to predict the 1-year, 2-year, and 5-year OS rates (Fig. 7A). The calibration plots showed that the nomogram had the best prediction accuracy for 1- year and 3-year OS rates in the entire cohort (Fig. 7B).
Validation of the predictive model
To verify the applicability and prognostic value of our predictive model, we randomly divided the ESCA data set from TCGA into two equal groups (named random set1 and random set2). We used the TCGA set, random set1, random set2, and the external dataset GSE53624 for assessing the predictive ability of our model. Considering the median risk score as the cutoff value, ESCA patients were classified into a high-risk group and a low-risk group. Figure 8A–H shows the distribution of risk scores and OS status of the four sets. The heat map indicated the association between the three NRGs and the risk scores (Fig. 8I–L). The expression of the three genes differed between the tissues of high-risk and low-risk groups. Furthermore, the receiver operating characteristic (ROC) curve proved the high prediction accuracy of the model for 1-, 3-, and 5-year survival in the four datasets (AUC values were > 0.5 for all). It is noteworthy that the AUC values of the model in the whole set, random set1, and random set2 for predicting 5-year survival were greater than 0.75 (Fig. 8). The prediction performance of the model in the external dataset GSE53624 was also better (Fig. 8). In addition, the results of the Kaplan–Meier plot revealed that the high-risk group had a worse OS compared with that of the low-risk group in the whole set (p = 0.0096) (Fig. 8Q) and random set1 (p = 0.0044) (Fig. 8R). Moreover, the assessment of the correlation between prognosis and the clinical factors (T stage, N stage, stage, age, gender, and risk score) by univariate Cox and multivariate Cox regression analysis revealed that stage I and different risk scores were associated with prognosis of ESCA patients (p < 0.05) in the whole set and random set1 (Fig. 8U, V). However, these findings were not the same as the results obtained for the random set2 and GSE53624 dataset (p > 0.05) (Fig. 8W, X). Overall, the model had a good predictive effect.
Subgroups analysis of clinical features in the predictive model
Next, we compared the risk scores among the subgroups of different clinical features (T stage, N stage, M stage, tumor stage, gender, and survival). Clinical features and risks for each patient are summarized in the heat map (Fig. 9A). Unfortunately, no difference between the risker score assessment between the subgroups was observed (p > 0.05) (Fig. 9B–G). Meanwhile, we found that the model exhibited good OS predictive ability in the male subgroup and stage I–II subgroup (p < 0.05) (Fig. 9H–Q).
Correlation between prognostic model and drug sensitivity
A combination of targeted drugs and chemotherapy is commonly used in the treatment of advanced ESCA. Consequently, we assessed the difference in efficacy of targeted drugs and chemotherapy drugs between the high-risk group and low-risk group. We predicted the efficacy of five chemotherapeutic drugs (Cisplatin, Paclitaxel, Docetaxel, Gemcitabine, and Methotrexate) and BIBW2992 against ESCA. The results indicated that the low-risk group was more sensitive to Paclitaxel and BIBW2992 than the high-risk group (p < 0.05) (Fig. 9R).
Protein expression analysis of the three NRGs in ESCA
We employed IHC to validate the expression of the proteins of the three NRGs in 20 pairs of ESCA tumor tissues and corresponding adjacent normal tissues. According to IHC staining analysis, the protein products of the three NRGs were located primarily in the cytoplasm, membrane, and mitochondria of cancer cells, with brown staining reflecting positive staining (Fig. 10A, C, E, G, I, K). In normal tissues, the proteins of three NRGs were only weakly expressed or not expressed (Fig. 10B, D, F, H, J, L). Based on statistical analysis, these three NRGs were expressed significantly higher in ESCA tissues than in adjacent non-tumor tissues (P < 0.001) (Fig. 10M–O).
NRGs are associated with tumor immune infiltration and immune checkpoints in ESCA
Tumor-infiltrating lymphocyte levels are an independent predictor of survival in cancers. Here, we proved that the expression of prognostic NRGs (PPIA, SLC25A5, and TNFRSF10B) correlated with immune infiltration in ESCA using the TCGA dataset. According to the CIBERSORT algorithm, we found that PPIA expression was significantly negatively associated with resting memory CD4 + T cell, activated mast cell, and memory B cell (p < 0.05), whereas it was significantly positively correlated with resting mast cell and macrophage M0 (p < 0.05) (Fig. 11A). SLC25A5 expression was significantly positively correlated with resting myeloid dendritic cell and activated mast cell (p < 0.05) (Fig. 11A). Meanwhile, TNFRSF10B expression was significantly positively correlated with resting memory CD4 + T cell, activated memory CD4 + T cell, neutrophil, and resting NK cell, and was significantly negatively correlated with resting myeloid dendritic cell and macrophage M2 (p < 0.05) (Fig. 11A).
To further confirm the results, we investigated the correlation between these three NRGs and immune infiltration cells in ESCA via the TIMER database. The results revealed a negative correlation between PPIA expression and the infiltration of B cells (p = 0.0284), CD8 + T cells (p = 0.0364), macrophages (p = 0.0118), and neutrophils (p = 0.0000982), whereas there was no correlation with tumor purity, CD4 + T cells, and dendritic cells (Fig. 11B). Furthermore, there was a negative correlation between SLC25A5 expression and the immune infiltration level of dendritic cells (Fig. 11B, p = 0.00989). However, SLC25A5 expression was not associated with B cells, CD4 + T cells, CD8 + T cells, neutrophils, and macrophages. Moreover, TNFRSF10B expression was positively associated with the infiltration of CD8 + T cells (Fig. 11B, p = 0.000435), while there was no association between TNFRSF10B levels and tumor purity, B cell, CD4 + T cell, macrophages, neutrophils, and dendritic cells (Fig. 11B). We further analyzed the correlation between the risk score of NRGs and immune infiltration in ESCA. As expected, the risk score of NRGs was significantly negatively correlated with the abundance of CD8 + T cells (p = 3.99e−4), neutrophils (p = 3.07e− 4), and myeloid dendritic cells (p = 1.141e −4) (Fig. 12A).
SIGLEC15, PDCD1LG2(PD-L2), TIGIT, PDCD1(PD-1), CD274(PD-L1), CTLA4, LAG3, and HAVCR2(TIM3) are immunological checkpoints that perform a vital function in tumor immune evasion. Considering that these three NRGs might be the predictive biomarkers in ESCA, the relationship of these NRGs with the above-mentioned checkpoints was assessed. Notably, the expression of TNFRSF10B was positively correlated with that of CTLA4 and SIGLEC15 (p < 0.05), while PPIA expression was negatively correlated with that of CD274, CTLA4, PDCD1, SIGLEC15, and TIGIT in ESCA (p < 0.05) (Fig. 12B).
Necroptosis is a newly discovered mechanism of programmed necrosis that plays an essential function in various cancers. Necroptosis is also the cell death mechanism in tumor cells resistant to apoptosis and is triggered by a multitude of different stimuli. The classical necroptosis pathway mediated by kinase-dependent RIPs is triggered when pro-apoptotic molecules fail to stimulate apoptotic bodies . Tumor multi-drug resistance characterized by apoptosis limits the clinical application of apoptosis inducers. Therefore, targeting necroptosis may be a novel strategy to bypass apoptotic resistance and eliminate cancer cells . Furthermore, the role of inflammation and immunity in necroptosis cannot be ignored. Biomarkers of necroptosis are significantly associated with the immunological profile in ESCA . It has been reported that necroptosis induced by chemotherapy can stimulate inflammatory responses and elicit immunogenic and anticancer effects . Necrotic tumor cells can hinder anti-tumor immunity; these cells, which are cleared by monocytes, neutrophils, and macrophages, can induce the release of inflammatory factors and ultimately trigger an adaptive immune response . The critical role of NRGs in the pathogenesis of cancer makes them a potential prognostic and therapeutic target in cancer. However, the clinical significance of NRGs in ESCA has not been elucidated to date and further study is still needed.
Thus, to investigate the role of NRGs in ESCA, we used the TCGA database to analyze the expression of 52 NRGs in ESCA and normal esophageal tissues. GO and KEGG functional enrichment analyses showed that these NRGs were mostly associated with the regulation of necroptosis, Influenza A, apoptosis, NOD-like receptor signaling pathway, and NF-Kappa B signaling pathway. In addition, immune-related functions were also enriched, such as cytokine-cytokine receptor interaction, T cell activation, and positive regulation of cell adhesion. As expected, the enriched functions were associated with the tumor immune microenvironment, carcinogenesis, progression, and necroptosis in ESCA. Preliminary research has indicated that NOD-like receptor signaling pathways are possible predictors of esophageal adenocarcinoma .
Furthermore, our results indicated that three differentially expressed NRGs (SLC25A5, PPIA, and TNFRSF10B) were significantly related to the prognosis of ESCA. TNFRSF10B was found to be a favorable prognostic gene, while SLC25A5 and PPIA were related to adverse prognosis in ESCA. We then performed LASSO Cox regression to construct a prognostic signature based on the three prognostic NRGs. According to the risk score, ESCA patients were divided into two groups (high-risk and low-risk groups), and the results revealed that the high-risk group had a significantly poorer OS than the low-risk group. Univariate analysis indicated that N1-3, M1, and high expression of SLC25A5 and TNFRSF10B were correlated with the prognosis of ESCA patients (Additional file 4: Table S3). Furthermore, according to Multivariate Cox analysis, N1-3 and M1 were independent predictors of prognosis. The ROC curve confirmed the prognostic signature to be an independent prognostic indicator by showing that it could predict the 1-year, 3-year, and 5-year OS rates relatively well compared with an ideal model in the entire cohort.
To test the applicability and prognostic value of our model, we randomly divided the ESCA data set from TCGA into random set1 and random set2. We used random set1, random set2, and the external dataset GSE53624 for assessing the predictive ability of our model. The model exhibited a good predictive effect in all three datasets. Subgroups analysis of clinical features in the predictive model indicated that the model had good OS predictive ability in the male subgroup and stage I-II subgroup. Correlation between the prognostic model and clinical treatment indicated that the low-risk group was more sensitive to Paclitaxel and BIBW2992 than the high-risk group. Subsequently, IHC staining analysis showed the three NRGs were expressed significantly higher in ESCA tissues than in adjacent non-tumor tissues. In summary, we, for the first time, constructed a necroptosis‑related prognostic gene signature for ESCA, which provides new clues for prognostic prediction in ESCA patients.
TNFRSF10B is a member of the TNF-receptor superfamily that is activated by TNF-related apoptosis-inducing ligand (TRAIL). TNFRSF10B/DR5 delivers apoptotic signals to the cell and induces apoptosis in cancer cells . The expression level of TNFRSF10B was associated with tumor progression and apoptosis . Furthermore, the high expression of DR5 mediated the extrinsic apoptotic pathway in various cancer cells . However, He et al. showed that the low expression of TNFRSF10B was associated with a poor prognosis in esophageal squamous cell carcinoma [26, 27], which is consistent with our findings. SLC25A5 (ANT2) is a by-product of nucleotide transferase which is specifically expressed in proliferating cells and participates in glycolytic metabolism. Hence, SLC25A5 is associated with cell growth and differentiation . Depletion of SLC25A5 can cause mitochondrial dysfunction and induce oxidative stress, leading to erythrocyte anemia and B-cell depletion . SLC25A5 promotes apoptosis through the regulation of bcl-2, caspase-3, and bax in prostate cancer . Studies have shown that high expression of SLC25A5 in cervical cancer could be an independent prognostic factor . PPIA (Cyclophilin A) is a member of the immunophilin family and acts as an immune inflammatory mediator that secretes oxidative stress-induced, which promotes the formation of foam cells by increasing the levels of ROS and pro-inflammatory cytokines. Furthermore, PPIA is involved in biological processes such as intracellular signaling, transcription, and apoptosis, therefore playing critical roles in microorganismal infections, inflammatory diseases, and tumor proliferation [32,33,34]. An increase in PPIA levels may lead to macrophage apoptosis through activation of mitochondrial death signaling pathways and caspase 3 cascade . Studies have confirmed that the upregulation of PPIA is associated with a poor OS in diseases such as lung adenocarcinoma . Some researchers have indicated that overexpression of PPIA is associated with decreased survival in esophageal squamous cell carcinoma and shown it to be an independent prognostic factor [37, 38]. Although these studies reveal the relationship between the NRGs and tumor progression, none of these studies have explored the expression and function of the NRGs in ESCA.
The tumor immune microenvironment is mainly composed of tumor-infiltrating lymphocytes and other immune cells such as dendritic cells, neutrophils, and macrophages. Tumor-infiltrating lymphocytes can inhibit or promote tumor progression . Our study confirmed that the expression levels of the three NRGs were significantly associated with immune cell infiltration in ESCA. Early findings suggested that abundant CD8+ T cells and CD4+ T cells could be prognostic indicators for the clinical outcome in cancers [40, 41]. Our results are consistent with the previous reports suggesting that TNFRSF10B expression associated with high CD8+ T cell abundance may indicate better clinical outcomes. Preliminary research revealed that the tumor microenvironment is associated with tumor growth and prognosis in ESCA, and increasing levels of immune infiltrates reduce the risk of distant metastasis and death [27, 42]. Our results showed that the expression levels of the three prognostic NRGs were significantly correlated with that of immunological checkpoints in ESCA. Our research might provide more clues for immunotherapy strategies in ESCA. All these outcomes illustrate that tumor immune evasion and antitumor immunity might be implicated in the three prognostic NRGs-mediated carcinogenic processes in ESCA.
It was interesting to note that different approaches yielded inconsistent results regarding the relationship between immune cells infiltration and the three prognostic NRGs. This inconsistency may be attributable to the following reasons. Although flow cytometry, immunohistochemistry staining, or single-cell sequencing can be used to estimate the immune cell status in a tumor sample, each has limitations that prevent them from being widely applied. Therefore, we used computational methods to evaluate immune-cell composition from bulk RNA-sequencing data. First of all, between the the actual situation and computer-based algorithms, there were some variations. Next, the mechanisms of immune cell infiltration in tumors are complex, and they are inevitably affected by intratumorally heterogeneity and small sample sizes. In the end, various algorithms are used in these methods, and they all have advantages and disadvantages.
Although the prognostic model exhibited good performance in the TCGA database and external datasets, our study still has some limitations in lacking mechanistic experiment and further experimental validation are needed.
In summary, we identified three prognosis-associated NRGs (SLC25A5, PPIA, and TNFRSF10B) and constructed a novel necroptosis-related prognostic gene signature. Furthermore, in this study, we demonstrated that SLC25A5, PPIA, and TNFRSF10B were correlated with tumor immune microenvironment infiltration and may be potential prognostic biomarkers for ESCA. The prognostic model exhibited good performance for the prediction of ESCA prognosis in TCGA database and external datasets. Our findings will benefit the treatment and diagnosis of ESCA. Nevertheless, the study results need to be validated in future fundamental research and extensive clinical trials.
Materials and methods
Datasets and data processing
The UCSC Xena database was utilized to obtain TCGA-ESCA RNA-Seq FPKM data and clinical information and survival information (https://toil-xena-hub.s3.us-east-1.amazonaws.com/download/TcgaTargetGtex_rsem_gene_tpm.gz; Full metadata) . One independent ESCA gene expression profiles (GSE53624) were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) and processed for analysis . There are 162 ESCA and 11 normal cases with clinical data were downloaded from the database (Table 1). 159 NRGs were downloaded from Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Additional file 1: Table S1) [45,46,47]. All statistical analyses were carried out using R version v4.0.3.
NRGs differential expression analysis
A total of 52 NRGs were extracted and showed in Additional file 2: Table S2, which was the most representative significance genes. The "limma" software packages were used to identify the differences in NRG expression between ESCA and normal esophageal tissues, and |Log2-fold change |> 1 and p-values < 0.05 were set as the filter conditions. Analyzing the mutation rates of 52 NRGs in ESCA patients and mapping waterfall and frequency plots of 52 NRGs in ESCA patients were generated using the "maftools" software package .
Functional enrichment analysis
Gene Ontology (GO) analyses, including biological processes, molecular functions, and cellular components, and KEGG pathway enrichment analysis were used to identify characteristic biological attributes via the R package “ClusterProfiler ”. Data visualization was carried out using R package "ggplot2 ". p < 0.05 was considered to be statistically significant.
Establishment of necroptosis risk scoring prognosis signature
The correlation between NRGs expression and survival was analyzed using univariate survival analysis. The Kaplan–Meier  was used to compare the survival rates of different NRGs expression. The risk score was calculated from the formula of the regression coefficients in the multiple Cox regression model. Thereafter we established the prognostic model for three prognostic NRGs by LASSO Cox regression analysis . According to the median risk score, ESCA patients were divided into two subgroups (high-risk and low-risk), then we compared OS between two subgroups by KM analysis. Time ROC analysis was employed to predict accuracy of the risk score. Based on the clinical characteristics (TNM stage, gender, age and three NRGs), we constructed a nomogram to quantitatively predict 1-, 3-, and 5-year OS. Each variable such as p-value, HR and 95% CI were visualized by a forest by R package "forestplot".
Validation of predictive model
We used TCGA ESCA database with external ESCA dataset GSE53624 to assess and verify the predictive ability of our prognostic model [53, 54]. The Clinical and survival information of the whole set was derived from the TCGA dataset, and the external dataset GSE53624 was obtained from the GEO database through the GEO query R package. TCGA was equally divided into two groups by random sampling, named Random set1, 2. Time ROC analysis was employed to predict accuracy of 1-, 3-, and 5-year survival. The correlation between risk score and survival time was examined via KM survival analysis. Prognostic and predictive value of Risk score was validated by Univariate Cox analysis. All analytical analyses were performed using the R packages "ggplot2", "coxph function", "pheatmap", "timeROC", "survival", and "survminer".
Correlation between prognostic model and drug sensitivity
We used the R package “pRRophetic” to predict drug sensitivity and the results were visualized by the expression matrix.
In order to further validate the expression of the three NRGs between ESCA and adjacent normal tissues, 20 paired ESCA tissues from Liuzhou People's Hospital were collected for immunohistochemical staining. The study was approved by the Ethics Committee of Liuzhou People's Hospital (Reference No. KY2021-026-01) and conducted according to the Declaration of Helsinki. Formalin fixed paraffin-embedded tissue were analyzed by Immunohistochemistry with PPIA antibody (1:100; Proteintech, China), SLC25A5(1:100; Proteintech, China), and TNFRSF10B (1:150; OriGene China) and horseradish peroxidase conjugated secondary antibodies (Maxim, China). For IHC quantification, the integrated optical density (IOD) for each slice was calculated using the Image-ProPlus6.0 software (Media Cybernetics, USA).
Tumor microenvironment estimation
We used CIBERSORT algorithm to confirm the association between three NRGs and the abundance of 22 infiltrating immune cells. Then we used spearman correlation analysis to further confirm the association between the risk scores and those NRGs expression. This study also explored the correlation between three NRGs and 8 immune checkpoint molecules. All statistical analyses information mentioned were visualized via R version 4.0.3.
Tumor immune estimation resource (TIMER) database
TIMER (https://cistrome.shinyapps.io/timer/) dataset comprise six tumor-infiltrating immune subsets. Calculated the levels of six subgroups for 10,897 tumors in 32 cancers using the TCGA. Based on the database, gene expression and tumor immune infiltration (B cells, CD4+ T cells, CD8+ T cells, Dendritic cells, Macrophages, and Neutrophils) were analyzed in several cancer types. Using the TIMER dataset, we examined the mRNA expression of these prognostic NRGs in patients with ESCA.
All Statistical analysis analyzed were conducted by the Log-rank test, such as fold-change (FC), HR, and p-values. The correlation between particular variables was validated using the Spearman's correlation analysis or Pearson correlation analysis. p-value or Log-rank p-value of < 0.05 was considered as having statistical significance.
Availability of data and materials
The datasets generated and analysed during the current study are available in the [UCSC Xena database] repository, https://xenabrowser.net/.
Differentially expressed genes
- UCSC Xena:
University Of Cingifornia Sisha Cruz
The Cancer Genome Atlas
Kyoto Encyclopedia of Genes and Genomes
Least absolute shrinkage and selection operator
Receiver operating characteristic
Tumor necrosis receptor 1
Receptor interacting proteins
Receptor-interacting serine/threonine protein kinase 1/3
Mixed lineage kinase domain-like protein
Reactive oxygen species
Tumor Immune Estimation Resource
- KM plotter:
Search Tool for the Retrieval of Interacting Genes
Area under the curves
- NK cell:
Natural killer cell
Sung H, Ferlay J, Siegel RL, et al. Global Cancer Statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2021;71(3):209–49.
Shah MA, Kennedy EB, Catenacci DV, et al. Treatment of locally advanced esophageal carcinoma: ASCO guideline. J Clin Oncol. 2020;38(23):2677–94.
di Pietro M, Canto MI, Fitzgerald RC. Endoscopic management of early adenocarcinoma and squamous cell carcinoma of the esophagus: screening, diagnosis, and therapy. Gastroenterology. 2018;154(2):421–36.
Pistritto G, Trisciuoglio D, Ceci C, Garufi A, D’Orazi G. Apoptosis as anticancer mechanism: function and dysfunction of its modulators and targeted therapeutic strategies. Aging (Albany NY). 2016;8(4):603–19.
Morrice JR, Gregory-Evans CY, Shaw CA. Necroptosis in amyotrophic lateral sclerosis and other neurological disorders. Biochim Biophys Acta Mol Basis Dis. 2017;1863(2):347–53.
Chen J, Kos R, Garssen J, Redegeld F. Molecular insights into the mechanism of necroptosis: the necrosome as a potential therapeutic target. Cells. 2019;8(12):1486.
Bertheloot D, Latz E, Franklin BS. Necroptosis, pyroptosis and apoptosis: an intricate game of cell death. Cell Mol Immunol. 2021;18(5):1106–21.
Saeed WK, Jun DW, Jang K, Koh DH. Necroptosis signaling in liver diseases: an update. Pharmacol Res. 2019;148:104439.
Guo D, Zhang W, Yang H, et al. Celastrol induces necroptosis and ameliorates inflammation via targeting biglycan in human gastric carcinoma. Int J Mol Sci. 2019;20(22):5716.
Ando Y, Ohuchida K, Otsubo Y, et al. Necroptosis in pancreatic cancer promotes cancer cell migration and invasion by release of CXCL5. PLoS ONE. 2020;15(1): e0228015.
Lomphithak T, Akara-Amornthum P, Murakami K, et al. Tumor necroptosis is correlated with a favorable immune cell signature and programmed death-ligand 1 expression in cholangiocarcinoma. Sci Rep. 2021;11(1):11743.
Xiang YK, Peng FH, Guo YQ, et al. Connexin32 activates necroptosis through Src-mediated inhibition of caspase 8 in hepatocellular carcinoma. Cancer Sci. 2021;112(9):3507–19.
Yan S, Li Q, Zhang D, et al. Necroptosis pathway blockage attenuates PFKFB3 inhibitor-induced cell viability loss and genome instability in colorectal cancer cells. Am J Cancer Res. 2021;11(5):2062–80.
Zeng F, Chen X, Cui W, et al. RIPK1 binds MCU to mediate induction of mitochondrial Ca(2+) uptake and promotes colorectal oncogenesis. Cancer Res. 2018;78(11):2876–85.
Xu B, Xu M, Tian Y, et al. Matrine induces RIP3-dependent necroptosis in cholangiocarcinoma cells. Cell Death Discov. 2017;3:16096.
He GW, Günther C, Thonn V, et al. Regression of apoptosis-resistant colorectal tumors by induction of necroptosis in mice. J Exp Med. 2017;214(6):1655–62.
Kreuzaler P, Watson CJ. Killing a cancer: what are the alternatives. Nat Rev Cancer. 2012;12(6):411–24.
Philipp S, Sosna J, Adam D. Cancer and necroptosis: friend or foe. Cell Mol Life Sci. 2016;73(11–12):2183–93.
Yamauchi T, Fujishima F, Hashimoto M, et al. Necroptosis in esophageal squamous cell carcinoma: an independent prognostic factor and its correlation with tumor-infiltrating lymphocytes. Cancers (Basel). 2021;13(17):4473.
Kroemer G, Galluzzi L, Kepp O, Zitvogel L. Immunogenic cell death in cancer therapy. Annu Rev Immunol. 2013;31:51–72.
Meng MB, Wang HH, Cui YL, et al. Necroptosis in tumorigenesis, activation of anti-tumor immunity, and cancer therapy. Oncotarget. 2016;7(35):57391–413.
Liu D. LYN, a key gene from bioinformatics analysis, contributes to development and progression of esophageal adenocarcinoma. Med Sci Monit Basic Res. 2015;21:253–61.
Wang S, El-Deiry WS. TRAIL and apoptosis induction by TNF-family death receptors. Oncogene. 2003;22(53):8628–33.
Hernandez-Cueto A, Hernandez-Cueto D, Antonio-Andres G, et al. Death receptor 5 expression is inversely correlated with prostate cancer progression. Mol Med Rep. 2014;10(5):2279–86.
Xiao Z, Nie K, Han T, et al. Development and validation of a TNF family-based signature for predicting prognosis, tumor immune characteristics, and immunotherapy response in colorectal cancer patients. J Immunol Res. 2021;2021:6439975.
He SY, Wang XB, Jiao YC. Data mining of esophageal squamous cell carcinoma from The Cancer Genome Atlas database. Zhonghua Zhong Liu Za Zhi. 2018;40(7):517–22.
Lu G, Chen L, Wu S, Feng Y, Lin T. Comprehensive analysis of tumor-infiltrating immune cells and relevant therapeutic strategy in esophageal cancer. Dis Markers. 2020;2020:8974793.
Han Y, Bi Y, Bi H, et al. miR-137 suppresses the invasion and procedure of EMT of human breast cancer cell line MCF-7 through targeting CtBP1. Hum Cell. 2016;29(1):30–6.
Chevrollier A, Loiseau D, Stepien G. What is the specific role of ANT2 in cancer cells? Med Sci (Paris). 2005;21(2):156–61.
Cho J, Seo J, Lim CH, et al. Mitochondrial ATP transporter Ant2 depletion impairs erythropoiesis and B lymphopoiesis. Cell Death Differ. 2015;22(9):1437–50.
Zhang H, Chen N, Deng Z, et al. Suppression of ANT2 by miR-137 inhibits prostate tumorigenesis. Front Genet. 2021;12:687236.
Nigro P, Pompilio G, Capogrossi MC. Cyclophilin A: a key player for human disease. Cell Death Dis. 2013;4(10):e888.
Li S, Han F, Qi N, et al. Determination of a six-gene prognostic model for cervical cancer based on WGCNA combined with LASSO and Cox-PH analysis. World J Surg Oncol. 2021;19(1):277.
Wu Y, Ma Z, Zhang Y, et al. Cyclophilin A regulates the apoptosis of A549 cells by stabilizing Twist1 protein. J Cell Sci. 2022;135(2):jcs259018.
Volker SE, Hedrick SE, Feeney YB, Clevenger CV. Cyclophilin A function in mammary epithelium impacts Jak2/Stat5 signaling, morphogenesis, differentiation, and tumorigenesis in the mammary gland. Cancer Res. 2018;78(14):3877–87.
Anandan V, Thankayyan Retnabai SK, Jaleel A, et al. Cyclophilin A induces macrophage apoptosis and enhances atherosclerotic lesions in high-fat diet-fed hyperglycemic rabbits. FASEB Bioadv. 2021;3(5):305–22.
Li Y, Guo H, Dong D, Wu H, Li E. Expression and prognostic relevance of cyclophilin A and matrix metalloproteinase 9 in esophageal squamous cell carcinoma. Diagn Pathol. 2013;8:207.
Wang S, Li M, Xing L, Yu J. High expression level of peptidylprolyl isomerase A is correlated with poor prognosis of liver hepatocellular carcinoma. Oncol Lett. 2019;18(5):4691–702.
Ju Q, Li XM, Zhang H, Zhao YJ. BRCA1-associated protein is a potential prognostic biomarker and is correlated with immune infiltration in liver hepatocellular carcinoma: a pan-cancer analysis. Front Mol Biosci. 2020;7:573619.
Li Y, Lu Z, Che Y, et al. Immune signature profiling identified predictive and prognostic factors for esophageal squamous cell carcinoma. Oncoimmunology. 2017;6(11):e1356147.
Hu Z, Xie J, Chen X, Tang J, Zhou K, Han S. Identification of an immune-related biomarker model based on the CircRNA-associated regulatory network for esophageal carcinoma. J Oncol. 2021;2021:1334571.
Baba Y, Nomoto D, Okadome K, et al. Tumor immune microenvironment and immune checkpoint inhibitors in esophageal squamous cell carcinoma. Cancer Sci. 2020;111(9):3132–41.
Goldman MJ, Craft B, Hastie M, et al. Visualizing and interpreting cancer genomics data via the Xena platform. Nat Biotechnol. 2020;38(6):675–8.
Barrett T, Wilhite SE, Ledoux P, et al. NCBI GEO: archive for functional genomics data sets–update. Nucl Acids Res. 2013;41(Database):D991-5.
Kanehisa M, Goto S. KEGG: Kyoto encyclopedia of genes and genomes. Nucl Acids Res. 2000;28(1):27–30.
Kanehisa M. Toward understanding the origin and evolution of cellular organisms. Protein Sci. 2019;28(11):1947–51.
Kanehisa M, Furumichi M, Sato Y, Ishiguro-Watanabe M, Tanabe M. KEGG: integrating viruses and cellular organisms. Nucl Acids Res. 2021;49(D1):D545–51.
Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 2018;28(11):1747–56.
Wu T, Hu E, Xu S, et al. clusterProfiler 4.0: a universal enrichment tool for interpreting omics data. Innovation (Camb). 2021;2(3):100141.
Ito K, Murphy D. Application of ggplot2 to Pharmacometric graphics. CPT Pharmacometrics Syst Pharmacol. 2013;2(10):e79.
Nagy Á, Munkácsy G, Győrffy B. Pancancer survival analysis of cancer hallmark genes. Sci Rep. 2021;11(1):6047.
Simon N, Friedman J, Hastie T, Tibshirani R. Regularization paths for Cox’s proportional hazards model via coordinate descent. J Stat Softw. 2011;39(5):1–13.
Lu M, Li J, Fan X, Xie F, Fan J, Xiong Y. Novel immune-related ferroptosis signature in esophageal cancer: an informatics exploration of biological processes related to the TMEM161B-AS1/hsa-miR-27a-3p/GCH1 regulatory network. Front Genet. 2022;13:829384.
Fan X, Ou Y, Liu H, et al. A ferroptosis-related prognostic signature based on antitumor immunity and tumor protein p53 mutation exploration for guiding treatment in patients with head and neck squamous cell carcinoma. Front Genet. 2021;12:732211.
This work was supported by grants from Science and Technology Program of Liuzhou (2021CBB0104), the Science and Technology Base and Talent Project of Guangxi (2021AC19414), the Research Fund of Liuzhou People's Hospital (lry202108), the Talent Introduction Scientific Research Projects Funded Start-Up Funds of Liuzhou People's Hospital (LRYGCC202114) and Guangxi Natural Science Foundation (2019JJA140660). We thank Bullet Edits Limited for the linguistic editing and proofreading of the manuscript.
This work was supported by grants from Science and Technology Program of Liuzhou (2021CBB0104), the Science and Technology Base and Talent Project of Guangxi (2021AC19414), the Research Fund of Liuzhou People's Hospital (lry202108), the Talent Introduction Scientific Research Projects Funded Start-Up Funds of Liuzhou People's Hospital (LRYGCC202114) and Guangxi Natural Science Foundation (2019JJA140660).
Ethics approval and consent to participate
This study was approved by the Ethics Committee of Liuzhou People's Hospital (Reference No. KY2022-026-01), and was performed according to the Declaration of Helsinki. Written informed consent forms were obtained from all subjects.
Consent for publication
All authors are consent for the publication of this work.
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
. 159 NRGs from KEGG in EC.
. 52 NRGs among the DEGs between EC and normal samples.
. Univariate survival analysis of differentially expressed NRGs with prognostic value in ESCA.
. Analysis of factors affecting the prognosis of patients with esophageal carcinoma.
About this article
Cite this article
Sun, K., Hong, Jj., Chen, Dm. et al. Identification and validation of necroptosis-related prognostic gene signature and tumor immune microenvironment infiltration characterization in esophageal carcinoma. BMC Gastroenterol 22, 344 (2022). https://doi.org/10.1186/s12876-022-02423-6
- Esophageal carcinoma
- Immune infiltration