The prognostic value of an immune-related gene signature and infiltrating tumor immune cells based on bioinformatics analysis in primary esophageal cancer.

Background: Immunotherapy is a promising novel treatment for esophageal cancer (ESCA). However, previous studies provide limited direct information about the prognostic significance of immune-related genes (IRGs) in primary ESCA development. This study explored the prognostic value of IRGs and infiltrating tumor immune cells in primary ESCA.
Methods: The ribonucleic acid (RNA)-sequencing data and clinical information of primary ESCA were downloaded from The Cancer Genome Atlas (TCGA) database. Which included the clinical factors and prognosis outcomes of the ESCA patients. The IRGs were downloaded from the ImmPORT database.
Results: We established the robust IRG prognostic signature of 4 IRGs (i.e., heat shock protein family A member 6, Oncostatin M, androgen receptor, and nuclear receptor subfamily 2 group F member 2) in primary ESCA, and divided the ESCA patients into high- and low-risk groups based on overall survival (OS). The Kaplan-Meier curves showed the high predictive ability of the prognostic signature in the training, testing, and full data sets (P=2.407e-03, P=1.044e-02, and P=2.535e-04, respectively). Multivariate Cox regression analyses were performed with age, grade, tumor stage, tumor type and the risk score as covariables. The risk score supports the use of a prognostic signature as an independent prognostic factor [training data set: hazard ratio (HR) =1.185, 95% confidence interval (95% CI): 1.013-1.388, P=0.034; testing data set: HR =2.056, 95% CI: 1.015-4.166, P=0.045; full data set: HR =1.197, 95% CI: 1.059-1.354, P=0.004]. The area under the curve (AUC) of the receiver operating characteristic curve validated the high predictive accuracy of the IRG signature in the training, testing, and full data set (AUCs =0.808, 0.657, and 0.751, respectively). The infiltration level of the activated mast cells was significantly higher in the high-risk group than the low-risk group; thus, infiltrating mast cells are associated with worse OS in ESCA patients. Conclusions: Our IRG prognostic signature provides a new direction to predict the survival of primary ESCA patients and has the potential ability to establish, promote, and improve personalized treatment procedures based on each patient's risk. 2022 Journal of Gastrointestinal Oncology. All rights reserved.

Introduction

Esophageal cancer (ESCA) is the 8th most common cancer type and the 6th leading cause of cancer death worldwide (1). Most ESCA patients are diagnosed in the advanced stages (2). Despite improvements in diagnostic methods and treatments, the overall prognosis for ESCA patients remains poor, with a 5-year overall survival (OS) rate of approximately 20% (3). The accurate prognosis and risk assessment of patients will enable the personalized treatment and improve the survival rate of patients (4). Thus, promoting early diagnosis and discovering novel prognostic markers and therapeutic targets are essential in benefiting patients with ESCA.

Due to the tumor heterogeneity in ESCA, a single biomarker has limited value in prognosis. Gene expression profiles play a crucial role in the analysis of cancer development and prognosis. Studies have indicated that gene signatures, which consist of several hub genes, might be a good choice for predicting cancer prognosis (5-8). Recently, circular ribonucleic acid (RNA), long non-coding RNA, and microRNA gene signatures have been identified in ESCA patients (9-11). However, these findings lack of specificity, there is still a long way to put these in practice. It is necessary to identify the most effective and robust biomarkers and optimize treatment strategies for ESCA patients.

There is increasing evidence that the immune system plays a crucial role in cancer development and progression (12,13). Cancer immunotherapy is regarded as a major milestone and represents a paradigm shift in the treatment of cancer. Over the past decades, immunotherapy has become a powerful clinical strategy for treating cancer, including lung cancer (14), bladder cancer (15), and skin cancer (16). However, the use of immunotherapy to treat ESCA has led to mixed results, which is partially due to a lack of reliable predictive markers of treatment response (17). Therefore, it is urgent to identify a significant risk score model based on immune-related genes (IRGs) to optimize immunotherapy treatment strategies for ESCA patients. We present the following article in accordance with the TRIPOD reporting checklist (available at https://jgo.amegroups.com/article/view/10.21037/jgo-22-576/rc).

Methods

The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013). The present study sought examine the potential of IRGs as biomarkers to predict the prognosis and immunotherapy outcomes of ESCA patients. Our findings provide a foundation for the subsequent research of IRGs and the use of personalized strategies in the treatment of ESCA. By integrating the expression profiles of IRGs with clinical information, we developed an IRG signature in ESCA and validated its predictive accuracy. We also examined the infiltrating immune cells in the high- and low-risk patient groups.

Data extraction from TCGA database

The transcriptome RNA-sequencing (RNA-seq) data and clinical information for ESCA patients were acquired from The Cancer Genome Atlas (TCGA; https://portal.gdc.cancer.gov/). RNA-seq was performed using the Fragments per Kilobase of transcript per Million mapped reads (HTSeq-FPKM version: July 19, 2019). The current study included 158 primary ESCA samples, including squamous cell carcinoma and adenocarcinoma samples, and 10 normal samples. A list of 2498 IRGs was acquired from the Immunology Database and Analysis Portal (ImmPort; https://www.immport.org/home) (18).

Identification of differentially expressed IRGs and the functional analysis

The “Limma” package in R software was used to calculate the differentially expressed genes (DEGs) between the tumor and non-tumor samples from patients with ESCA (|logFold change (FC)| >1; P value <0.05). The “org.Hs.eg.db” package was used to match the “Entrez ID” of each IRG. Differentially expressed IRGs were IRGs included in the DEGs. Functional enrichment analyses were conducted to explore the potential molecular mechanisms of the differentially expressed IRGs via the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis pathways. Functional categories with an adjusted P value <0.05 were considered significant pathways.

TF and IRG regulatory network

The clinical data of TCGA-ESCA cohort was acquired and diagnoses of “0.days_to_death” for the dead patients and “0.days_to_last_follow_up” for the alive patients (as defined in the clinical data set) were defined as OS. The univariate regression analysis was performed using the “survival” package of R to select the prognostic IRGs and transcription factors (TFs) with a P value <0.01. The correlation between the prognosis-related IRGs and TFs was analyzed using the “cor.test” function of R. The cutoff criteria were set as a correlation coefficient >0.3 and a P value <0.05. Cytoscape 3.7.2 was used to construct and visualize the regulatory network.

The establishment and validation of the IRG prognostic signature

The ESCA patients in the TCGA-ESCA cohort were randomly divided into a training and a testing data sets using the “caret” package of R. A univariate Cox regression analysis was performed to selected IRGs correlated to OS, following a multivariate Cox regression analysis to identify the IRGs in the training data set that independently predicted OS. A prognostic signature was constructed according to the expression of the IRGs. A risk score was calculated using a linear combination of the gene expression levels weighted by the regression coefficients from the IRG signature. The model estimated the survival of every patient both in the training and testing data sets. The median risk score was used as a cutoff to separate patients to high-risk or low-risk groups. The predictive power of the IRG signature was evaluated using the area under the curve (AUC) of the receiver operating characteristic (ROC) curve at 1-year via the “survivalROC” package.

Immune cell infiltration

The association between the IRG signature and infiltrating immune cells was first analyzed in the Tumor Immune Estimation Resource (TIMER; https://cistrome.shinyapps.io/timer/), which included 6 types of immune cells (i.e., dendritic cells, macrophages, B cells, neutrophils, cluster of differentiation (CD)8+ T cells, and CD4+ T cells). The relative levels of the different immune cell types were quantified using CIBERSORT in a complex gene expression mixture. The preparation data was used in a subsequent analysis to assess each sample’s immune infractions through the CIBERSORT algorithm (R script v1.03) to estimate each cell type’s abundance in a mixed cell population using the gene expression data.

Statistical analysis

Multivariate Cox regression analyses were performed with clinical variables such as age, grade, tumor stage, tumor type and the risk score as covariables. The “survival” package was used for the Cox regression analysis. The normalization and the differential expression analyses were performed by the “Limma” package. The Wilcoxon rank-sum test, which is a non-parametric statistical hypothesis test, was used to compare 2 groups. All the statistical analyses were conducted using the R software (Version 3.6.3). A P value <0.05 was considered statistically significant. Threshold values for AUC between 0.5 and 1 are considered as better than random classifiers.

Results

Differentially expressed IRGs in ESCA

A total of 3,327 DEGs were identified between the ESCA tumor and non-tumor tissue samples using the Wilcoxon signed-rank test with cutoff of a |log2 FC| >1 and a P value <0.05 (see ). From the DEGs, 243 IRGs, including 40 downregulated genes and 203 upregulated genes, were identified. The expression patterns of the IRGs were visualized by the “Limma” package of R (see ). Further, a gene functional enrichment analysis revealed that “leukocyte migration,” “external side of plasma membrane,” and “receptor ligand activity” were the most frequent biological terms among the biological processes, cellular components, and molecular functions, respectively (see ). Additionally, the “cytokine-cytokine receptor interaction” was the most enriched pathway in the KEEG analysis (see ).

Figure 1

The DEGs. (A) A heatmap and (B) a volcano plot showing the DEGs in the ESCA tumor and non-tumor samples. (C) A heatmap and (D) a volcano plot showing the differentially expressed IRGs. The red dots indicate the highly expressed genes, and the green dots indicate the lowly expressed genes. N, normal samples; T, tumor samples; FC, fold change; DEG, differentially expressed gene; ESCA, esophageal cancer; IRG, immune-related gene.

Figure 2

The functional analysis of the differentially expressed IRGs. (A) GO analysis. The outer circles represent a scatter plot for each term of the logFC for the assigned genes. The red circles denote upregulation, and the blue circles denote downregulation. (B) The top 10 most significant genes according to the KEEG pathways. GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; FC, fold change.

TF regulatory network

A TF-involved network was used to identify hub TFs that might regulate the prognostic IRGs in ESCA. A total of 51 TFs were found to be differentially expressed between the ESCA tumor samples and non-cancerous samples (see ). Among the 51 TFs, 20 were correlated with ESCA patients’ OS (see ). These TFs and the prognostic IRGs with a correlation score greater than 0.3 was used to construct a regulatory network (see ).

Figure 3

A TF-involved regulatory network. (A) Volcano plot showing the differentially expressed TFs between the ESCA tumor samples and non-cancerous samples. The red dots indicate the highly expressed TFs, the green dots indicate the lowly expressed TFs, and the black dots indicate no difference expressed TFs. (B) The differentially expressed TFs correlated with OS. (C) The regulatory network established according to the TFs and IRGs. FC, fold change; TF, transcription factors; ESCA, esophageal cancer; OS, overall survival; IRG, immune-related gene.

The construction of the IRG prognostic signature and the predictability assessment in the training data set for ESCA

The entire group (N=158) was randomly assigned to the training data set (N=80; see Table S1) and the testing data set (N=78; see Table S2) to construct a IRG prognostic signature. The clinical data of the patients are summarized in . The univariate Cox regression analysis revealed that 10 IRGs were significantly associated with OS (P<0.05). These 10 IRGs were then included in a multivariate Cox stepwise analysis. Finally, an OS prediction gene signature was developed based on the 4 IRGs, including heat shock protein family A member 6 (HSPA6), oncostatin M (OSM), androgen receptor (AR), and nuclear receptor subfamily 2 group F member 2 (NR2F2).

Table 1

Clinical characteristics of patients with ESCA

Characteristic	Training data set (N, %)	Testing data set (N, %)	P
Age (years)			0.235
≤65	46, 57.50	52, 66.67
>65	34, 42.50	26, 33.33
Gender			0.771
Male	69, 86.25	66, 84.62
Female	11, 13.75	12, 15.38
Grade			0.716
G1–2	41, 51.25	40, 51.28
G3	20, 25.00	23, 29.49
Unknown	19, 23.75	15, 19.23
Tumor type			0.153
ADC	35, 43.75	43, 55.13
SCC	45, 56.25	35, 44.87
T			0.364
T1–2	33, 41.25	31, 39.74
T3–4	37, 46.25	42, 53.85
Unknown	10, 12.50	5, 6.41
N			0.478
N0	30, 37.50	35, 44.87
N1–3	40, 50.00	37, 47.44
Unknown	10, 12.50	6, 7.69
M			0.474
M0	57, 71.25	62, 79.49
M1	5, 6.25	3, 3.85
Unknown	18, 22.50	13, 16.66
Stage			0.351
I–II	41, 51.25	43, 55.13
III–IV	27, 33.75	29, 37.18
Unknown	12, 15.00	6, 7.69

ESCA, esophageal cancer; ADC, adenocarcinoma; SCC, squamous cell carcinoma.

Additionally, a risk score was calculated based on gene expression and the Cox regression coefficient using the following formula: (0.0079 × expression value of HSPA6) + (0.2933 × expression value of OSM) + (−3.3788 × expression value of AR) + (0.0122 × expression value of NR2F2). All the ESCA patients were divided into the high- or low-risk groups according to the median risk score. As the Kaplan-Meier curves show, the high-risk group had a worse OS than the low-risk group in the training data set (P=2.407e−03, log-rank test; see ). The AUC of the ROC curve for the IRG signature (see ) for 1-year survival was 0.808. The distribution of patient risk scores and survival times in the training data set was also plotted (see ).

Figure 4

The IRG prognostic signature for the patients with ESCA in the training data set. (A) Kaplan-Meier curve plot stratifying patients in the high-risk group with those in the low-risk group. (B) ROC curve for the OS-related prognostic signature. (C) Heatmap showing expression of IRGs in different risk group. (D) The number of patients two risk groups (green indicates low; red indicates high). (E) The scatterplots of ESCA patients with another survival status (green indicates alive; red indicates dead). IRG, immune-related gene; ESCA, esophageal cancer; ROC, receiver operating characteristic; OS, overall survival; AUC, area under the curve.

Validation of the IRG prognostic signature

The prognostic value of the 4-IRG signature was further validated in the testing data set and the full data set. Each patient was determined to be high-risk or low-risk by comparing their risk score to the median risk score calculated from the training data set. The Kaplan-Meier survival curves differed significantly between the 2 predicted groups in the testing data set (P=1.044e-02; see ) and the full data set (P=2.535e-04; see ). The ROC curve in the testing data set had an AUC of 0.657, and that in the entire data set had an AUC of 0.751 (see ). The heatmap, distribution of risk score, and survival status for the testing data set and the full data set were also illustrated (see ).

Figure 5

The validation of the IRG prognostic signatures. (A) Kaplan-Meier curves for the testing data set. (B) Kaplan-Meier curves for the full data set. (C) The AUCs of the ROC curves for predicting 1-year survival for the testing data set. (D) The AUCs of the ROC curves for predicting 1-year survival for the full data set. (E) The distribution heatmap for different risk groups for the testing data set. (F) The distribution heatmap in different risk groups for the full data set. (G) The number of patients in the other risk groups (green indicates low; red indicates high) for the testing data set. (H) The number of patients in the other risk groups (green indicates low; red indicates high) for the full data set. (I) Scatterplots of ESCA patients with different survival status (green indicates alive; red indicates dead) for the testing data set. (J) Scatterplots of ESCA patients with different survival status (green indicates alive; red indicates dead) for the full data set. ESCA, esophageal cancer; ROC, receiver operating characteristic; OS, overall survival; AUC, area under the curve.

The independence of the prognostic value of the IRG signature and nomogram of 4 IRGs

The independence of the prognostic value of the IRG signature was evaluated using other clinical variables. Multivariate Cox regression analyses were performed with age, grade, stage (T: primary tumor, N: lymph node status, M: distant metastasis), tumor type (adenocarcinoma or squamous cell carcinoma), and the risk score as covariables. The results showed that the risk score was significantly correlated with OS in both data sets (see : Training data set; : Testing data set; : Full data set), which supports its potential as a prognostic marker for ESCA. A further analysis indicated that the AUCs of the ROC curves were 0.796 (risk score), 0.594 (age), 0.583 (grade), 0.655 (stage), 0.611 (T), 0.676 (N), 0.509 (M), and 0.450 (tumor type) (see ). To establish a quantitative method for predicting ESCA prognosis, we drew a horizontal line to determine the point of each variable and calculated the total score by summing up all the points. The points were normalized to a 0 to 100 distribution. This nomogram enabled prediction of the 1-, 2-, 3-year survival according to a total score (see ).

Figure 6

Multivariate Cox regression analyses of the relationship between different clinical features and the overall survival in ESCA patients. (A) Training data set. (B) Testing data set. (C) Full data set. ESCA, esophageal cancer.

Figure 7

The nomogram and the IRG signature used to predict the accuracy of patient survival prognosis. (A) ROC curves comparing the accuracy of survival prediction by risk score and clinical characteristics. (B) Nomogram for predicting OS developed in ESCA. IRG, immune-related gene; ROC, receiver operating characteristic; AUC, area under the curve; ESCA, esophageal cancer; OS, overall survival.

Tumor immune cell infiltration in ESCA

The association between the IRG signature and tumor immune cell infiltration was explored with the TIMER database. The results showed that the 4 IRGs in the IRG signature had a strong association with the 6 types of infiltrating immune cells included in the TIMER database (see ). Using the “CIBERSORT” package, we further compared the profiles of the infiltrating immune cells between the low- and high-risk groups. We found that the infiltration levels of resting CD4 memory T cells, activated mast cells, and neutrophils in the high-risk group were significantly higher than those in the low-risk group. Conversely, the infiltration levels of resting mast cells in the high-risk group were substantially lower than those in the low-risk group (see ). However, only the infiltration levels of the activated mast cells were significantly correlated with OS in ESCA (see ). Notably, the infiltration levels of the resting CD4 memory T cells, neutrophils, and resting mast cells were not significantly correlated with OS in ESCA (see ).

Figure 8

The association between the 4 IRGs from the prognostic signature and immune cell infiltration based on data from the TIMER database. IRG, immune-related gene.

Figure 9

The landscape of immune infiltration of different risk statuses in ESCA. (A) Comparisons of the 22 significant immune fractions between the high- and low-risk groups. (B-E) Differential immune cell infiltration is significantly associated with OS in ESCA. ESCA, esophageal cancer; OS, overall survival.

Discussion

ESCA is extraordinarily malignant, has a poor prognosis, and caused approximately 544076 deaths worldwide in 2020 and its incidence continues to increase (19). Cancer immunotherapy is a novel therapeutic strategy that has made remarkable advances in recent years (20). The significance of IRGs in cancer progression and immunotherapy has been well established; however, their roles in ESCA remains poorly defined. It has been reported that infection and inflammation account for approximately 25% of cancer-causing factors (21). Esophagus inflammation, such as idiopathic achalasia and Barrett’s esophagus, has been shown to increase the risk of ESCA (22).

In this study, the functional enrichment analysis revealed that the different IRGs in primary ESCA were the most frequently implicated in inflammatory pathways, including cell chemotaxis, leukocyte chemotaxis, and leukocyte migration. We speculated that the initiation of ESCA often occurs in densely infiltrated inflammatory environments, which are a consistent hallmark of cancer (23). We studied the changes in the immune genome profile in primary ESCA and uncovered its possible underlying molecular mechanisms. The KEEG analysis revealed that these IRGs were enriched in the cytokine-cytokine receptor interaction pathway, which reflects previous findings that these pathways play a crucial role in the proliferation, immune responses, and progression in several cancers (24-26). The dysregulation of cytokine interactions is also involved in the pathogenesis of ESCA (27). We constructed an IRG-TF regulatory network to examine the underlying mechanisms of primary ESCA development and found 20 TFs related to the prognostic IRGs. Our study provides evidence of an inflammatory response during the initiation of ESCA, which correlates with the activation of immune-related pathways.

Due to tumor heterogeneity, there are individual differences in patients responses to immunotherapy. Thus, identifying robust gene signatures that correspond to cancer patients’ immune status is essential in identifying reliable prognostic biomarkers and stratifying patients into high- or low-risk groups to optimize their responsiveness to immunotherapy. To date, IRG signatures have been explored for several cancers, including lung cancer (28), ovarian cancer (29), papillary thyroid cancer (30), and renal papillary cell carcinoma (31). Recently, Wang et al. examined the IRG signature for ESCA (32); however, metastatic ESCA was not excluded from that study. The immunogenomic analysis of primary ESCA remains rare. Our study integrated the analysis of IRGs in primary ESCA to explore the clinical significance and molecular characteristics of IRGs in ESCA. We developed a robust IRG prognostic signature, consisting of 4 IRGs, for primary ESCA. This prognostic signature was then validated using testing and full data sets. The Kaplan-Meier curves and ROC curves revealed that the prognostic signature had accurate predictive performance. Based on the multivariate Cox regression analyses, our study indicated that the 4-IRG signature could serve as an independent prognostic biomarker for patients with primary ESCA.

All 4 IRGs in the prognostic signature are involved in cancer development. NR2F2 and OSM have been shown to be related to metastasis in gastric cancer and breast cancer, respectively (33,34); HSPA6 has been shown to play a critical role in the recurrence of hepatocellular carcinoma (35); AR has been shown to play a crucial role in cancer progression and continues to be a primary therapeutic target in prostate cancer (36). However, these previous studies (33-36) included limited direct information about the IRGs of primary ESCA. Thus, these genes’ immunomodulatory roles in monitoring the progression of primary ESCA and predicting prognosis need to be further investigated.

A well-known hallmark of cancer is tumor cells’ ability to escape immune destruction. Immune destruction not only affects cancer cells but also affects the infiltrating immune cells (37). We used the TIMER database to examine the relationship between the IRG signature and the infiltrating immune cells of primary ESCA patients. Our results showed that our IRG signature had a relationship with several immune cells. However, the database was limited to 6 types of immune cells. We used the CIBERSORT algorithm to further evaluate the immune status and correlate clinical outcomes in primary ESCA patients. Our study showed significantly different infiltrating immune cell profiles for the high- and low-risk groups of ESCA patients. The infiltration levels of the activated mast cells were significantly correlated with OS, and the high-risk group had a much higher level than the low-risk group.

Previous studies have revealed that mast cells had pro- or anti-tumorigenic roles depending on the tumor types, the stage, and localization within the tumor. Activated mast cells have been shown to be associated with poor prognosis in thyroid cancer (38), gastric cancer (39), hepatic cellular carcinoma (40), and colorectal cancer (41). Mast cell density is increased in gastric cancer, and correlated with angiogenesis, the number of metastatic lymph nodes, and patient survival (42). Conversely, it is associated with a favorable prognosis in prostate cancer (43). A previous study revealed that interleukin-33 mediates mast cell activation and further promotes cancer progression in gastric cancer (44), while in pancreatic cancer, mast cells secrete tryptase and promote the growth of cancer via the activation of angiopoietin-1 (45). In this study, we found that infiltrating activated mast cells in the primary ESCA tumor environment is associated with a poor prognosis. However, the underlying mechanism remains unknown.

Our research had several limitations. We verified the prediction model with testing and full data sets; however, further in-vitro and in-vivo studies and studies using data from other primary ESCA patient cohorts need to be conducted to confirm the model’s accuracy. Additionally, other potential prognostic variables, such as lymphovascular invasion and the neutrophil-to-lymphocyte ratio reported to be correlated with OS in ESCA, were not included in our study. Finally, the action and mechanisms of the 4 IRGs in primary ESCA should be elucidated in the future.

Conclusions

We constructed and verified a novel 4-IRG signature as a prognostic biomarker for primary ESCA patients. Our study also identified novel targets for immunotherapy and individualized therapy in ESCA.

The article’s supplementary files as