Construction of the Interaction Network of Hub Genes in the Progression of Barrett’s Esophagus to Esophageal Adenocarcinoma

Introduction

Esophageal cancer includes esophageal squamous cell carcinoma (ESCC) and esophageal adenocarcinoma (EAC), of which primary EAC accounts for only 5%-10%. However, the prognosis of EAC is poor, 5-year survival rate of which is less than 20%, and the incidence rate is increasing year by year.¹ EAC is more prevalent in the United States, Australia, and Western European countries. Studies have shown that 2% to 7% of adults in Western countries have Barrett’s esophagus (BE), a major risk factors for EAC, resulting in a more than seven-fold increase in the number of patients with EAC in the United States over the past few decades.²

Studies have confirmed that BE is a precancerous lesion of EAC, and majority of EAC is derived from BE.³ Epidemiological studies have found that the risk of cancerous changes in BE is 25 to 30 times higher than that of normal people.⁴ BE is thought to evolve into EACs through a series of genetic and epigenetic alterations that accumulate in metaplastic cells, giving them an advantage in progressive growth, ultimately leading to malignancy.⁵ Most EAC patients are unaware of the medical history of BE, or do not have timely intervention for BE, resulting in advanced tumor progression when the patients are diagnosed with EAC. Therefore, finding biomarkers for early diagnosis of EAC and high-risk BE patients has important preventive significance and clinical value. Currently, endoscopy-based tissue biopsy is the primary means of assessing whether BE will progress to EAC.⁶ However, numerous reports suggest that patients with BE may develop incurable esophageal cancer despite adherence to a regular endoscopic surveillance program, seriously challenging the efficacy of surveillance endoscopy for cancer prevention.^7–9 The researchers are exploring other approaches that may improve risk prediction, such as the p53 protein and flow cytometric DNA content analysis, both associated with malignant transformation of BE.^10–12 In addition, some researchers have proposed risk prediction models incorporating clinical information, histological or molecular features.^13,14 However, none of these risk prediction models have sufficient evidence to be justified in clinical practice.

Epithelial metaplasia in BE is widely considered to be a consequence of gastroesophageal reflux disease (GERD).¹⁵ BE forms when GERD damages the squamous epithelium of the esophagus, exposing the basal layer of pluripotent stem cells to the refluxed gastric juice and thus stimulating abnormal differentiation. The study found that GERD may induce the expression of CDX gene through BMP4 and epidermal growth factor receptor, causing esophageal squamous cells to express markers of columnar cells, and finally mediate the development of Barrett’s metaplasia.^16–18 Esophageal cells exposed to acid develop DNA double-strand breaks, which may contribute to tumorigenesis.¹⁷ Acid and bile can also resist apoptosis by activating the NF-κB pathway and regulate the expression of inflammation and proliferation-related genes.¹⁹ In Barrett-associated esophageal adenocarcinoma, bile acid activated tuberous sclerosis complex 1/rapamycin pathway.²⁰ Currently, there are few studies focusing on the malignant transformation of BE, and its molecular mechanism is largely unknown.

In this study, we screened the RNA-seq data of 94 NE tissues, 113 BE tissues, and 147 EAC tissues in the GEO database. Then, the hub genes in the development of normal esophageal squamous epithelium (NE)-BE-EAC were explored through bioinformatics combined with immunohistochemistry. Furthermore, we also analyzed the probable mechanisms of these hub genes and their potential as biomarkers for EAC prognosis, risk appraisal of disease progression and clinical treatment. Our study is based on a large number of clinical samples and have high reliability, which is helpful for revealing the probable molecular mechanism of malignant transformation of BE and formulating clinical treatment strategies.

Materials and Methods

Data Collection and Experiment Grouping

In the GEO database, we screened 7 datasets containing RNA-seq data from 94 cases of NE tissues, 113 cases of BE tissues and 147 cases of EAC tissues. In detail, the GSE13083 dataset contains 7 cases of NE tissues and 7 cases of BE tissues; the GSE26886 dataset contains 19 cases of NE tissues, 20 cases of BE tissues and 21 cases of EAC tissues; the GSE42363 dataset contains 14 cases of EAC tissues; GSE13898 dataset contains 28 cases of NE tissues, 15 cases of BE tissues and 75 cases of EAC tissues; GSE39491 dataset contains 40 cases of NE tissues and 40 cases of BE tissues; GSE37201 dataset contains 22 cases of EAC tissues; GSE37200 dataset contains 31 cases of BE tissues and 15 cases of EAC tissues.

Data Preprocessing and Analysis of Differentially Expressed Genes (DEGs)

The seven datasets in this study are all from the GEO database. For the datasets that have not been normalized, log2 processing is performed uniformly. If the dataset data is not standardized, the normalize.quantiles function in the R software preprocessCore package is used to perform data normalization. According to the annotation information of the corresponding platform, statistics convert the probe ID to gene symbol, remove the probes of multiple genes, and calculate the average value of the genes corresponding to multiple probes. The removeBatchEffect function in the R software limma package was used to remove batch effects. Differentially expressed mRNA was investigated using the limma package of R software. “P < 0.05 and log2(fold change) > 0.58 or log2(fold change) < −0.58” was defined as the threshold for the screening of DEGs. Expression heatmaps are displayed by the R package pheatmap. To better understand the probable mechanism of DEGs, the ClusterProfiler package in R software was used to perform GO analysis and KEGG analysis of DEGs.

Screening of Hub Genes

The overlapping DEGs between NE and BE, NE and EAC, BE and EAC were analyzed based on the Venny online tool, and then the protein-protein interaction network of these overlapping DEGs was constructed using the STRING website. The original data of the protein interaction network was imported into the Cytoscape software, and the cytohubba plug-in was used to screen the top 7 as hub genes according to degree. Then, based on the 7 datasets of the GEO database, we utilized boxplot to visualize the expression levels of hub genes in NE tissue, BE tissue and EAC tissue. A p-value less than 0.05 was considered statistically significant.

Correlation of Hub Genes with Immune Cell Infiltration, Immunotherapy Response, Immune Checkpoint Genes, Tumor Mutational Burden (TMB)

The RNA-seq data (level 3) and corresponding clinical information of EAC were obtained from the TCGA database (https://portal.gdc.com). We used the R software package immunedeconv, including four algorithms TIMER, MCP-counter, EPIC and quanTIseq, to analyze the correlation between hub genes and different immune cell infiltration, and the analysis results were displayed by the R software package pheatmap. Spearman correlation analysis was used to describe correlations between quantitative variables without a normal distribution. A p-value less than 0.05 was considered statistically significant.

The TIDE algorithm uses a set of gene expression signatures to assess two distinct tumor immune escape mechanisms, including tumor-infiltrating cytotoxic T lymphocyte (CTL) dysfunction and CTL rejection by immunosuppressive factors.²¹ High TIDE score, poor response to immune checkpoint blockade (ICB), and short survival after receiving ICB. Therefore, immunotherapy responses of different groups with high or low expression of hub genes can be predicted using the TIDE algorithm. The results are displayed by the R packages ggplot2 and ggpubr. A p-value less than 0.05 was considered statistically significant.

Based on the 7 datasets in the GEO database, we extracted the expression values of immune checkpoint genes CD274, CTLA4, HAVCR2, LAG3, PDCD1, PDCD1LG2, TIGIT and SIGLEC15, and then used the R software package ggstatsplot to analyze the relationship between hub genes and immune checkpoint genes. A p-value less than 0.05 was considered statistically significant.

Based on EAC RNA-seq data (level 3) in the TCGA dataset (https://portal.gdc.com), we used Spearman to analyze the correlation between hub genes and TMB. The results are presented by the R package ggstatsplot. A p-value less than 0.05 was considered statistically significant.

Construction of Hub Genes ceRNA Network

Using the ENCORI database (https://starbase.sysu.edu.cn/index.php), we predicted miRNAs that could target both hub genes and their ceRNA molecules. Then, the list of candidate hub genes, lncRNAs and miRNAs was imported into Cytoscape software to construct a ceRNA network of hub genes.

Correlation of Hub Genes and EAC Prognosis

The online tool Kaplan-Meier Plotter (http://kmplot.com/analysis/index.php?p=serviceandcancer=pancancer_rnaseq) was used to analyze the correlation between hub genes and overall survival (OS) and recurrence-free survival (RFS) of EAC patients. The EAC samples were divided into two groups based on the expression of hub genes by the best cut-off value automatically and were analyzed using Kaplan-Meier method and Log rank-P test (P<0.05 as significant).

Correlation Between Hub Genes and Signaling Pathways

The RNA-seq data (level 3) and corresponding clinical information of EAC were obtained from the TCGA database (https://portal.gdc.com). The key genes contained in the corresponding pathways were collected and analyzed by the R software GSVA package. The parameter method=“ssgsea” was selected, and finally the correlation between genes and pathway scores was analyzed by Spearman correlation. p<0.05 was considered statistically significant.

Construction of Predictive Model for Cancer Diagnosis

Firstly, the cancer diagnosis performance of 7 hub genes expression was evaluated by using receiver operating characteristic (ROC) curves, and the overall diagnostic accuracy of the test was assessed by the area under the ROC curve (AUC). Then the logistic regression test was used to identify the genes that were significantly associated with cancer diagnosis based on the 7 hub genes. According to the output results of the 7 genes predictive model, the final logistic regression model for cancer diagnostic was constructed by TGFBI, MMP1 and CXCL1 showing the most significant P values and the exponentiated values of the coefficients, and the ROC curve was utilized to judge discrimination of the 3 genes predictive model. All analyses were performed using SPSS.

Drug Sensitivity Analysis

Using the “drug sensitivity analysis” module on GSCALite database (http://bioinfo.life.hust.edu.cn/web/GSCALite/), we analyzed the correlation between hub genes and 265 small molecules from Genomics of Drug Sensitivity in Cancer (GDSC). The correlation of hub genes expression with the small molecule/drug sensitivity (IC50) was performed by Spearman correlation analysis. The negative correlation means that the gene high expression is sensitive to the drug, vise verse. p<0.05 was considered statistically significant.

Immunohistochemistry of NE Tissue, BE Tissue and EAC Tissue

We verified the protein expression levels of hub genes in NE (n=20), BE (n=20) and EAC (n=20) tissues by immunohistochemistry. These tissues came from Nanyang Central Hospital and were diagnosed by the pathology department. The experimental steps of immunohistochemistry refer to our previous article.²² To compare the protein expression levels of hub genes among NE, BE and EAC tissues, the intensity of immunohistochemistry staining was then quantitatively scored. The criteria for quantitative scoring refer to our previous article.²³

Results

Exploration of the DEGs Between NE Tissue and BE Tissue, NE Tissue and EAC Tissue, BE Tissue and EAC Tissue, and their Functional Enrichment Analysis

Based on the GEO database (http://www.ncbi.nih.gov/geo), we analyzed RNA-seq data of 94 NE tissues, 113 BE tissues, and 147 EAC tissues. There were 1173 upregulated genes and 1051 downregulated genes in BE tissues relative to NE tissues (Figure 1A, Table S1). The functions of these genes in the development of BE have not yet been revealed. Therefore, we performed KEGG and GO analyses on these DEGs, to preliminarily reveal their probable mechanisms in the occurrence and development of BE. The results of functional enrichment analysis showed that these DEGs may be related to immune response, gastric acid secretion, protein digestion and absorption, adhesion junctions, tight junctions and TGF-β signaling pathway, mTOR signaling pathway, MAPK signaling pathways, etc. (Figure 1B, Table S2). Relative to the NE tissues, there were 970 upregulated genes and 944 downregulated genes in the EAC tissues (Figure 1C, Table S3). Similarly, we performed functional enrichment analysis on these DEGs, and the results showed that their functions were similar to the DEGs in BE, such as immune response, gastric acid secretion, adhesion junctions, tight junctions, TGF-β signaling pathway, MAPK signaling pathway, etc. (Figure 1D, Table S4). By comparing RNA-seq data of BE tissue and EAC tissue, we found 94 upregulated genes and 296 downregulated genes in EAC relative to BE (Figure 1E, Table S5). Then, we performed functional enrichment analysis of these DEGs. The results show that these DEGs are closely related to inflammatory processes such as immune response, TLR signaling pathway, TNF signaling pathway, and NF-κB signaling pathway (Figure 1F, Table S6).

Figure 1 Analysis of DEGs during the development of NE-BE-EAC. (A) The limma package of R software was used to identify DEGs between NE and BE based on RNA-seq data of 94 NE tissues and 113 BE tissues from the GEO database. (B) The ClusterProfiler package in R software was used to perform GO analysis and KEGG analysis of up- and down-regulated DEGs. Similarly, DEGs between NE and EAC were analyzed based on RNA-seq data of 94 NE tissues and 147 EAC tissues in the GEO database (C), And the GO and KEGG analysis were performed using ClusterProfiler package in R software (D). Based on the RNA-seq data of 113 BE tissues and 147 EAC tissues in the GEO database, the DEGs between BE and EAC were analyzed (E), and the GO and KEGG analysis were performed using ClusterProfiler package in R software (F).

7 Hub Genes Were Identified During the Transition from BE to EAC

To explore oncogenes and tumor suppressor genes that may play a role in the development of NE-BE-EAC, we firstly identified the overlapping DEGs between BE (1173 upregulated genes and 1051 downregulated genes) and EAC (970 upregulated genes and 944 downregulated genes) by Wayne analysis. The results showed that the DEGs in BE and EAC had great similarity, sharing 813 overlapping upregulated genes and 841 overlapping downregulated genes (Figure 2A and B, Table S7). This result further proves the close connection between BE and EAC. In this study, we aim to explore the genes whose expression levels were progressively upregulated or downregulated during the development of NE-BE-EAC. Therefore, we proceeded to compare the above overlapping DEGs with those in EAC (relative to BE) by Wayne analysis. Through analysis, we identified 12 upregulated genes (COL1A1, IFI30, AEBP1, AGT, CXCL1, NID2, COL4A1, UBD, SERPINA3, TGFBI, MMP1, MMP12) and 7 downregulated genes (MAL, DSC2, S100A14, MALL, DHRS9, FGFBP1, ALDH3A1), whose expression levels were progressively upregulated or downregulated during the development of NE-BE-EAC (Figure 2C–E). To further identify the hub genes involved in the conversion of BE to EAC, we constructed the protein-protein interaction network of these 19 genes base on the STRING database (Figure 2F). Then based on this network, we finally screened 7 hub genes using Cytoscape software, including COL1A1, TGFBI, MMP1, COL4A1, NID2, MMP12, CXCL1 (Figure 2G). Finally, we further verified the expression levels of these genes in 94 NE tissues, 113 BE tissues and 147 EAC tissues, and the results showed that relative to NE tissues, COL1A1, IFI30, AEBP1, AGT, CXCL1, NID2, COL4A1, UBD, SERPINA3, TGFBI, MMP1, MMP12 were progressively upregulated in BE and EAC, while MAL, DSC2, S100A14, MALL, DHRS9, FGFBP1, ALDH3A1 were progressively downregulated (Figure 3A–G, Figure S1A–L). These results suggest that COL1A1, TGFBI, MMP1, COL4A1, NID2, MMP12, CXCL1 may act as key oncogenes in the development of NE-BE-EAC.

Figure 2 Exploration of hub genes during the development of NE-BE-EAC. (A and B) Relative to NE, up- and down-regulated DEGs overlapping in BE and EAC were identified by Wayne analysis. (C and D) Compared to BE, the upregulated and downregulated genes in EAC are labeled as “URGs in EAC-BE” and “DRGs in EAC-BE”, respectively. Relative to NE, overlapping up- and down-regulated genes in BE and EAC were labeled as “URGs in EAC+BE” and “DRGs in EAC+BE”, respectively. Overlapping genes between “URGs in EAC-BE” and “URGs in EAC+BE”, “DRGs in EAC-BE” and “DRGs in EAC+BE” were identified by Wayne analysis. (E) These overlapping DEGs include 12 upregulated genes and 7 downregulated genes. (F) Analysis of protein interactions between these 19 DEGs using the STRING database. (G) Finally, 7 hub genes were screened based on protein interactions using the cytohubba plugin in Cytoscape software. The linkages between the proteins indicate the predicted interactions.

Figure 3 Analysis of mRNA expression levels of hub genes in NE, BE and EAC. Based on the RNA-seq data in the GEO database, the mRNA expression levels of COL1A1 (A), TGFBI (B), MMP1 (C), COL4A1 (D), NID2 (E), MMP12 (F), CXCL1 (G) in 94 NE tissues, 113 BE tissues and 147 EAC tissues were analyzed. ***(P < 0.001) or ****(P < 0.0001) represents statistical significance.

The Protein Expression Levels of Hub Genes are Progressively Upregulated in the Development of NE-BE-EAC

Based on RNA-seq data, we have determined the mRNA expression levels of hub genes in NE, BE and EAC tissues. Furthermore, NE, BE and EAC tissues were collected to verify the protein expression levels of hub genes by immunohistochemistry. According to the quantitative scores of immunohistochemistry, we confirmed that the expression levels of COL1A1, TGFBI, MMP1, COL4A1, NID2, MMP12, CXCL1 in BE tissues were significantly higher than that in NE tissues, but lower than that in EAC tissues, which was consistent with RNA-seq analysis (Figure 4). These results further suggested that these hub genes may play important roles in the development of NE-BE-EAC.

Figure 4 The protein expression levels of hub genes in NE, BE and EAC. The protein expression level of COL1A1, TGFBI, MMP1, COL4A1, NID2, MMP12, CXCL1 was detected and quantified by immunohistochemistry in NE (n=20), BE (n=20) and EAC (n=20) tissues. **(P < 0.01) or ***(P < 0.001) represents statistical significance.

The Hub Genes Were Closely Associated with Immune Cell Infiltration and Tumor Immunotherapy

We used TIMER, MCP-counter, EPIC and quanTIseq to analyze the correlation between hub genes and immune cell infiltration (Figure 5A–D). The results of the four algorithms are mostly consistent, such as 7 oncogenic hub genes were positively correlated with the infiltration of B cells, endothelial cells, macrophages, monocytes, NK cells, neutrophils, Treg cells and other inflammatory cells. In some individual cases, the results of the four algorithms are inconsistent. For example, the EPIC algorithm shows that most of the oncogenic hub genes are negatively correlated with the infiltration of CD8⁺ T cells, while the results of other three algorithms are opposite or not statistically significant. In this study, we also analyzed the correlation between hub genes and immune checkpoint molecules (CD274, CTLA4, HAVCR2, LAG3, PDCD1, PDCD1LG2, TIGIT, SIGLEC15) (Figure 6A–G). The results showed that these hub genes, especially NID2, were positively correlated with the expression of immune checkpoint molecules. The Tumor Immune Dysfunction and Exclusion (TIDE) score can be used to predict patient response to ICB therapy.²¹ A high TIDE score predicts poor ICB efficacy and short survival after ICB treatment. In this study, we found that the expression levels of COL1A1 and NID2 were positively correlated with the TIDE score, suggesting that these three hub genes can be used as biomarkers for predicting response to ICB therapy (Figure S2A–G). In addition to the TIDE score, tumor mutational burden (TMB) is also one of the important markers for predicting the response to ICB therapy.²⁴ Higher TMB indicates that patients may be more sensitive to ICB treatment. The results showed that COL4A1 was negatively correlated with TMB score (Figure S3A–G). Therefore, COL1A1, NID2 and COL4A1 may also serve as biomarkers for predicting the response to ICB therapy.

Figure 5 Correlation analysis of hub genes with immune cell infiltration. Based on the RNA-seq data of EAC in TCGA database, we used the R package immunedeconv, including four latest algorithms TIMER (A), MCP-counter (B), EPIC (C) and quanTIseq (D), to analyze the correlation between hub genes and different immune cell infiltrations, and the analysis results were presented by the R package pheatmap. Spearman correlation analysis was used to characterize the correlation between hub genes and different types of immune cells. *(P < 0.05) or **(P < 0.01) represents statistical significance.

Figure 6 Correlation analysis of hub genes with immune checkpoint molecules. Based on the RNA-seq data of NE, BE and EAC in the GEO database, we extracted the expression values of immune checkpoint genes CD274, CTLA4, HAVCR2, LAG3, PDCD1, PDCD1LG2, TIGIT and SIGLEC15, and then analyzed the correlation between COL1A1 (A), TGFBI (B), MMP1 (C), COL4A1 (D), NID2 (E), MMP12 (F), CXCL1 (G) and immune checkpoint genes using the R package ggstatsplot. *(P < 0.05), **(P < 0.01) or ***(P < 0.001) represents statistical significance.

Construction of the Downstream Regulatory Network of Hub Genes

The molecular mechanisms of these hub genes in the development of EAC are largely unknown. In this study, we divided the EAC tissues in the TCGA database into two groups based on the mean expression level of hub gene, then screened out the DEGs between these two groups, and performed functional enrichment analysis (Figure 7A–G, Table S8–14). The analysis showed that there is a great degree of similarity in the functions of these hub genes in EAC, for example, they are all associated with intercellular junctions, cytoskeleton regulation, migration of inflammatory cells, epithelial cell proliferation, lymphocyte activation, vascular smooth muscle contraction, etc. The probable molecular mechanisms are mainly related to PI3K/Akt signaling pathway, JAK/STAT signaling pathway, NF-κB signaling pathway, TGF-β signaling pathway, TNF signaling pathway, Wnt signaling pathway, TLR signaling pathway, etc. These mechanisms have been shown to be closely related to the development of inflammation and tumorigenesis. In order to further reveal the probable mechanism of hub genes in the development of EAC, based on the RNA-seq data of EAC in the TCGA database, we analyzed the correlation of hub genes with tumor inflammation signature (TIS), cellular response to hypoxia, tumor proliferation signature, EMT markers, ECM-related genes, angiogenesis, apoptosis, DNA repair, G2M checkpoint, inflammatory response, PI3K/AKT/mTOR pathway, P53 pathway, MYC targets, TGF-β, IL-10 anti-inflammatory signaling pathway, genes upregulated by reactive oxygen species (ROS), DNA replication, collagen formation, degradation of ECM, ferroptosis. The TIS score was positively correlated with the sensitivity of ICB therapy. The analysis showed that the expression level of COL1A1 was positively correlated with the TIS score. In addition, COL1A1 was positively correlated with EMT, ECM-related genes, angiogenesis, apoptosis, inflammatory response, TGF-β signaling pathway, IL-10 anti-inflammatory signaling pathway, collagen formation, ECM degradation, and ferroptosis (Figure S4). TGFBI was positively correlated with EMT, ECM-related genes, angiogenesis, apoptosis, inflammatory response, PI3K/AKT/mTOR pathway, TGF-β signaling pathway, IL-10 anti-inflammatory signaling pathway, collagen formation, ECM degradation, and ferritin degradation (Figure S5). MMP1 was positively correlated with EMT, angiogenesis, inflammatory response, IL-10 anti-inflammatory signaling pathway, collagen formation, and ECM degradation (Figure S6). COL4A1 was positively correlated with TIS score, EMT, ECM-related genes, angiogenesis, apoptosis, inflammatory response, TGF-β signaling pathway, IL-10 anti-inflammatory signaling pathway, ROS pathway, collagen formation, ECM degradation, and ferroptosis (Figure S7). NID2 was positively correlated with TIS score, EMT, ECM-related genes, angiogenesis, apoptosis, inflammatory response, PI3K/AKT/mTOR pathway, TGF-β signaling pathway, IL-10 anti-inflammatory signaling pathway, ROS pathway, collagen formation, ECM degradation, and ferroptosis (Figure S8). MMP12 was positively correlated with TIS score, hypoxic response, angiogenesis, apoptosis, inflammatory response, IL-10 anti-inflammatory signaling pathway, degradation of ECM, and ferritin deposition (Figure S9). CXCL1 was positively correlated with angiogenesis, apoptosis, inflammatory response, IL-10 anti-inflammatory signaling pathway, ROS pathway, and ferroptosis (Figure S10). By constructing downstream regulatory network hub genes, we further validated that COL1A1, TGFBI, MMP1, COL4A1, NID2, MMP12, CXCL1 may play important oncogenic roles during NE-BE-EAC development.

Figure 7 Construction of a downstream regulatory network for hub genes in EAC. Based on the RNA-seq data of EAC in TCGA database, EAC was divided into high and low expression groups by the mean of COL1A1 (A), TGFBI (B), MMP1 (C), COL4A1 (D), NID2 (E), MMP12 (F), CXCL1 (G) expression levels, and then screened for DEGs between these two groups. The ClusterProfiler package in R software was used to perform GO analysis and KEGG analysis of up- and down-regulated DEGs.

Construction of Upstream Regulatory Network of Hub Genes Based on ceRNA Mechanism

In the present study, we have confirmed that the mRNA expression levels of hub genes are progressively up- or down-regulated during the development of NE-BE-EAC. However, the molecular mechanisms regulating the expression of these hub genes are not known. In this study, based on the ENCORI database, we constructed the upstream regulatory network of hub genes (Figure 8A–F). For example, in the ceRNA regulatory network of COL1A1, miR-128-3p is able to target both COL1A1 and multiple lncRNAs (OLMALINC, MEG8, NEAT1, AC118758.3), then these lncRNAs may be the ceRNAs of COL1A1. Similarly, we found that KRT18P13 may be a ceRNA for TGFBI. ATP2B1 and SP2-AS1 may be ceRNAs of MMP1. AC009487.2, AC067904.2, MSC-AS1 and AC022306.2 may be the ceRNAs of COL4A1. NEAT1, SNHG8, MALAT1 and other genes may be ceRNAs of NID2. MTND6P4, NEAT1, TUBAP2 and other genes may be ceRNAs of CXCL1.

Figure 8 Construction of the upstream regulatory network of hub genes. Base on the ENCORI database, we predicted miRNAs that simultaneously target COL1A1 (A), TGFBI (B), MMP1 (C), COL4A1 (D), NID2 (E), CXCL1 (F) and their ceRNA molecules. Then, the list of these candidate hub genes, lncRNAs and miRNAs were imported into Cytoscape software to construct the ceRNA network of hub genes. The red color represents hub genes, the pink color represents miRNAs, and the fruit green color represents lncRNAs. The connecting lines between different genes represent the predicted interactions.

The Expression of TGFBI Was Negatively Correlated with the OS of EAC

In the present study, our analysis revealed that hub genes are progressively up- or down-regulated during the development of NE-BE-EAC and may be involved in the development of the disease through multiple mechanisms. However, the correlation between these hub genes and the prognosis of EAC is unclear. We analyzed the correlation of hub genes with OS and RFS of EAC using Kaplan Meier plotter (Figure 9A–N). The analysis showed that the expression level of TGFBI was negatively correlated with the OS of EAC patients. In addition, although the correlation between COL1A1, COL4A1, NID2, MMP12 and overall survival of EAC patients was not statistically significant, there was a trend that the higher the expression levels of COL1A1, COL4A1, NID2, MMP12, the shorter the overall survival of patients. Similarly, we found a trend that the higher the expression levels of COL1A1, TGFBI, COL4A1, NID2, and MMP12, the shorter the RFS of the patients. In conclusion, our findings suggested that TGFBI could be used as biomarkers to predict the prognosis of EAC patients.

Figure 9 Correlation analysis of hub genes with overall survival (OS) and relapse-free survival (RFS) of EAC patients. Based on the RNA-seq data of EAC tissues in the TCGA database, the online tool Kaplan-Meier Plotter was used to analyze the correlation between hub genes and OS (A–G) and RFS (H–N) of EAC patients.

3 Genes Model Was Constructed to Predict the Risk of Disease Development in NE-BE-EAC

The development of NE-BE-EAC is a long process.²⁵ Although, BE is a precancerous lesion of EAC, not all BE will progress to EAC. Then, finding biomarkers that can predict whether BE will become cancerous or not has a very important diagnostic value. Firstly, we analyzed the predictive power of the single hub gene based on Area Under the Curve (AUC). In general, the closer the AUC is to 1, the better the predictive power of the model. According to the ROC curves analysis, all 7 hub genes had a comparatively good diagnostic value, with all of the genes acting as risk factors for the diagnosis of cancer (Figure 10A). Furthermore, we performed logistic regression with all 7 hub genes and the results showed that only TGFBI, MMP1 and CXCL1 could be significantly associated with cancer diagnosis and risk in this model based on the selection of P values (P < 0.1), and the top exponentiated values of the coefficients, although this 7 genes model exhibited excellent diagnostic performance (AUC=0.794) (Figure 10B, Table 1). In this case, the final 3 genes model was constructed with TGFBI, MMP1 and CXCL1 expression data, indicating a satisfactory diagnostic value (AUC=0.787) (Figure 10C, Table 2). Therefore, we believe that the 3 genes model is an optimal model to predict the risk of disease development in NE-BE-EAC.

Table 1 Variables in the Equation Based on 7 Genes Model

Table 2 Variables in the Equation Based on 3 Genes Model

Figure 10 Construction of disease risk prediction model for NE-BE-EAC based on hub genes. Based on RNA-seq data from 94 NE tissues, 113 BE tissues and 147 EAC tissues in the GEO database, the disease risk prediction based on single hub gene was assessed by ROC curves, and the accuracy of the model was assessed by AUC (A). Similarly, we constructed the risk prediction models of disease progression based on 7 genes (B) and 3 genes (TGFBI, MMP1 and CXCL1) (C). The AUC values are presented separately in the ROC curves.

Potential Drugs Were Identified to Inhibit the Development of NE-BE-EAC Based on the Expression of Hub Genes

In the present study, our analysis showed that these hub genes play an important role in the development of NE-BE-EAC, which suggests that they may be promising therapeutic targets. We analyzed the correlation between IC50 drug data of 265 small molecules from Genomics of Drug Sensitivity in Cancer (GDSC) and the expression levels of hub genes, where a positive correlation indicates that a certain drug treatment leads to an upregulation of the expression levels of hub genes (Figure S11). Therefore, through drug sensitivity analysis, we can select drugs that are negatively correlated with the expression levels of hub genes, and these drugs will inhibit the expression levels of hub genes to exert anti-tumor therapeutic effects. Our analysis showed that PI3K inhibitor TGX221, bleomycin, PKC inhibitor Midostaurin, Bcr-Abl inhibitor Dasatinib, HSP90 inhibitor 17-AAG, and Docetaxel were negatively correlated with the expression levels of oncogenic hub genes, suggesting the potential of these drugs to inhibit disease progression in NE-BE-EAC.

Discussion

There is growing evidence that the progression of BE to EAC is a multistep process that includes the development of non-developmental abnormal BE to low-grade atypical hyperplasia, then to high-grade atypical hyperplasia, and ultimately to EAC.^25–27 However, the molecular mechanism by which BE becomes cancerous is largely unknown. In view of the high lethality of EAC, it is of great clinical value to study the carcinogenic mechanisms of BE and search for biomarkers that predict disease progression. From the perspective of tumor prevention, timely detection and treatment of BE is an important strategy to prevent EAC.

BE is the only recognized precancerous lesion of EAC, and almost all EAC originates from BE.²⁸ However, majority of EAC patients are unaware of their history of BE. In this study, we screened RNA-seq data of 94 NE epithelial tissues, 113 BE tissues and 147 EAC tissues in the GEO database, and analyzed the DEGs between normal and BE, normal and EAC, BE and EAC by bioinformatics methods, respectively. Surprisingly, we found that 69% of upregulated genes in BE overlapped with 84% of upregulated genes in EAC relative to NE tissue, while 80% of downregulated genes in BE overlapped with 89% of downregulated genes in EAC. Such a high overlap rate is surprising, suggesting that BE and EAC share highly similar genetic characteristics and coincidentally demonstrating the close association between BE and EAC. It is well known that activation of oncogenes and/or inactivation of oncogenes is the central mechanism leading to tumor development. In particular, during tumorigenesis, the generation of the first cancer cell is a progressive process involving a multi-stage response and accumulation of mutations, requiring the involvement of many oncogenes and tumor suppressor genes.²⁹ Precancerous lesions are an intermediate stage in the transformation of normal cells to cancer cells, which have similar and different genetic characteristics from tumors.^30,31 Moreover, they differ in their biological behavior and clinical manifestations, such as cell proliferation ability, motility, etc., which may be associated with progressive abnormal expression of oncogenes and tumor suppressor genes. In the present study, we identified progressive upregulated oncogenes and downregulated oncogenes in the disease process, which may be key factors influencing the process of NE-BE-EAC.

Dysregulation of immune function is one of the important factors that promote tumor development.³² Tumor cells evade recognition and attack by the immune system through various mechanisms, among which immunosuppression induced by immune checkpoint molecules such as PD-L1 and CTLA is one of the main mechanisms leading to tumor immune escape.^33,34 Under normal conditions, immune checkpoint molecules can maintain immune tolerance by regulating the strength of autoimmune response. However, during tumor development, activation of immune checkpoint PD-L1 can suppress immune function and help tumor cells evade the defense and attack of immune system. In the present study, we found that oncogenic hub genes positively correlated with the expression levels of PD-L1, CTLA and other immune checkpoint molecules, suggesting that oncogenic hub genes may inhibit the surveillance function of the immune system by regulating the expression of immune checkpoint molecules during the transformation of BE to EAC, which ultimately promotes the formation and rapid proliferation of EAC cells. In addition, we found that the expression levels of oncogenic hub genes positively correlated with the infiltration of inflammatory cells such as B cells, endothelial cells, macrophages, monocytes, NK cells, neutrophils, and Treg cells in the EAC tumor microenvironment. Under normal circumstances, B cells, NK cells and other immune cells are involved in the clearance of tumor cells by the immune system, but this requires a sufficient number of immune cells that are functioning properly.^35,36 Upregulated immune checkpoint molecules inhibit the activation of immune cells, so oncogenic hub genes may lead to inactivation of immune cells such as B cells, NK cells, and ultimately promote the proliferation of EAC cells.

Most people with EAC are unaware of whether they have BE or not. In addition, not all BE progresses to cancer, so many people may miss the best time to treat BE because they ignore it. These discouraging results are mainly due to the lack of biomarkers to predict whether BE is cancerous or not. In the present study, we evaluated the potential of hub genes as biomarkers for EAC prognosis, diagnosis and treatment. First, our results suggest that TGFBI can be used as biomarkers to predict the prognosis of EAC patients. In addition, although the correlation of other hub genes with OS and RFS of EAC was not statistically significant, there was some degree of correlation, which may be due to the insufficient sample size of EAC in the TCGA database. It is essential to increase the sample size to recalculate the correlation between hub genes and EAC prognosis. Next, we constructed a risk model for predicting disease progression of NE-BE-EAC, which is one of the core results of this study. By comparing multiple risk models, we finally determined that a 3 genes model based on CXCL1, MMP1 and TGFBI was optimal. According to this model, we can predict the possibility of BE to be cancerous by detecting the expression levels of CXCL1 and TGFBI in tissues. Specifically, when CXCL1, MMP1 and TGFBI are progressively upregulated, it suggests a high risk of disease progression. For the exploration of hub genes as biomarkers for clinical treatment, we found that COL1A1, NID2 and COL4A1 were significantly correlated with TMB, TIDE scores, suggesting that these hub genes can be used as biomarkers to predict the effect of tumor immunotherapy. Finally, based on drug sensitivity analysis, we screened PI3K inhibitor TGX221, bleomycin, PKC inhibitor Midostaurin, Bcr-Abl inhibitor Dasatinib, HSP90 inhibitor 17-AAG, and Docetaxel as potential candidates for inhibiting the progression of BE to EAC.

Currently, there are few studies focusing on the dynamic process of BE conversion to EAC. In the present study, based on RNA-seq data from a large number of clinical samples and experimental validation, we identified hub genes in the disease progression of NE-BE-EAC. Through bioinformatics analysis, we preliminarily revealed the probable molecular mechanisms of these hub genes and constructed the upstream regulatory network of hub genes. More importantly, we examined the potential of hub genes as biomarkers for EAC prognosis, prediction of disease progression and tumor immunotherapy, and proposed potential chemotherapeutic agents to inhibit disease progression in NE-BE-EAC, which have potential clinical applications.