Stratification From Heterogeneity of the Cell-Death Signal Enables Prognosis Prediction and Immune Microenvironment Characterization in Esophageal Squamous Cell Carcinoma
Yiyuan Zhang1,2*† 
Yanxing Chen3†- 1Department of Radiation Oncology, Shandong First Medical University and Shandong Academy of Medical Sciences, Jinan, China
- 2Department of Radiation Oncology, Shandong Cancer Hospital and Institute, Shandong First Medical University and Shandong Academy of Medical Sciences, Jinan, China
- 3State Key Laboratory of Oncology in South China, Department of Medical Oncology, Sun Yat-sen University Cancer Center, Collaborative Innovation Center for Cancer Medicine, Guangzhou, China
Esophageal squamous cell carcinoma (ESCC) is the primary subtype of esophageal cancer (EC) characterized by a high incidence rate and extremely poor prognosis worldwide. Previous studies suggested that the specific cell death signal was linked to different immune subtypes in multiple cancers, while a comprehensive investigation on ESCC is to be performed yet. In the current study, we dissected different cell death signals in ESCC tumors and then integrated that functional information to stratify ESCC patients into different immunogenic cell death (ICD) subtypes. By systematically analyzing the transcriptomes of 857 patients and proteomic profile of 124 patients, we found that the signals of necroptosis, pyroptosis, and ferroptosis are positively associated with activated immunity in ESCC. We identified two ICD pattern terms, namely, ICD-high and ICD-low subtypes that positively correlated to both progression-free survival and overall survival. In addition, cell fraction deconvolution analysis revealed that more infiltrated leukocytes were enriched in ICD-high types, especially antigen-presenting cells, such as dendritic cells and macrophages. With the XGBoost algorithm, we further developed a 14-gene signature which can simplify the subtyping for allocating new samples, by which we validated the prognosis value of the signature and proved that the ICD score scheme could serve as a promising biomarker for stratifying patients with immunotherapy in several immune checkpoint blockade treatment cohorts. Collectively, we successfully constructed the ICD scheme, which enables predicting of the prognosis or immunotherapy efficacy in ESCC patients and uncovered the critical interplay between cell death signals and immune status in ESCC.
Introduction
Esophageal cancer (EC) is the 10th common cancer worldwide, while its mortality ranks sixth among all types of cancers (Sung et al., 2021). Esophageal squamous cell carcinoma (ESCC) accounts for nearly 90% of all the incident esophageal cancers each year, and regions of high incidence of ESCC include Eastern to Central Asia (Abnet et al., 2018). Even though great effort has been devoted to investigating its molecular characteristic (Cancer Genome Atlas Research, 2017) effectively, a therapeutic strategy against metastasis ESCCs was still limited until the development of immunotherapy. Recently, immunotherapy combined with chemotherapy has been proven effective in the first-line treatment of ESCC (Luo et al., 2021; Sun et al., 2021), which indicates that the therapeutic strategy of metastasis ESCC is stepping into a new era. Interestingly, different chemotherapy methods were chosen in these two clinical trials, among which Luo et al. (Luo et al., 2021) used Taxol plus platinum-based chemotherapy, while Sun et al. (Sun et al., 2021) used fluoropyrimidine plus platinum-based chemotherapy. As the efficacy of the former used in all patients is comparable to the latter used in high-PD-L1 patients, how to choose optimal chemotherapy for combinations based on the mechanism remains to be clarified.
As up to 60% of solid cancer patients do not respond to single-agent ICI therapy, the combination strategy breaks new ground for the use of immunotherapy (Gandhi et al., 2018; Paz-Ares et al., 2018; Galsky et al., 2020; Shitara et al., 2020). Despite huge success made in clinical trials, the mechanisms of how chemotherapy enhancing immunotherapy are to be fully understood yet, which would make choosing the optimal combination strategy a huge challenge.
Recently, scientists found that chemotherapy can induce immunogenic cell death (ICD) of tumor cells (Galluzzi et al., 2020), which could significantly enhance the anti-tumor immunity in a tumor microenvironment. However, different types of cell death have different effects on anti-tumor immunity (Legrand et al., 2019). For example, intrinsic apoptosis induced by mitochondrial outer membrane permeabilization allows the host to quickly and efficiently clear away dead cell corpses without triggering an immune response (Arandjelovic and Ravichandran, 2015). Autophagy could suppress antigen processing and presentation of tumor, and inhibition of autophagy may enhance immunotherapy (Deng et al., 2021; Thorburn and Towers, 2021). On the contrary, necroptosis provoked by TNF signaling of IFN stimulation is strongly correlated with activated immunity (Newton et al., 2016), which ensures that a powerful alert message is sent to the immune system. Similarly, pyroptosis is also an immunogenic cell death since its hyperinflammatory nature through IL-1 secretion and DAMP release could trigger inflammatory responses in the tumor microenvironment and increase leukocyte infiltration (Legrand et al., 2019). A recent research study also showed that tumor ferroptosis could be promoted by CD8+ T cell (Wang et al., 2019a), which proved that ferroptosis is another type of cell death that is positively associated with activated anti-tumor immunity. As the types and levels of cell death do have huge difference across different tumors, subtyping tumors with a cell death signal could help us characterize cancer patients concerning cell death and subsequently uncover the connection across cell death, immune microenvironment, and clinical features. Moreover, the subtyping would further help us choose specific agents that could induce specific cell death with high immunogenicity, which is an important step for making the most of an immunotherapy combination strategy.
Specific cell death signatures and cell death subtyping such as pyroptosis-related signature (Ye et al., 2021), pyroptosis-related subtyping (Song et al., 2021), and immunogenomic gene signature of cell-death-associated genes (Ahluwalia et al., 2021) have been proven that they could stratify cancer patients into different subgroups that are strongly associated with clinical features and immune microenvironment. However, none of this research combined different cell death signals together to characterize the tumors and systematically subtype the tumors concerning all immunogenic cell death in different omics levels. Here, with the transcriptome data of 857 patients and proteomic data of 124 patients, we found that signals of necroptosis, pyroptosis, and ferroptosis were positively associated with activated immunity in ESCC and subsequently defined ICD-high and ICD-low subtypes of the ESCCs with these three signals in both transcriptomic and proteomic levels. The subtyping could not only help us systematically understand the connection between cell death and immunity in ESCC but also help predicting the prognosis of ESCC patients and immunotherapy efficacy.
Methods
ESCC Data Source and Preprocessing
A total of seven ESCC microarray cohorts (GSE53625 (Li et al., 2014), GSE23400 (Su et al., 2011), GSE38129 (Hu et al., 2015), GSE47404 (Sawada et al., 2015), GSE69925 (Tanaka et al., 2015), GSE121931 (Li et al., 2020), and GSE44021 (Yang et al., 2020)) were collected from the GEO database (https//www.ncbi.nlm.nih.gov/geo/) (Table 1). A total of three GEO datasets (GSE53625, GSE23400, and GSE38129) were enrolled to perform a comparison between tumor and normal samples, with the standard that the number of tumor/normal sample pairs from the identical dataset and the identical platform was no less than 30. In total, 262 paired tumor and normal samples were merged. The tumor transcriptomics cohort consists of 775 tumor samples from the aforementioned GEO datasets, with the standard that the number of tumor samples from the sample dataset and the identical platform was no less than 50. Non-biological technical biases were corrected by “Combat” function of the sva R package (Leek et al., 2012). Mutation, copy number variation, and transcriptome data of ESCC from TCGA database were downloaded from the UCSC Xena datahub (https://xenabrowser.net/datapages/). The proteomic data and corresponding clinical data for 124 ESCC samples were obtained from Liu et al.‘s study(Liu et al., 2021). In total, four published immunotherapy datasets were collected to support the predictive value of ICD-score in immune checkpoint blockade treatment, including three SKCM cohorts and 1 urothelial carcinoma. (Van Allen et al., 2015; Mariathasan et al., 2018; Gide et al., 2019; Liu et al., 2019).
TABLE 1. Summary of all datasets in this study.
Immunogenic Cell-Death Pathway and Consensus Cluster
We summarized five cell death pathways which are apoptosis, autophagy, ferroptosis, necroptosis, and pyroptosis from the GO database (http://geneontology.org/), KEGG database (https://www.genome.jp/kegg/pathway.html), and Reactome database (https://reactome.org/). Gene set variation analysis (GSVA) was used to evaluate the relative activity of these cell death pathways (Hänzelmann et al., 2013). By comparing the correlation between immune-associated pathways and five cell death pathways, we selected the three most related pathways (ferroptosis, necroptosis, and pyroptosis), termed as the immunogenic cell death (ICD) pathway, to perform the clustering analysis. An unsupervised clustering analysis was applied to identify distinct cell death pathway activity patterns based on the GSVA score of the ICD pathway and classify patients for further analysis (Monti et al., 2003). The number of clusters and their stability was determined by the consensus clustering algorithm. The ConsensusCluserPlus package was utilized to perform the aforementioned steps, and 1,000 times repetitions were set to guarantee the robustness of clustering. The PCA algorithm was used to assess heterogeneity between different clusters.
Estimation of TME Cell Infiltration
We quantify the relative abundance of different immune cells in the ESCC microenvironment by GSVA. The marking gene set of each infiltration immune cell was obtained from Dvir Aran et al.‘s study (Aran et al., 2017). Quantification of tumor microenvironment signals including lymphocyte infiltration, macrophages, TGF-beta, IFN-gamma, and wound healing was performed on transcriptomic data of ESCC samples according to Thorsson et al. (Thorsson et al., 2018a).
Construction of the ICD-Score
On the basis of the most important genes in the ICD pathway, we constructed the ICD-score. Specifically, we first take an intersection of the ICD gene between transcriptomics data and protein data. Next, we sort the overlapping genes by the Extreme Gradient Boosting (XGBoost) algorithm according to the importance of distinguishing ESCC patients into ICD-high or ICD-low clusters. XGBoost is one of the most widely used tree-based boosting algorithms, in which a set of weak classifiers is combined to form a strong classifier sequentially (Chen and Guestrin, 20162016). In each iteration, misclassification errors of a previous classifier were corrected by the next classifier to perform more accurately. More importantly, it can avoid overfitting by effective ways. Importance is applied to access the contribution of each variable to the classification, which is measured for a single decision tree by the amount that each attribute split point improves the performance measure. The importance of each variable is the average of all decision trees. We selected top 30% genes ranked by importance to construct the ICD-score by calculating their arithmetic mean. We build the protein–protein interaction network of 14 genes with the STRING (https://www.string-db.org/) database and IMEx database (Breuer et al., 2013) in NetworkAnalyst (Zhou et al., 2019).
Statistical Analysis
The correlation coefficient between two numeric variables’ cell death was calculated using the Spearman coefficient. The Wilcoxon rank-sum test was used to compare the differences of numeric variables between two groups. Univariate Cox regression was applied to evaluate the relationship between different subgroups and prognosis. Multivariate Cox regression was used to identify independent prognostic factors. Multi-factor Cox regression results were visualized by the ‘forestplot’ R package. All statistics p values were two-tailed, and p < 0.05 was regarded as statistically significance. All data were processed using R 4.1.0.
Result
Identification of the Immunogenic Cell Death Pathway in ESCC
In total, five cell death pathways and the related genes were collected from a public database for the analysis, including apoptosis (234 genes), autophagy (169 genes), ferroptosis (40 genes), necroptosis (40 genes), and pyroptosis (50 genes) (Supplementary Table S1). The activity of apoptosis, ferroptosis, necroptosis, and pyroptosis in tumor tissue was significantly higher than that in normal tissue, while the difference of apoptosis was moderate. More interestingly, the autophagy pathway showed a reverse trend (Figure 1A). The result revealed that differences exist between normal tissue and tumor tissue concerning the types of cell death. In the transcriptomic level, a significantly positive correlation was observed between ferroptosis, necroptosis, pyroptosis, and immune-associated pathways, such as immune effector process and activation of immune response (Figure 1B). On the contrary, apoptosis and autophagy showed no significantly positive correlation with immune response signals. The result was further validated in the proteomics level (Figure 1C). Accordingly, only ferroptosis, necroptosis, and pyroptosis have a strong correlation with immune response signals, and we consequently defined the three cell death pathways as immune-associated cell death (ICD) pathways.
FIGURE 1. Identification of the immune cell death pathway. (A) GSVA score of cell death signaling between tumor and normal tissue in the transcriptomics level; (B,C) correlation of five cell death signaling and immune responding-related signaling in the transcriptomics level (B) and proteomics level (C).
The Landscape of Genetic Variation of Cell Death Pathways in ESCC
In order to investigate the genomic alteration in cell death pathways, we summarized the incidence of somatic mutation and copy number variation (CNV) concerning the corresponding genes in TCGA-ESCC cohort. Among 93 ESCC patients, genetic mutation of ICD pathways occurred in 86 patients, with a frequency of 92.47%, primarily including missense mutation, nonsense mutation, and splice site. Ranked by mutation frequency, TP53, a well-known tumor suppressor gene (Olivier et al., 2010; Fromentel and Soussi, 1992), is the most prevalently mutated genes among all the ICD genes. PIK3CA and PLEC were the following two genes, with frequencies of 13 and 5%, respectively (Figure 2A). As ESCC is a tumor with a high copy number variation (CNV) burden (Cui et al., 2020), CNV also occurs frequently on ICD genes. Specifically, TP63 in pyroptosis, TFRC in ferroptosis, and FADD in necroptosis presented the highest frequency of amplification, while CHMP2B, GPX4, and ELANE in pyroptosis, ACSL3 in ferropotosis, and BOK in necroptosis had the highest frequency of deletion (Figures 2B–D). In addition, CNV of different genes had different effects on the activity of corresponding cell death pathways. CHMP2B, TF, and ITPK1 were three genes whose copy numbers were most positively correlated with the respective pathways, while NLRP3, GCLM, and FASLG were three genes whose copy numbers were most positively correlated with the respective pathways (Supplementary Figure S1A). The result showed that genomic alteration of ICD genes was common in ESCC, which implied that the irregulation of ICD pathways may contribute to the tumorigenesis and development ESCC.
FIGURE 2. Landscape of genetic alterations of cell death genes in ESCC. (A) Mutation frequency of cell death genes in 91 patients with ESCC from the TCGA-ESCC cohort. Upper: the bar plot shows the number of mutant genes in individual patients. Right: the bar plot indicates the proportion of each variant type. (B–D). CNV frequency of pyroptosis, (B) ferroptosis, and (C) necroptosis (D) in the TCGA-ESCC cohort. The height of the column represents the alteration frequency. Variation color indicates the class of CNV (red for amplifications and blue for deletions).
Subtyping ESCCs With the Signal of ICD Pathways
A consensus clustering analysis based on ICD pathway GSVA scores was performed on the integrated dataset after batch effect removal (Supplementary Figures S2A,B). After selecting the optimal k value (Figure 3A, Supplementary Figure S3A), two distinct ICD patterns were eventually identified in the transcription level (Figure 3A). PCA analysis also revealed the heterogeneity of the two groups concerning cell death signals (Figure 3B). An identical phenomenon was also observed in the proteomics dataset (Supplementary Figures S4A,B). By comparing the relative level of the cell death signal, we found that ferroptosis, necroptosis, and pyroptosis signals were extremely variant between the two groups (Figure 3C, Supplementary Figure 4C). Thus, we defined the group with higher ICD pathway signals as an ICD-high cluster, and the other group as an ICD-low cluster.
FIGURE 3. Patterns of the immune cell death pathway and prognosis characteristics of each pattern. (A) Unsupervised clustering of three immune cell death pathways’ GSVA scores in the transcriptomics level. (B) PCA analysis for two ICD clusters in the transcriptomics level. (C) Cell death pathway patterns of two ICD clusters in the transcriptomics level. Top: the bar plot shows the cell death pathway GSVA score of two ICD clusters. Bottom: the heatmap reveals the GSVA score pattern of two ICD clusters (D–F). Kaplan–Meier analysis of overall survival (OS) or progression-free survival (PFS) comparing the ESCC patients with two ICD clusters. (D) Kaplan–Meier analysis of OS in the transcriptomics level. (E) Kaplan–Meier analysis of PFS in the transcriptomics level. (F) Kaplan–Meier analysis of OS in the proteomics level (red for the ICD-high cluster; blue for the ICD-low cluster).
By integrating the clinical data, we found that ESCC patients in the ICD-high cluster had a better overall survival than those in the ICD-low cluster (Figure 3D, HR = 0.62, p-value = 0.01) in the transcriptomic cohort. Similarly, longer progression-free survival (Figure 3E, HR = 0.62, p-value = 0.03) and overall survival (Figure 3F, HR = 0.62, p-value = 0.08) were also detected in ICD-high patients than in ICD-low patients in the proteomic level, which further validated the prognostic relevance of the ICD subtyping.
The Immune Characteristic of Two ICD-Associated Patterns
As better prognosis was detected in ICD-high patients, we speculated that the strong correlation with anti-tumor immunity might be one of the explanations for the clinical relevance. To further clarify the tumor immune microenvironment between the two groups in detail, we assessed the level of infiltrated immune cell, immune subtyping signals, and immunomodulator genes in these two groups. As expected, the abundance of most infiltrated immune cells was significantly higher in the ICD-high cluster than in the ICD-low cluster (Figure 4A), which was validated with the proteomics dataset (Supplementary Figure S5A). Specifically, we noticed that myeloid cells, such as dendritic cells, neutrophils, monocytes, and macrophages were the most upregulated cell types in the ICD-high cluster (Figures 4B,C). It is well known that dendritic cells and macrophages are the main antigen-presenting cells among the immune reaction. As ferroptosis, necroptosis, and pyroptosis were reportedly associated with immune activation in the tumor microenvironment, we assumed that the higher level of neoantigen presentation activity might have occurred in the ICD-high cluster, probably due to more release of the immunogenic tumor neoantigen. Vesteinn Thorsson et al. (Thorsson et al., 2018b) categorized cancer immunity into six subtypes by five signatures (lymphocyte infiltration, macrophages, TGF-beta, IFN-gamma, and wound healing). As expected, signals directly concerning the immune reaction were significantly higher in the ICD-high cluster, especially macrophage signals (Figures 4D,E). In addition, we also assessed the expression level of immunomodulator genes between the two groups. Both immunomodulatory ligands and receptors tend to have a higher expression level in the ICD-high group than in the ICD-low group (Figures 4F,G). In addition, the expressions of markers of cytotoxic activity were also prevailing in the ICD-high cluster (Figure 4H), suggesting a more activated anti-tumor immunity in ICD-high patients. Conclusively, these results implied that the ICD-high cluster has a relatively hotter tumor immune environment.
FIGURE 4. Immune infiltration characteristics of two ICD clusters. (A) Heatmap presents the filtration of immune cells in two ICD clusters (B,C). The enrichment of major antigen-presenting cells of two clusters in the transcriptomics level (B) and proteomics level (C–E). GSVA scores of five immune signatures of two ICD clusters in the transcriptomics level (D) and proteomics level (E,F). The mean expression of immunoregulator ligands of two ICD clusters in the transcriptomics level. (G) Mean expression of immunoregulator receptors of two ICD clusters in the transcriptomics level. (H) Mean expression of cytotoxic activity signature of two ICD-score clusters in the transcriptomics level. (Mann–Whitney U test: NS: no statistical difference, *: p < 0.05, **: p < 0.01, ***: p < 0.0001) (red for the ICD-score high cluster; blue for the ICD-score low cluster).
Construction of ICD-Score and Investigation of Its Biological Relevance
Our findings have already shown the potential role of the ICD pathway in prognosis and immune infiltration modulation. To further push forward the use of the subtyping of the ICD signal, denoting the subtyping with a brief gene signature would be more helpful. Thus, on the basis of the ICD genes, we constructed a set of scoring systems to quantify the ICD pathway pattern of individual ESCC patients, termed as the ICD-score. In order to select optimal features to construct a more concise model, we used the XGBoost algorithm, which has been proved to be better than other machine learning algorithms on classification and feature selection in several research studies (Huang et al., 2020; Le et al., 2020; Chen et al., 2021). Accordingly, we obtained 14 genes to profile the ICD-score (Figure 5A), among which PRNP was the most important gene, followed by MAP3K5 and SCAF11. A protein–protein intersection network built by different databases revealed the tight connection between the 14 proteins (Figure 5B, Supplementary Figure S6A). Similar to subtyping, the ICD-score was also positively correlated with infiltration of myeloid cells in charge of neoantigen presentation, such as dendritic cells, neutrophils, monocytes, and macrophages (Supplementary Figure S6B). As expected, the ICD-score low cluster had a worse overall survival rate (Figures 5C,D). Adjusting the influence of the TNM stage by multi-variable Cox regression, we found that the ICD-score can act as an independent factor to predict the prognosis of ESCC patients (transcriptomics: log-rank p = 0.02, HR = 1.58 [1.07-2.34]; proteomics: log-rank p = 0.02, HR = 2.08[1.12-3.85]) (Figure 5E, Supplementary Figure S6C). Our analysis suggested that the ICD-score model can represent ICD patterns and act as an independent factor in predicting the prognosis of ESCC patients.
FIGURE 5. Construction of the ICD score and prognosis characteristics. The importance of 14 genes calculated by the XGBoost algorithm. (A) Higher importance value of a gene indicated higher importance. (B) Protein–protein interaction network of 14 proteins built by the STRING database. The nodes represent the proteins; the line represents the relationship of proteins. (C,D) Kaplan–Meier analysis of overall survival (OS) comparing the ESCC patients with the two ICD score clusters. (C) Kaplan–Meier analysis of OS in the transcriptomics level. (D) Kaplan–Meier analysis of OS in the proteomics levels (red for the ICD-high cluster; blue for the ICD-low cluster). (E) Forest plot: Cox multivariate mode of the ICD-score and the TNM stage in the transcriptomics level.
Correlation Between Oncogenic Alteration and ICD-Score
The correlation between ICD-score and hallmark pathways was systematically estimated with the transcriptomic cohort. Apoptosis, protein secretion, TP53, complement, and TNF-alpha signaling via NF-kappaB had a strong positive association with the ICD-score, while KRAS signaling, WNT beta catenin signaling, hedgehog signaling, and spermatogenesis were negatively associated with the ICD-score (Figure 6A). Among the 10 oncogenic signaling pathways reported previously (Sanchez-Vega et al., 2018), we found that the TP53 pathway alteration and RTK-RAS alteration were both significantly positively associated with the ICD-score (Figure 6B). Chromosomes 11q13.3 and 9q21.3 were the most two common chromosome alterations in the development of ESCC (Chang et al., 2017; Zhang et al., 2018). In our study, the chromosome 11q13 amplification group exhibited higher ICD-scores than the normal group (p = 0.012), while 9q21 deletion showed no statistical difference with the amplicon and normal groups (Figures 6C,D). The result enables us to connect the oncogenic genomic alteration and the ICD-score in ESCC, which would help us to further understand the possible regulated upstream of the ICD subtypes.
FIGURE 6. Association between the ICD-score and hallmark pathways. (A) Correlation of the ICD score and hallmark pathways (red for positive correlation; gray for negative correlation). (B) ICD score of mutated and wild-type groups of hallmark pathways in the TCGA-ESCC cohort. (C,D). ICD score comparison of patients with common copy number alteration in ESCC and normal patients. (C) 11q13.3 amplification group vs. normal group. (D) 9p21.3 deletion group vs. non-deletion group.
The Predictive Value of ICD-Scores in Immune-Oncology Therapy
ICI treatment represented by PD-1/CTLA4 inhibitors has undoubtedly caused a major breakthrough in oncology therapy. Considering the strong correlation between ICD subtyping and tumor immune microenvironment, we further investigated whether ICD-score could predict patients’ response to immune checkpoint therapy based on four immunotherapy cohorts. As expected, patients with a higher ICD-score derived significantly longer OS than those with a lower ICD-score (Figures 7A–D), revealing that the ICD-score might be another promising biomarker in immunotherapy.
FIGURE 7. Survival analysis of ICD-score high/low clusters in four immunotherapy cohorts. (A–D) Kaplan–Meier analysis of overall survival in Tuba N Gide et al.‘s cohort (A), David Liu et al.‘s cohort (B), Eliezer M Van Allen et al.‘s cohort (C), and Sanjeev Mariathasan et al.‘s cohort (D) (red for the ICD-score high cluster; blue for the ICD-score low cluster).
Discussion
Accurate molecular typing is an important part of precision medicine in the future. In the past few years, with the development of high-throughput sequencing technology, molecular typing of esophageal cancer based on transcriptome expression data has made great progress. For example, Fengjing Wang et al. obtained two ESCC groups with different clinical survival outcomes from the transcriptome data of 179 ESCC cases and identified the corresponding biomarkers such as FOXA1 and EYA2 for subtype I and LAMC2 and KRT14 for subtype II (Wang et al., 2019b). Meng Liu et al. applied consensus clustering methods on transcriptome data from 125 Xinjiang patients with ESCC to divide the patients into three types and verified the subtypes in two independent cohorts including the TCGA dataset (Liu et al., 2020). It is worth noting that the aforementioned studies were mainly based on the analysis from the whole transcriptome, with limited biological insights, and the discovery dataset was from a single dataset. The main features of our study are that our clustering was inspired by the heterogeneity of ICD in ESCC patients, and the subtypes were obtained by an unsupervised cluster analysis of the patients’ ICD GSVA score from six studies, and these subtypes were further validated and confirmed in an extraproteomic dataset. To our knowledge, the cohort of our study is so far the largest dataset for building ESCC molecular subtypes. Moreover, to facilitate the usage of the ICD scheme in the future, the XGBoost was applied to shrink the gene list representing each ICD group, making a simplified subtyping system. Therefore, our study provides an important immune-associated cell death signal molecular subtyping scheme for ESCC classification in the era of precision medicine.
Regulated cell death (RCD) is an important mechanism for organisms to protect themselves from various diseases and even cancers (Galluzzi et al., 2018). Previous studies have revealed several types of RCDs which were proved to cope with different kinds of adverse conditions. The most common RCDs include intrinsic and extrinsic apoptosis, necroptosis, proptosis, and ferroptosis, which have also been shown to be related to the immunogenic response of tumors (Legrand et al., 2019). In the cancer basic research field, there have been many studies that have explored abnormal different cell death signals in tumors and their crucial roles in tumor immunity. However, in ESCC, the systematic study of interactions between these RCDs and the immune microenvironment has not been carried out yet. In the present study, by performing comprehensive bioinformatics analysis of six published transcriptomic datasets, one proteomic dataset and four ICC treatment cohorts, we characterized all known cell death signals status in ESCCs. With machine learning approaches, we constructed an applicable ICD scoring scheme to classify the ESCC into subtypes, which could be further characterized by distinct genomic features as well as immune profiles. We expected this scheme could be utilized not only in research on the molecular mechanisms governing ESCC but also as a potential prognostic and immunotherapy biomarker for this deadly cancer.
Data Availability Statement
Publicly available datasets were analyzed in this study. These data can be found at: we have submitted the source data supporting our article to www.jianguoyun.com and shared the download link in Download link of Integrated Data.docx as a Supplementary Material. Information of all published datasets used in the research study is listed in Table 1 in the article. The intermediate files and code are deposited in https://github.com/chenyx47/ESCC_ICD.git/,and the code is also attached with the filename: Code.zip. All the material is available for supporting our article.
Author Contributions
YZ and YC contributed to the conception and design of the study. YC organized the database and performed the statistical analysis. YZ and YC wrote the first draft of the manuscript. All authors contributed to manuscript revision and read and approved the submitted version.
Funding Statements
The research study was funded by the Innovation Project of Shandong Academy of Medical Sciences (2019-04), the Academic Promotion Program of Shandong First Medical University (2019ZL002), the National Key Research and Development Program of China (2018YFE0114100), the Key Research and Development Program of Shandong province, China (2019GGX101057), and the Science Technology Program of Jinan (201805051).
Conflict of Interest
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.
Publisher’s Note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2022.855404/full#supplementary-material
Supplementary Table S1 | Detail of the genes of the selected cell death pathways and their sources.
References
Abnet, C. C., Arnold, M., and Wei, W.-Q. (2018). Epidemiology of Esophageal Squamous Cell Carcinoma. Gastroenterology 154 (2), 360–373. doi:10.1053/j.gastro.2017.08.023
Ahluwalia, P., Ahluwalia, M., Mondal, A. K., Sahajpal, N., Kota, V., Rojiani, M. V., et al. (2021). Immunogenomic Gene Signature of Cell-Death Associated Genes with Prognostic Implications in Lung Cancer. Cancers (Basel) 13 (1). doi:10.3390/cancers13010155
Aran, D., Hu, Z., and Butte, A. J. (2017). xCell: Digitally Portraying the Tissue Cellular Heterogeneity Landscape. Genome Biol. 18 (1), 220. doi:10.1186/s13059-017-1349-1
Arandjelovic, S., and Ravichandran, K. S. (2015). Phagocytosis of Apoptotic Cells in Homeostasis. Nat. Immunol. 16 (9), 907–917. doi:10.1038/ni.3253
Breuer, K., Foroushani, A. K., Laird, M. R., Chen, C., Sribnaia, A., Lo, R., et al. (2013). InnateDB: Systems Biology of Innate Immunity and Beyond-Rrecent Updates and Continuing Curation. Nucleic Acids Res. 41 (Database issue), D1228–D1233. doi:10.1093/nar/gks1147
Cancer Genome Atlas Research, N. (2017). Integrated Genomic Characterization of Oesophageal Carcinoma. Nature 541 (7636), 169–175. doi:10.1038/nature20805
Chang, J., Tan, W., Ling, Z., Xi, R., Shao, M., Chen, M., et al. (2017). Genomic Analysis of Oesophageal Squamous-Cell Carcinoma Identifies Alcohol Drinking-Related Mutation Signature and Genomic Alterations. Nat. Commun. 8, 15290. doi:10.1038/ncomms15290
Chen, C., Yang, D., Gao, S., Zhang, Y., Chen, L., Wang, B., et al. (2021). Development and Performance Assessment of Novel Machine Learning Models to Predict Pneumonia after Liver Transplantation. Respir. Res. 22 (1), 94. doi:10.1186/s12931-021-01690-3
Chen, T., and Guestrin, C. (20162016). “Xgboost: A Scalable Tree Boosting System,” in Proceedings of the 22nd Acm Sigkdd International Conference on Knowledge Discovery and Data Mining, 785–794.
Cui, Y., Chen, H., Xi, R., Cui, H., Zhao, Y., Xu, E., et al. (2020). Whole-genome Sequencing of 508 Patients Identifies Key Molecular Features Associated with Poor Prognosis in Esophageal Squamous Cell Carcinoma. Cell Res 30 (10), 902–913. doi:10.1038/s41422-020-0333-6
Deng, J., Thennavan, A., Dolgalev, I., Chen, T., Li, J., Marzio, A., et al. (2021). ULK1 Inhibition Overcomes Compromised Antigen Presentation and Restores Antitumor Immunity in LKB1-Mutant Lung Cancer. Nat. Cancer 2 (5), 503–514. doi:10.1038/s43018-021-00208-6
Fromentel, C. C. D., and Soussi, T. (1992). TP53 Tumor Suppressor Gene: a Model for Investigating Human Mutagenesis. Genes Chromosom. Cancer 4 (1), 1–15. doi:10.1002/gcc.2870040102
Galluzzi, L., Vitale, I., Aaronson, S. A., Abrams, J. M., Adam, D., Agostinis, P., et al. (2018). Molecular Mechanisms of Cell Death: Recommendations of the Nomenclature Committee on Cell Death 2018. Cell Death Differ 25 (3), 486–541. doi:10.1038/s41418-017-0012-4
Galluzzi, L., Humeau, J., Buqué, A., Zitvogel, L., and Kroemer, G. (2020). Immunostimulation with Chemotherapy in the Era of Immune Checkpoint Inhibitors. Nat. Rev. Clin. Oncol. 17 (12), 725–741. doi:10.1038/s41571-020-0413-z
Galsky, M. D., Arija, J. Á. A., Bamias, A., Davis, I. D., De Santis, M., Kikuchi, E., et al. (2020). Atezolizumab with or without Chemotherapy in Metastatic Urothelial Cancer (IMvigor130): a Multicentre, Randomised, Placebo-Controlled Phase 3 Trial. The Lancet 395 (10236), 1547–1557. doi:10.1016/s0140-6736(20)30230-0
Gandhi, L., Rodríguez-Abreu, D., Gadgeel, S., Esteban, E., Felip, E., De Angelis, F., et al. (2018). Pembrolizumab Plus Chemotherapy in Metastatic Non-small-cell Lung Cancer. N. Engl. J. Med. 378 (22), 2078–2092. doi:10.1056/nejmoa1801005
Gide, T. N., Quek, C., Menzies, A. M., Tasker, A. T., Shang, P., Holst, J., et al. (2019). Distinct Immune Cell Populations Define Response to Anti-PD-1 Monotherapy and Anti-PD-1/Anti-CTLA-4 Combined Therapy. Cancer Cell 35 (2), 238–255. doi:10.1016/j.ccell.2019.01.003
Hänzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: Gene Set Variation Analysis for Microarray and RNA-Seq Data. BMC Bioinformatics 14, 7. doi:10.1186/1471-2105-14-7
Hu, N., Wang, C., Clifford, R. J., Yang, H. H., Su, H., Wang, L., et al. (2015). Integrative Genomics Analysis of Genes with Biallelic Loss and its Relation to the Expression of mRNA and Micro-RNA in Esophageal Squamous Cell Carcinoma. BMC Genomics 16 (1), 732. doi:10.1186/s12864-015-1919-0
Huang, Z., Hu, C., Chi, C., Jiang, Z., Tong, Y., and Zhao, C. (2020). An Artificial Intelligence Model for Predicting 1-Year Survival of Bone Metastases in Non-small-cell Lung Cancer Patients Based on XGBoost Algorithm. Biomed. Res. Int. 2020, 3462363. doi:10.1155/2020/3462363
Le, N. Q. K., Do, D. T., Chiu, F. Y., Yapp, E. K. Y., Yeh, H. Y., and Chen, C. Y. (2020). XGBoost Improves Classification of MGMT Promoter Methylation Status in IDH1 Wildtype Glioblastoma. J. Pers Med. 10 (3). doi:10.3390/jpm10030128
Leek, J. T., Johnson, W. E., Parker, H. S., Jaffe, A. E., and Storey, J. D. (2012). The Sva Package for Removing Batch Effects and Other Unwanted Variation in High-Throughput Experiments. Bioinformatics 28 (6), 882–883. doi:10.1093/bioinformatics/bts034
Legrand, A. J., Konstantinou, M., Goode, E. F., and Meier, P. (2019). The Diversification of Cell Death and Immunity: Memento Mori. Mol. Cel 76 (2), 232–242. doi:10.1016/j.molcel.2019.09.006
Li, J., Chen, Z., Tian, L., Zhou, C., He, M. Y., Gao, Y., et al. (2014). LncRNA Profile Study Reveals a Three-lncRNA Signature Associated with the Survival of Patients with Oesophageal Squamous Cell Carcinoma. Gut 63 (11), 1700–1710. doi:10.1136/gutjnl-2013-305806
Li, M., Zhao, J., Li, X., Chen, Y., Feng, C., Qian, F., et al. (2020). HiFreSP: A Novel High-Frequency Sub-pathway Mining Approach to Identify Robust Prognostic Gene Signatures. Brief Bioinform 21 (4), 1411–1424. doi:10.1093/bib/bbz078
Liu, D., Schilling, B., Liu, D., Sucker, A., Livingstone, E., Jerby-Arnon, L., et al. (2019). Integrative Molecular and Clinical Modeling of Clinical Outcomes to PD1 Blockade in Patients with Metastatic Melanoma. Nat. Med. 25 (12), 1916–1927. doi:10.1038/s41591-019-0654-5
Liu, M., An, H., Zhang, Y., Sun, W., Cheng, S., Wang, R., et al. (2020). Molecular Analysis of Chinese Oesophageal Squamous Cell Carcinoma Identifies Novel Subtypes Associated with Distinct Clinical Outcomes. EBioMedicine 57, 102831. doi:10.1016/j.ebiom.2020.102831
Liu, W., Xie, L., He, Y.-H., Wu, Z.-Y., Liu, L.-X., Bai, X.-F., et al. (2021). Large-scale and High-Resolution Mass Spectrometry-Based Proteomics Profiling Defines Molecular Subtypes of Esophageal Cancer for Therapeutic Targeting. Nat. Commun. 12 (1), 4961. doi:10.1038/s41467-021-25202-5
Luo, H., Lu, J., Bai, Y., Mao, T., Wang, J., Fan, Q., et al. (2021). Effect of Camrelizumab vs Placebo Added to Chemotherapy on Survival and Progression-free Survival in Patients with Advanced or Metastatic Esophageal Squamous Cell Carcinoma. JAMA 326 (10), 916–925. doi:10.1001/jama.2021.12836
Mariathasan, S., Turley, S. J., Nickles, D., Castiglioni, A., Yuen, K., Wang, Y., et al. (2018). TGFβ Attenuates Tumour Response to PD-L1 Blockade by Contributing to Exclusion of T Cells. Nature 554 (7693), 544–548. doi:10.1038/nature25501
Monti, S., Tamayo, P., Mesirov, J., and Golub, T. (2003). Consensus Clustering: a Resampling-Based Method for Class Discovery and Visualization of Gene Expression Microarray Data. Machine Learn. 52 (1), 91–118. doi:10.1023/a:1023949509487
Newton, K., Wickliffe, K. E., Maltzman, A., Dugger, D. L., Strasser, A., Pham, V. C., et al. (2016). RIPK1 Inhibits ZBP1-Driven Necroptosis during Development. Nature 540 (7631), 129–133. doi:10.1038/nature20559
Olivier, M., Hollstein, M., and Hainaut, P. (2010). TP53 Mutations in Human Cancers: Origins, Consequences, and Clinical Use. Cold Spring Harbor Perspect. Biol. 2 (1), a001008. doi:10.1101/cshperspect.a001008
Paz-Ares, L., Luft, A., Vicente, D., Tafreshi, A., Gümüş, M., Mazières, J., et al. (2018). Pembrolizumab Plus Chemotherapy for Squamous Non-small-cell Lung Cancer. N. Engl. J. Med. 379 (21), 2040–2051. doi:10.1056/nejmoa1810865
Sanchez-Vega, F., Mina, M., Armenia, J., Chatila, W. K., Luna, A., La, K. C., et al. (2018). Oncogenic Signaling Pathways in the Cancer Genome Atlas. Cell 173 (2), 321–e10. doi:10.1016/j.cell.2018.03.035
Sawada, G., Niida, A., Hirata, H., Komatsu, H., Uchi, R., Shimamura, T., et al. (2015). An Integrative Analysis to Identify Driver Genes in Esophageal Squamous Cell Carcinoma. PLoS One 10 (10), e0139808. doi:10.1371/journal.pone.0139808
Shitara, K., Van Cutsem, E., Bang, Y.-J., Fuchs, C., Wyrwicz, L., Lee, K.-W., et al. (2020). Efficacy and Safety of Pembrolizumab or Pembrolizumab Plus Chemotherapy vs Chemotherapy Alone for Patients with First-Line, Advanced Gastric Cancer. JAMA Oncol. 6 (10), 1571–1580. doi:10.1001/jamaoncol.2020.3370
Song, W., Ren, J., Xiang, R., Kong, C., and Fu, T. (2021). Identification of Pyroptosis-Related Subtypes, the Development of a Prognosis Model, and Characterization of Tumor Microenvironment Infiltration in Colorectal Cancer. Oncoimmunology 10 (1), 1987636. doi:10.1080/2162402x.2021.1987636
Su, H., Hu, N., Yang, H. H., Wang, C., Takikita, M., Wang, Q.-H., et al. (2011). Global Gene Expression Profiling and Validation in Esophageal Squamous Cell Carcinoma and its Association with Clinical Phenotypes. Clin. Cancer Res. 17 (9), 2955–2966. doi:10.1158/1078-0432.ccr-10-2724
Sun, J.-M., Shen, L., Shah, M. A., Enzinger, P., Adenis, A., Doi, T., et al. (2021). Pembrolizumab Plus Chemotherapy versus Chemotherapy Alone for First-Line Treatment of Advanced Oesophageal Cancer (KEYNOTE-590): a Randomised, Placebo-Controlled, Phase 3 Study. The Lancet 398 (10302), 759–771. doi:10.1016/s0140-6736(21)01234-4
Sung, H., Ferlay, J., Siegel, R. L., Laversanne, M., Soerjomataram, I., Jemal, A., et al. (2021). Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA A. Cancer J. Clin. 71 (3), 209–249. doi:10.3322/caac.21660
Tanaka, Y., Aoyagi, K., Minashi, K., Komatsuzaki, R., Komatsu, M., Chiwaki, F., et al. (2015). Discovery of a Good Responder Subtype of Esophageal Squamous Cell Carcinoma with Cytotoxic T-Lymphocyte Signatures Activated by Chemoradiotherapy. PLoS One 10 (12), e0143804. doi:10.1371/journal.pone.0143804
Thorburn, A., and Towers, C. G. (2021). Enhancing Anti-tumor Immunity by Autophagy Inhibition. Nat. Cancer 2 (5), 484–486. doi:10.1038/s43018-021-00214-8
Thorsson, V., Gibbs, D. L., Brown, S. D., Wolf, D., Bortone, D. S., Ou Yang, T. H., et al. (2018). The Immune Landscape of Cancer. Immunity 48 (4), 812–e14. doi:10.1016/j.immuni.2018.03.023
Thorsson, V., Gibbs, D. L., Brown, S. D., Wolf, D., Bortone, D. S., Ou Yang, T.-H., et al. (2018). The Immune Landscape of Cancer. Immunity 48 (4), 812–830. doi:10.1016/j.immuni.2018.03.023
Van Allen, E. M., Miao, D., Schilling, B., Shukla, S. A., Blank, C., Zimmer, L., et al. (2015). Genomic Correlates of Response to CTLA-4 Blockade in Metastatic Melanoma. Science 350 (6257), 207–211. doi:10.1126/science.aad0095
Wang, F., Yan, Z., Lv, J., Xin, J., Dang, Y., Sun, X., et al. (2019). Gene Expression Profiling Reveals Distinct Molecular Subtypes of Esophageal Squamous Cell Carcinoma in Asian Populations. Neoplasia 21 (6), 571–581. doi:10.1016/j.neo.2019.03.013
Wang, W., Green, M., Choi, J. E., Gijón, M., Kennedy, P. D., Johnson, J. K., et al. (2019). CD8+ T Cells Regulate Tumour Ferroptosis during Cancer Immunotherapy. Nature 569 (7755), 270–274. doi:10.1038/s41586-019-1170-y
Yang, H., Su, H., Hu, N., Wang, C., Wang, L., Giffen, C., et al. (2020). Integrated Analysis of Genome-wide miRNAs and Targeted Gene Expression in Esophageal Squamous Cell Carcinoma (ESCC) and Relation to Prognosis. BMC Cancer 20 (1), 388. doi:10.1186/s12885-020-06901-6
Ye, Y., Dai, Q., and Qi, H. (2021). A Novel Defined Pyroptosis-Related Gene Signature for Predicting the Prognosis of Ovarian Cancer. Cell Death Discov. 7 (1), 71. doi:10.1038/s41420-021-00451-x
Zhang, W., Hong, R., Li, L., Wang, Y., Du, P., Ou, Y., et al. (2018). The Chromosome 11q13.3 Amplification Associated Lymph Node Metastasis Is Driven by miR-548k through Modulating Tumor Microenvironment. Mol. Cancer 17 (1), 125. doi:10.1186/s12943-018-0871-4
Keywords: cell-death signals, esophageal squamous cell carcinoma, heterogeneity, prognosis, immune microenvironment
Citation: Zhang Y and Chen Y (2022) Stratification From Heterogeneity of the Cell-Death Signal Enables Prognosis Prediction and Immune Microenvironment Characterization in Esophageal Squamous Cell Carcinoma. Front. Cell Dev. Biol. 10:855404. doi: 10.3389/fcell.2022.855404
Received: 15 January 2022; Accepted: 14 March 2022;
Published: 12 April 2022.
Edited by:
Yong Yu, Johannes Kepler University of Linz, AustriaReviewed by:
Huiyun Liu, Tufts University School of Medicine, United StatesYueyuan Zheng, Cedars Sinai Medical Center, United States
Di Peng, Huazhong University of Science and Technology, China
Copyright © 2022 Zhang and Chen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Yiyuan Zhang, zhangyydahu@163.com
†These authors have contributed equally to this work