Abstract

Background: Esophageal cancer (ESCA) is one of the most commonly diagnosed cancers in the world. Tumor immune microenvironment is closely related to tumor prognosis. The present study aimed at analyzing the competing endogenous RNA (ceRNA) network and tumor-infiltrating immune cells in ESCA.

Methods: The expression profiles of mRNAs, lncRNAs, and miRNAs were downloaded from the Cancer Genome Atlas database. A ceRNA network was established based on the differentially expressed RNAs by Cytoscape. CIBERSORT was applied to estimate the proportion of immune cells in ESCA. Prognosis-associated genes and immune cells were applied to establish prognostic models basing on Lasso and multivariate Cox analyses. The survival curves were constructed with Kaplan–Meier method. The predictive efficacy of the prognostic models was evaluated by the receiver operating characteristic (ROC) curves.

Results: The differentially expressed mRNAs, lncRNAs, and miRNAs were identified. We constructed the ceRNA network including 23 lncRNAs, 19 miRNAs, and 147 mRNAs. Five key molecules (HMGB3, HOXC8, HSPA1B, KLHL15, and RUNX3) were identified from the ceRNA network and five significant immune cells (plasma cells, T cells follicular helper, monocytes, dendritic cells activated, and neutrophils) were selected via CIBERSORT. The ROC curves based on key genes and significant immune cells all showed good sensitivity (AUC of 3-year survival: 0.739, AUC of 5-year survival: 0.899, AUC of 3-year survival: 0.824, AUC of 5-year survival: 0.876). There was certain correlation between five immune cells and five key molecules.

Conclusion: The present study provides an effective bioinformatics basis for exploring the potential biomarkers of ESCA and predicting its prognosis.

Introduction

Esophageal cancer (ESCA) is one of the leading triggers for cancer-related death in the world, and its incidence rate continues to rise [1]. The 5-year survival rate of ESCA is less than 15%, so it has poor prognosis and high mortality [2]. Due to the lack of effective early diagnosis, ESCA usually appears in the advanced stage. Patients can no longer swallow solid foods, and the clinical prognosis is poor [3]. Accurately assessing the prognosis of patients can provide them with personalized treatment and reduce the mortality of patients [4]. Therefore, further research is necessary to discover potential biomarkers to improve the diagnosis and prognosis.

Nowadays, immunotherapy and targeted treatment are widely used, which have become a significant way to improve the prognosis of ESCA patients [5]. The tumor microenvironment (TME) plays an important role in the formation, development, and treatment of tumors [6]. Immune cells are the major component of TME and closely associated with tumor occurrence. In addition, it has been proved that the prognosis and malignant degree of tumors are associated with the TME immune cells to a certain extent. The competing endogenous RNA (ceRNA) network was composed of mRNAs, lncRNAs, and miRNAs. The complex crosstalk of the ceRNA network has been verified in numerous diseases including cancer [7]. For example, it is reported that down-regulated lncRNA UCA1 acted as ceRNA to adsorb microRNA-498 to suppress proliferation and invasion of esophageal cancer cells by inhibiting ZEB2 expression [8]. Moreover, the crosstalk between the tumor cells and tumor-infiltrating immune cells is usually regulated via ceRNA networks [9].

In the present study, we obtained the expression profiles of mRNAs, lncRNAs, and miRNAs from the Cancer Genome Atlas (TCGA) database and established a ceRNA network on the basis of differentially expressed genes (DEGs) between tumor and normal samples. The fraction of different immune cells in ESCA was estimated by CIBERSORT. Then, we constructed two prognostic models according to the results of the ceRNA network and CBERSORT analysis. We further evaluated the relationship between the immune cells and the hub genes involved in the ceRNA network to find potential molecular pathways for the immunotherapy of ESCA patients.

Materials and methods

Data source

Gene expression profiles and the related clinical data of TCGA-ESCA patients were downloaded from the TCGA database (data release 23.0- April 7, 2020), including mRNA, lncRNA, and miRNA. We downloaded the FPKM value of gene expression of 160 ESCA samples and 11 normal samples. After excluding the samples that had missing clinical information, 159 samples with complete clinical information were finally obtained. Data acquisition was in accordance with TCGA publication guidelines.

Differential gene expression analysis and construction of the ceRNA network

Using the edgeR package in R software, we identified differentially expressed mRNA (DEmRNA), miRNA (DEmiRNA), and lncRNA (DElncRNA) [10]. Genes with |logFC| > 1 and P.adjust < 0.05 were regarded as DEGs. MiRcode database was used to predict lncRNA–miRNA interactions [11]. Three databases, including miRTarBase, TargetScan, and miRDB, were performed to find the potential mRNAs for DEmiRNA [12–14]. Only the genes included in all three datasets were selected for further analysis. By crossing the predicted mRNAs with DEmRNAs, we obtained the mRNAs in the network. According to the results above, a ceRNA network was constructed and further visualized by Cytoscape [15].

Functional and pathway enrichment analyses

To explore potential biological processes related to the DEmRNAs involved in the ceRNA network, we performed gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis using the clusterProfiler R package [16]. The P.adjust < 0.05 displays statistical significance. The most significant GO terms and pathways were selected for visualization.

Construction of a five-gene prognostic model based on the ceRNA network

The DEGs correlated with overall survival (OS) were screened via univariate Cox regression analysis with the survival R package. Then, lasso regression was employed to reduce dimensionality, and the most robust markers were selected to construct the prognostic model. The 10-fold cross validation was applied to define the optimal value of the penalty parameter λ [17]. Each patient’s risk score was finally calculated by multiplying the expression level of each gene by its Cox regression coefficient and adding these values.

Immune infiltrate analysis

CIBERSORT is a calculation method for characterizing cell subsets of interest in high-dimensional genomic data derived from large tissue samples [18]. To explore the difference of infiltrating immune cells, we used CIBERSORT to evaluate the fraction of 22 kinds of immune cells in the ESCA samples and normal samples with a deconvolution algorithm [19]. Only samples with CIBERSORT P<0.05 could be used for further analysis.

Construction of a prognostic model based on immune cells

We used univariate Cox regression analysis to determine the prognostic value of the 22 immune cells. Then, the key members of the immune cells were identified via Lasso and multivariate Cox analyses. The risk score of each patient was finally calculated via multiplying the abundance of each immune cell by its Cox regression coefficient and adding these values.

The Co-expression and TIMER analysis

Pearson correlation analysis was performed to analyze the co-expression patterns between the immune cells and the key molecules involved in the ceRNA network. As a public resource, TIMER uses a deconvolution method to infer the abundance of immune infiltration from gene expression profiles [20]. Based on the TIMER website, we evaluated the correlation between the key molecules involved in the ceRNA network and six TME infiltrating cells including CD4+ T cells, dendritic cells, B cells, CD8+ T cells, neutrophils, and macrophages.

Statistical analysis

Kaplan–Meier method was used to construct survival curves. The receiver operating characteristic (ROC) curve was plotted with the SurvivalROC R package. Statistical P<0.05 was considered statistically significant. R software (version 3.6.3) was used to perform all data analyses.

Results

Identification of DEGs in ESCA

A total of 3211 DEmRNAs, 381 DElncRNAs, and 144 DEmiRNAs were identified using the edgeR package in R software. Among them, 1461 lncRNAs, 19 miRNAs, and 1511 mRNAs were up-regulated and 554 lncRNAs, 28 miRNAs, and 803 mRNAs were down-regulated. The cut-off criteria was |logFC| >1 and P.adjust<0.05. Figure 1A–F showed the heatmaps and volcano plots.

Differentially expressed genes between ESCA samples and normal samples

Figure 1
Differentially expressed genes between ESCA samples and normal samples

(A) Heatmap and (B) volcano plot of differentially expressed mRNAs between normal and tumor samples. (C) Heatmap and (D) volcano plot of differentially expressed lncRNAs between normal and tumor samples. (E) Heatmap and (F) volcano plot of differentially expressed miRNAs between normal and tumor samples. Red points represent up-regulated genes. Green points represent down-regulated genes. Black points represent genes with no significant difference; ESCA, esophageal cancer.

Figure 1
Differentially expressed genes between ESCA samples and normal samples

(A) Heatmap and (B) volcano plot of differentially expressed mRNAs between normal and tumor samples. (C) Heatmap and (D) volcano plot of differentially expressed lncRNAs between normal and tumor samples. (E) Heatmap and (F) volcano plot of differentially expressed miRNAs between normal and tumor samples. Red points represent up-regulated genes. Green points represent down-regulated genes. Black points represent genes with no significant difference; ESCA, esophageal cancer.

The ceRNA network building and functional enrichment analysis

With miRcode database, we predicted 115 pairs of interacting lncRNA–miRNA. Then, 1043 potential mRNAs included in all three datasets (miRTarBase, TargetScan, and miRDB) were predicted for DEmiRNAs. Finally, 147 target mRNAs were gained by intersecting the predicted mRNAs with DEmRNAs. To investigate how lncRNAs mediate mRNA by binding to miRNA in ESCA, we constructed the ceRNA network according to the results above, including 23 lncRNA nodes, 19 miRNA nodes, and 147 mRNA nodes (Figure 2A).

The ceRNA network and functional annotation

Figure 2
The ceRNA network and functional annotation

(A) The ceRNA network of differentially expressed mRNAs, differentially expressed miRNAs, and differentially expressed lncRNAs; Red, up-regulated; Blue, down-regulated. Circle represents mRNA; Diamond represents lncRNA; Rectangle represents miRNA. (B) GO enrichment analyses of the mRNAs involved in the ceRNA network. (C) KEGG pathway enrichment analyses of the mRNAs involved in the ceRNA network; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Figure 2
The ceRNA network and functional annotation

(A) The ceRNA network of differentially expressed mRNAs, differentially expressed miRNAs, and differentially expressed lncRNAs; Red, up-regulated; Blue, down-regulated. Circle represents mRNA; Diamond represents lncRNA; Rectangle represents miRNA. (B) GO enrichment analyses of the mRNAs involved in the ceRNA network. (C) KEGG pathway enrichment analyses of the mRNAs involved in the ceRNA network; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

To gain further insight into the biological functions and potential pathways bound up with the ceRNA network, GO and KEGG pathway analyses were performed by clusterProfiler R package. The biological processes of GO term enrichment demonstrated that the DEmRNAs involved in the ceRNA network were particularly enriched in mesenchymal cell differentiation, epithelial cell proliferation, and cell cycle checkpoint (Figure 2B and Supplementary Table S1). In addition, the results of KEGG analysis suggested that the DEG-associated pathways were cellular senescence, p53 signaling pathway, and cell cycle (Figure 2C and Supplementary Table S2).

Construction of a five-gene prognostic model based on the ceRNA network

Univariate Cox regression analysis was applied to screen the prognostic-related DEmRNAs involved in the ceRNA network. The results demonstrated that 8 of the 147 DEmRNAs were significantly correlated with OS. Then, Lasso Cox regression was used to identify hub survival-associated genes by reducing the dimensionality (Figure 3B,C). Finally, multivariate Cox regression analysis was performed to calculate the relative coefficients of the genes with the best prognostic value. The result indicated that HMGB3, HSPA1B, KLHL15, and RUNX3 were independent prognostic factors (P<0.05) (Figure 3A). Supplementary Table S3 listed the coefficients for each gene. The prognostic risk score for each patient was calculated. Then, patients were separated into high and low-risk groups basing on the median risk score. The survival analysis showed worse prognosis for patients in high-risk group (Figure 3F). ROC curve was performed to determine the predictive accuracy of this prognostic signature, the AUC value was 0.739 for 3-year survival and 0.899 for 5-year survival (Figure 3G). This revealed that this prognostic signature had good prediction accuracy for ESCA patients. The distribution of survival state showed that the number of deaths increased with the increase of risk score (Figure 3D). Heatmap showed the trend of gene expression in high- and low-risk groups (Figure 3E).

Construction of a five-gene prognostic model

Figure 3
Construction of a five-gene prognostic model

(A) The results of the multivariate Cox regression. (B) Least absolute shrinkage and selection operator (LASSO) coefficient profiles of the five key molecules. (C) Tuning parameter selection by 10-fold cross-validation in the LASSO model. The partial likelihood deviance was plotted against log(Lambda/λ), and λ was the tuning Parameter. (D) Risk score distribution and survival status. The black dotted line is the optimum cut-off dividing patients into low-risk and high-risk groups. The red curve represents high risk and the green curve represents low risk. The dots indicate the survival status, the red dot indicates the death of the patient, and the green dot indicates alive. (E) Expression heat map. (F) Survival curve for the low-risk and high-risk groups. (G) ROC curve analysis in the TCGA-ESCA cohort; ESCA, esophageal cancer; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas.

Figure 3
Construction of a five-gene prognostic model

(A) The results of the multivariate Cox regression. (B) Least absolute shrinkage and selection operator (LASSO) coefficient profiles of the five key molecules. (C) Tuning parameter selection by 10-fold cross-validation in the LASSO model. The partial likelihood deviance was plotted against log(Lambda/λ), and λ was the tuning Parameter. (D) Risk score distribution and survival status. The black dotted line is the optimum cut-off dividing patients into low-risk and high-risk groups. The red curve represents high risk and the green curve represents low risk. The dots indicate the survival status, the red dot indicates the death of the patient, and the green dot indicates alive. (E) Expression heat map. (F) Survival curve for the low-risk and high-risk groups. (G) ROC curve analysis in the TCGA-ESCA cohort; ESCA, esophageal cancer; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas.

Composition of immune cells in ESCA

The abundance of 22 kinds of immune cells was estimated to explore the difference of tumor microenvironment between ESCA samples and normal samples (Supplementary Table S5). Figure 4A,B showed the ratio of 22 immune cells detected via the CIBERSORT algorithm. The violin plot was performed by the Wilcoxon test. The results suggested that the infiltration levels of macrophages M0 (P=0.008) and dendritic cells activated (P=0.032) were obviously higher in ESCA samples, while CD8 T cells (P=0.039), CD4 T cells memory resting (P=0.023), monocytes (P=0.01), and mast cells resting (P<0.001) were significantly lower (Figure 4C).

Evaluation of proportions of tumor-infiltrating immune cells based on CIBERSORT

Figure 4
Evaluation of proportions of tumor-infiltrating immune cells based on CIBERSORT

The (A) composition and (B) heat map of 22 subtypes of immune cells in ESCA. X axis represents ESCA samples. Y axis represents the proportion and abundance of immune cells. (C) Differences of 22 tumor-infiltrating immune cells between normal and tumor groups; ESCA, esophageal cancer.

Figure 4
Evaluation of proportions of tumor-infiltrating immune cells based on CIBERSORT

The (A) composition and (B) heat map of 22 subtypes of immune cells in ESCA. X axis represents ESCA samples. Y axis represents the proportion and abundance of immune cells. (C) Differences of 22 tumor-infiltrating immune cells between normal and tumor groups; ESCA, esophageal cancer.

Construction of a prognostic model on the base of immune cells

Univariate Cox regression analysis was used to screen the immune cells with prognostic value. Five immune cells including Plasma cells, T cells follicular helper, monocytes, dendritic cells activated, and neutrophils were identified by Lasso and multivariate Cox analyses (Figure 5A–C). A new prognostic signature was built basing on the five immune cells. The relative coefficients were calculated by multivariate Cox regression analysis and listed in Supplementary Table S4. Then, we calculated the patient’s risk score and separated them into high- and low-risk groups according to the median risk score. The survival analysis showed worse prognosis for patients in high-risk group (Figure 4E). ROC curve indicated good prediction accuracy with the AUC of 0.824 for 3-year survival and 0.876 for 5-year survival (Figure 4F). The survival status distribution graph indicated that the number of deaths increased with the increase of risk score (Figure 5D).

Construction of a prognostic model based on the immune cells

Figure 5
Construction of a prognostic model based on the immune cells

(A) The results of the multivariate Cox regression. (B) Least absolute shrinkage and selection operator (LASSO) coefficient profiles of the five immune cells. (C) Tuning parameter selection by 10-fold cross-validation in the LASSO model. The partial likelihood deviance was plotted against log(Lambda/λ), and λ was the tuning Parameter. (D) Risk score distribution and survival status. The black dotted line is the optimum cut-off dividing patients into low-risk and high-risk groups. The red curve represents high risk and the green curve represents low risk. The dots indicate the survival status, the red dot indicates the death of the patient, and the green dot indicates alive. (E) Survival curve and (F) ROC curve analyses; ROC, receiver operating characteristic.

Figure 5
Construction of a prognostic model based on the immune cells

(A) The results of the multivariate Cox regression. (B) Least absolute shrinkage and selection operator (LASSO) coefficient profiles of the five immune cells. (C) Tuning parameter selection by 10-fold cross-validation in the LASSO model. The partial likelihood deviance was plotted against log(Lambda/λ), and λ was the tuning Parameter. (D) Risk score distribution and survival status. The black dotted line is the optimum cut-off dividing patients into low-risk and high-risk groups. The red curve represents high risk and the green curve represents low risk. The dots indicate the survival status, the red dot indicates the death of the patient, and the green dot indicates alive. (E) Survival curve and (F) ROC curve analyses; ROC, receiver operating characteristic.

The Co-expression analysis

We assessed possible correlations between 22 kinds of immune cells (Figure 6A), and Figure 6B showed the important co-expression patterns between the five immune cells and the five key molecules involved in the ceRNA network. The results showed that HMGB3 was positively correlated with follicular helper T cells (R = 0.16, P=0.049). HOXC8 was negatively associated with plasma cells (R = −0.18, P=0.024). RUNX3 was positively associated with monocytes (R = 0.18, P=0.021). HSPA1B was negatively associated with monocytes (R = −0.17, P=0.032) and plasma cells (R = −0.18, P=0.022) (Supplementary Figure S1). Besides, the correlations obtained from the TIMER database are generally consistent with the results above (Figure 6C–G).

Correlation analysis of the key genes and immune cells

Figure 6
Correlation analysis of the key genes and immune cells

(A) Correlation heat map of all immune infiltration cells. (B) Correlation heat map of key genes and immune cells involved in the prognostic models. (C–G) Correlation between key genes expression and six types of immune cells in the TIMER database.

Figure 6
Correlation analysis of the key genes and immune cells

(A) Correlation heat map of all immune infiltration cells. (B) Correlation heat map of key genes and immune cells involved in the prognostic models. (C–G) Correlation between key genes expression and six types of immune cells in the TIMER database.

Discussion

ESCA is one of the most commonly diagnosed cancers in the world which is associated with poor prognosis and high mortality rate [21]. Thus, exploration of effective molecules and prediction models of ESCA could be helpful for its diagnosis and treatment.

In the present study, we focused on the ceRNA network and TME infiltrating immune cells. We identified 3211 DEmRNAs, 381 DElncRNAs, and 144 DEmiRNAs based on the TCGA database. A ceRNA network was constructed. Five key molecules were identified from the ceRNA network and five significant immune cells were identified via CIBERSORT. Two risk models were established on the base of the significant TME infiltrating immune cells and mRNAs. The high AUC values of both prognostic models indicated the advantage in predicting OS of ESCA patients. GO term and KEGG pathway enrichment analysis suggested that the DEmRNAs involved in the ceRNA network were mostly involved in cell cycle checkpoint and p53 signaling pathway. Previous studies revealed that dysregulation of cell cycle process played a vital role in the progression of tumors [22]. P53 signaling pathway is a typical tumor suppressor pathway, which is related to many types of cancer [23]. It has been reported that SASS6 can promote the proliferation of esophageal squamous carcinoma cells by inhibiting the p53 signaling pathway [24]. Among the ceRNA network, five significant molecules including HMGB3, HOXC8, HSPA1B, KLHL15, and RUNX3 could be considered as independent prognostic factors according to the results of multivariate Cox regression. HMGB3, high mobility group-box 3, is a member of the high mobility group box subfamily [25]. It has been proved that HMGB3 can affect the occurrence and development of many tumors. A recent study showed that HMGB3 was up-regulated in breast cancer, and silencing HMGB3 can inhibit breast cancer cell proliferation and tumor growth [26]. HMGB3 has also been reported as a potential prognostic marker in colorectal cancer. It can promote growth and migration of colorectal cancer via activating WNT/β-catenin pathway [27]. In addition, the high expression of HMGB3 is also an important prognostic indicator of low survival rate in patients with esophageal cancer [28]. HOXC8 is one of the HOX family genes. The HOX genes are significant factors in the regulation of embryogenesis. They encode a set of transcription factors and regulate the expression of downstream target genes through specific DNA binding [29]. In the past decade, HOX genes have been proved to be dysregulated in many solid tumors [30]. Furthermore, over-expression of HOX gene was related to poor prognosis [31–33]. HSPA1B is a member of the heat shock protein 70 (HSP70) family. It is reported that the expression of HSPs was higher in various tumors and closely associated with tumor progression [34]. HSPA1B could be considered as a promising target for prognostic prediction in esophageal cancer [35]. RUNX3, runt-related transcription factor 3, is a tumor suppressor gene involved in the TGF-β signaling pathway [36]. Previous studies have indicated that RUNX3 is a tumor suppressor gene for several human cancers including esophageal cancer [37]. The low expression level of RUNX3 showed worse prognosis for ESCA patients [38]. The protein encoded by KLHL15 may be involved in protein ubiquitination and cytoskeletal organization. Nevertheless, the specific role of KLHL15 in ESCA needs to be further explored.

Then, to explore the difference of tumor microenvironment between ESCA samples and normal samples, we observed the immune cell composition between normal and tumor groups. We found that Macrophages M0 and dendritic cells activated were obviously higher in ESCA samples, while CD8 T cells and CD4 T cells memory resting were obviously higher in normal samples. Tumor-associated macrophages, a major immune component of a variety of tumors, could promote tumor progression via secreting pro-angiogenic and growth factors [39,40]. Macrophages M0 can be polarized into different phenotypes, including macrophages M1 and macrophages M2 [41]. Some studies have indicated that macrophages M1 participated in pathogen clearance and proinflammatory response, while macrophages M2 were anti-inflammatory and correlated with tumor progression [42]. Furthermore, it is reported that high infiltration of tumor-associated macrophages is associated with worse prognosis of ESCA patients [43,44]. Dendritic cells are potent antigen-presenting cells, and the infiltration of dendritic cells in a tumor is considered the host’s immune response to the tumor [45]. They provide costimulatory molecules and cytokines, which provide necessary signals for the activation and differentiation of T cells, thus form the immune response [46]. In addition, they can activate anti-tumor responses via interacting with other immune cells, including NK cells and B cells [47,48].

T cells, including CD4 and CD8, were important components of TME and closely correlated with the development and prognosis of tumors [49]. Cytotoxic T lymphocytes (CTLs), which mainly express T cell co-receptor CD8, are closely associated with anti-tumor immune response [50]. It has been reported that CD8 T cells in esophageal cancer were related to positive OS [51]. Moreover, a recent study showed that high CD8 T cells and CD4 T cells infiltration have the potential to be prognostic markers and predict better OS in ESCA patients [52].

However, our study is limited because we need a larger sample size to validate our results, and the two prognostic models need more samples to further confirm its predictive ability. Furthermore, the present study is only a multi-dimensional correlation research, thus further experiments are necessary for better understanding the functions of the hub genes in ESCA.

Conclusions

On the basis of tumor-infiltrating immune cells and the ceRNA network, two predictive models were constructed to predict the prognosis of ESCA patients. Five key molecules including HMGB3, HOXC8, HSPA1B, KLHL15, and RUNX3 could be considered as independent prognostic factors. In addition, the significant immune cells described in the study might play a vital role in the occurrence and development of ESCA. In conclusion, the present study provides an effective bioinformatics basis for exploring the molecular mechanism of esophageal cancer and predicting its prognosis.

Data Availability

All data were downloaded from TCGA (https://portal.gdc.cancer.gov/).

Competing Interests

The authors declare that there are no competing interests associated with the manuscript.

Funding

This work is supported by the Science and Technology Project of Nantong City [grant number JC2020012].

Author Contribution

Yuhua Chen and Hao Zhou conceived and designed the study. Zhendong Wang, Zhanghao Huang, Jinjie Wang, and Miaosen Zheng performed the data analysis. Hao Zhou wrote the manuscript. Xuejun Ni and Lei Liu revised the manuscript. All authors read and approved the final manuscript. Yuhua Chen and Hao Zhou contributed equally to this work.

Abbreviations

     
  • ceRNA

    competing endogenous RNA

  •  
  • DEG

    differentially expressed gene

  •  
  • ESCA

    esophageal cancer

  •  
  • LASSO

    least absolute shrinkage and selection operator

  •  
  • OS

    overall survival

  •  
  • ROC

    receiver operating characteristic

  •  
  • TCGA

    The Cancer Genome Atlas

References

1.
Siegel
R.L.
,
Miller
K.D.
and
Jemal
A.
(
2020
)
Cancer statistics, 2020
.
CA Cancer J. Clin.
70
,
7
30
[PubMed]
2.
Domper Arnal
M.J.
,
Ferrandez Arenas
A.
and
Lanas Arbeloa
A.
(
2015
)
Esophageal cancer: risk factors, screening and endoscopic treatment in Western and Eastern countries
.
World J. Gastroenterol.
21
,
7933
7943
[PubMed]
3.
Dong
Z.
,
Zhang
H.
,
Zhan
T.
and
Xu
S.
(
2018
)
Integrated analysis of differentially expressed genes in esophageal squamous cell carcinoma using bioinformatics
.
Neoplasma
65
,
523
531
[PubMed]
4.
Li
J.
,
Li
X.
,
Zhang
C.
,
Zhang
C.
and
Wang
H.
(
2020
)
A signature of tumor immune microenvironment genes associated with the prognosis of nonsmall cell lung cancer
.
Oncol. Rep.
43
,
795
806
[PubMed]
5.
Wang
L.
,
Wei
Q.
,
Zhang
M.
et al.
(
2020
)
Identification of the prognostic value of immune gene signature and infiltrating immune cells for esophageal cancer patients
.
Int. Immunopharmacol.
87
,
106795
[PubMed]
6.
Wu
T.
and
Dai
Y.
(
2017
)
Tumor microenvironment and therapeutic response
.
Cancer Lett.
387
,
61
68
[PubMed]
7.
Kristensen
L.S.
,
Hansen
T.B.
,
Veno
M.T.
and
Kjems
J.
(
2018
)
Circular RNAs in cancer: opportunities and challenges in the field
.
Oncogene
37
,
555
565
[PubMed]
8.
Wang
P.
,
Liu
X.
,
Han
G.
et al.
(
2019
)
Downregulated lncRNA UCA1 acts as ceRNA to adsorb microRNA-498 to repress proliferation, invasion and epithelial mesenchymal transition of esophageal cancer cells by decreasing ZEB2 expression
.
Cell Cycle
18
,
2359
2376
[PubMed]
9.
Salmena
L.
,
Poliseno
L.
,
Tay
Y.
,
Kats
L.
and
Pandolfi
P.P.
(
2011
)
A ceRNA hypothesis: the Rosetta Stone of a hidden RNA language?
Cell
146
,
353
358
[PubMed]
10.
Robinson
M.D.
,
McCarthy
D.J.
and
Smyth
G.K.
(
2010
)
edgeR: a bioconductor package for differential expression analysis of digital gene expression data
.
Bioinformatics
26
,
139
140
[PubMed]
11.
Jeggari
A.
,
Marks
D.S.
and
Larsson
E.
(
2012
)
miRcode: a map of putative microRNA target sites in the long non-coding transcriptome
.
Bioinformatics
28
,
2062
2063
[PubMed]
12.
Chou
C.H.
,
Shrestha
S.
,
Yang
C.D.
et al.
(
2018
)
miRTarBase update 2018: a resource for experimentally validated microRNA-target interactions
.
Nucleic Acids Res.
46
,
D296
D302
[PubMed]
13.
Wong
N.
and
Wang
X.
(
2015
)
miRDB: an online resource for microRNA target prediction and functional annotations
.
Nucleic Acids Res.
43
,
D146
D152
[PubMed]
14.
Agarwal
V.
,
Bell
G.W.
,
Nam
J.W.
and
Bartel
D.P.
(
2015
)
Predicting effective microRNA target sites in mammalian mRNAs
.
Elife
4
,
e05005
15.
Shannon
P.
,
Markiel
A.
,
Ozier
O.
et al.
(
2003
)
Cytoscape: a software environment for integrated models of biomolecular interaction networks
.
Genome Res.
13
,
2498
2504
[PubMed]
16.
Yu
G.
,
Wang
L.G.
,
Han
Y.
and
He
Q.Y.
(
2012
)
clusterProfiler: an R package for comparing biological themes among gene clusters
.
OMICS
16
,
284
287
[PubMed]
17.
Gao
J.
,
Kwan
P.W.
and
Shi
D.
(
2010
)
Sparse kernel learning with LASSO and Bayesian inference algorithm
.
Neural. Netw.
23
,
257
264
[PubMed]
18.
Chen
B.
,
Khodadoust
M.S.
,
Liu
C.L.
,
Newman
A.M.
and
Alizadeh
A.A.
(
2018
)
Profiling tumor infiltrating immune cells with CIBERSORT
.
Methods Mol. Biol.
1711
,
243
259
[PubMed]
19.
Newman
A.M.
,
Liu
C.L.
,
Green
M.R.
et al.
(
2015
)
Robust enumeration of cell subsets from tissue expression profiles
.
Nat. Methods
12
,
453
457
[PubMed]
20.
Li
T.
,
Fan
J.
,
Wang
B.
et al.
(
2017
)
TIMER: A Web Server for comprehensive analysis of tumor-infiltrating immune cells
.
Cancer Res.
77
,
e108
e110
[PubMed]
21.
Bray
F.
,
Ferlay
J.
,
Soerjomataram
I.
et al.
(
2018
)
Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries
.
CA Cancer J. Clin.
68
,
394
424
[PubMed]
22.
Tripathi
V.
,
Shen
Z.
,
Chakraborty
A.
et al.
(
2013
)
Long noncoding RNA MALAT1 controls cell cycle progression by regulating the expression of oncogenic transcription factor B-MYB
.
PLos Genet.
9
,
e1003368
[PubMed]
23.
Labuschagne
C.F.
,
Zani
F.
and
Vousden
K.H.
(
2018
)
Control of metabolism by p53 - Cancer and beyond
.
Biochim. Biophys. Acta Rev. Cancer
1870
,
32
42
[PubMed]
24.
Xu
Y.
,
Zhu
K.
,
Chen
J.
et al.
(
2020
)
SASS6 promotes proliferation of esophageal squamous carcinoma cells by inhibiting the p53 signaling pathway
.
Carcinogenesis
42
,
254
262
25.
Agresti
A.
and
Bianchi
M.E.
(
2003
)
HMGB proteins and gene expression
.
Curr. Opin. Genet. Dev.
13
,
170
178
[PubMed]
26.
Gu
J.
,
Xu
T.
,
Huang
Q.H.
,
Zhang
C.M.
and
Chen
H.Y.
(
2019
)
HMGB3 silence inhibits breast cancer cell proliferation and tumor growth by interacting with hypoxia-inducible factor 1alpha
.
Cancer Manag. Res.
11
,
5075
5089
[PubMed]
27.
Zhang
Z.
,
Chang
Y.
,
Zhang
J.
et al.
(
2017
)
HMGB3 promotes growth and migration in colorectal cancer by regulating WNT/beta-catenin pathway
.
PLoS ONE
12
,
e0179741
[PubMed]
28.
Gao
J.
,
Zou
Z.
,
Gao
J.
et al.
(
2015
)
Increased expression of HMGB3: a novel independent prognostic marker of worse outcome in patients with esophageal squamous cell carcinoma
.
Int. J. Clin. Exp. Pathol.
8
,
345
352
[PubMed]
29.
Grier
D.G.
,
Thompson
A.
,
Kwasniewska
A.
et al.
(
2005
)
The pathophysiology of HOX genes and their role in cancer
.
J. Pathol.
205
,
154
171
[PubMed]
30.
Abate-Shen
C.
(
2002
)
Deregulated homeobox gene expression in cancer: cause or consequence?
Nat. Rev. Cancer
2
,
777
785
[PubMed]
31.
Morgan
R.
,
Simpson
G.
,
Gray
S.
et al.
(
2016
)
HOX transcription factors are potential targets and markers in malignant mesothelioma
.
BMC Cancer
16
,
85
[PubMed]
32.
Kelly
Z.
,
Moller-Levet
C.
,
McGrath
S.
et al.
(
2016
)
The prognostic significance of specific HOX gene expression patterns in ovarian cancer
.
Int. J. Cancer
139
,
1608
1617
[PubMed]
33.
Shen
L.Y.
,
Zhou
T.
,
Du
Y.B.
,
Shi
Q.
and
Chen
K.N.
(
2019
)
Targeting HOX/PBX dimer formation as a potential therapeutic option in esophageal squamous cell carcinoma
.
Cancer Sci.
110
,
1735
1745
[PubMed]
34.
Lianos
G.D.
,
Alexiou
G.A.
,
Mangano
A.
et al.
(
2015
)
The role of heat shock proteins in cancer
.
Cancer Lett.
360
,
114
118
[PubMed]
35.
Sun
H.
,
Cai
X.
,
Zhou
H.
et al.
(
2018
)
The protein-protein interaction network and clinical significance of heat-shock proteins in esophageal squamous cell carcinoma
.
Amino Acids
50
,
685
697
[PubMed]
36.
Levanon
D.
,
Negreanu
V.
,
Bernstein
Y.
et al.
(
1994
)
AML1, AML2, and AML3, the human members of the runt domain gene-family: cDNA structure, expression, and chromosomal localization
.
Genomics
23
,
425
432
[PubMed]
37.
Wang
Y.
,
Qin
X.
,
Wu
J.
et al.
(
2014
)
Association of promoter methylation of RUNX3 gene with the development of esophageal cancer: a meta analysis
.
PLoS ONE
9
,
e107598
[PubMed]
38.
Sugiura
H.
,
Ishiguro
H.
,
Kuwabara
Y.
et al.
(
2008
)
Decreased expression of RUNX3 is correlated with tumor progression and poor prognosis in patients with esophageal squamous cell carcinoma
.
Oncol. Rep.
19
,
713
719
[PubMed]
39.
Solinas
G.
,
Germano
G.
,
Mantovani
A.
and
Allavena
P.
(
2009
)
Tumor-associated macrophages (TAM) as major players of the cancer-related inflammation
.
J. Leukoc. Biol.
86
,
1065
1073
[PubMed]
40.
Wynn
T.A.
,
Chawla
A.
and
Pollard
J.W.
(
2013
)
Macrophage biology in development, homeostasis and disease
.
Nature
496
,
445
455
[PubMed]
41.
Miao
X.
,
Leng
X.
and
Zhang
Q.
(
2017
)
The current state of nanoparticle-induced macrophage polarization and reprogramming research
.
Int. J. Mol. Sci.
18
,
42.
Li
J.
,
Xie
Y.
,
Wang
X.
et al.
(
2019
)
Prognostic impact of tumor-associated macrophage infiltration in esophageal cancer: a meta-analysis
.
Future Oncol.
15
,
2303
2317
[PubMed]
43.
Sugimura
K.
,
Miyata
H.
,
Tanaka
K.
et al.
(
2015
)
High infiltration of tumor-associated macrophages is associated with a poor response to chemotherapy and poor prognosis of patients undergoing neoadjuvant chemotherapy for esophageal cancer
.
J. Surg. Oncol.
111
,
752
759
[PubMed]
44.
Yagi
T.
,
Baba
Y.
,
Okadome
K.
et al.
(
2019
)
Tumour-associated macrophages are associated with poor prognosis and programmed death ligand 1 expression in oesophageal cancer
.
Eur. J. Cancer
111
,
38
49
[PubMed]
45.
Ikeguchi
M.
,
Ikeda
M.
,
Tatebe
S.
,
Maeta
M.
and
Kaibara
N.
(
1998
)
Clinical significance of dendritic cell infiltration in esophageal squamous cell carcinoma
.
Oncol. Rep.
5
,
1185
1189
[PubMed]
46.
Tran Janco
J.M.
,
Lamichhane
P.
,
Karyampudi
L.
and
Knutson
K.L.
(
2015
)
Tumor-infiltrating dendritic cells in cancer pathogenesis
.
J. Immunol.
194
,
2985
2991
[PubMed]
47.
Vivier
E.
,
Tomasello
E.
,
Baratin
M.
,
Walzer
T.
and
Ugolini
S.
(
2008
)
Functions of natural killer cells
.
Nat. Immunol.
9
,
503
510
[PubMed]
48.
Batista
F.D.
and
Harwood
N.E.
(
2009
)
The who, how and where of antigen presentation to B cells
.
Nat. Rev. Immunol.
9
,
15
27
[PubMed]
49.
Shankaran
V.
,
Ikeda
H.
,
Bruce
A.T.
et al.
(
2001
)
IFNgamma and lymphocytes prevent primary tumour development and shape tumour immunogenicity
.
Nature
410
,
1107
1111
[PubMed]
50.
Leclerc
M.
,
Voilin
E.
,
Gros
G.
et al.
(
2019
)
Regulation of antitumour CD8 T-cell immunity and checkpoint blockade immunotherapy by Neuropilin-1
.
Nat. Commun.
10
,
3345
[PubMed]
51.
Schumacher
K.
,
Haensch
W.
,
Roefzaad
C.
and
Schlag
P.M.
(
2001
)
Prognostic significance of activated CD8(+) T cell infiltrations within esophageal carcinomas
.
Cancer Res.
61
,
3932
3936
[PubMed]
52.
Gao
Y.
,
Guo
W.
,
Geng
X.
et al.
(
2020
)
Prognostic value of tumor-infiltrating lymphocytes in esophageal cancer: an updated meta-analysis of 30 studies with 5,122 patients
.
Ann. Transl. Med.
8
,
822
[PubMed]

Author notes

*

These authors contributed equally to this work.

This is an open access article published by Portland Press Limited on behalf of the Biochemical Society and distributed under the Creative Commons Attribution License 4.0 (CC BY).

Supplementary data