To explore DEGs between EC and normal tissues, gene expression data were retrieved from the esophagus dataset (20160128) from The Cancer Genome Atlas (TCGA) database (http://tcga-data.nci.nih.gov; accessed, November 10, 2018). The ‘RTCGAToolbox’ package (version 1.0) in the R software environment was used to download quantitative gene expression data and clinical characteristics of patients with EC in the TCGA database (6). The running time was ‘20160128’ (updated version of esophageal cancer data in TCGA database when analyzing data). A total of 1,413 datasets of human ECs were retrieved from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). Following reading and screening, the three gene expression datasets GSE23400 (7), GSE20347 (8) and GSE5364 (9) were selected for further analysis. Among these, GSE23400 and GSE5364 were based on the GPL96 platform [(HG-U133A) Affymetrix Human Genome U133A Array; Affymetrix; Thermo Fisher Scientific, Inc.], and GSE20347 was based on the Agilent GPL571 platform [(HG-U133A_2) Affymetrix Human Genome U133A 2.0 Array; Affymetrix; Thermo Fisher Scientific, Inc.].

The ‘RTCGAToolbox’ was used to evaluate DEGs between EC and normal tissue samples using data from the TCGA database (6), and genes that met the cutoff criteria [adjusted P-value =0.05; |log fold change (FC)| ≥2.0] were considered to be significant DEGs.

The GEO2R online analysis tool (https://www.ncbi.nlm.nih.gov/geo/geo2r/) was used to evaluate DEGs between EC and normal tissue samples in data from the GEO database. Due to the small sample size, the parameters used to indicate significant DEGs were adjusted P-value =0.05 and |log FC| ≥1.0.

Finally, a Venn diagram web tool (http://bioinformatics.psb.ugent.be/webtools/Venn/) was used to identify the intersection of DEGs between the TCGA and GEO datasets. The present study first identified the intersection of the three GEO datasets, and then cross-identified the DEGs common to the TCGA and GEO datasets. These intersecting genes were considered to be the DEGs between EC and normal tissue samples.

GO (10) and KEGG (11) pathway enrichment analysis of DEGs was performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID; http://david.ncifcrf.gov/) (12). P<0.05 was considered to indicate a statistically significant difference.

PPI information was analyzed using the STRING database (http://string-db.org/) (13). PPI pairs were extracted based on the condition of medium confidence >0.4. The PPI network was then visualized using Cytoscape software (www.cytoscape.org/; version 22.03.6.1) (14). In the present study, genes with a connectivity >40 were considered to be hub genes. Plug-in molecular complex detection (MCODE) was used to screen PPI networks (15). The maximum depth (value =100), the degree cutoff (value =10), the node score (value =0.2) and the k score (value =2) were set as the cutoff criteria.

To explore the association between hub gene expression and prognosis of EC, the present study used a Kaplan-Meier plotter tool (http://kmplot.com/analysis/) to assess the impact of differential expression of the hub genes on the prognosis of EC (16). The Kaplan-Meier plotter database data are mainly derived from GEO (Affymetrix microarray only), European Genome-phenome Archive (EGA) and TCGA. Overall survival (OS) and relapse-free survival (RFS) were selected as the main indicators of prognostic assessment. The mean expression of all probes of the same gene was calculated as the expression of each gene. All possible cutoff values between the lower and upper quartiles were computed and the best performing threshold was used as the cutoff. The log-rank test was used to analyze OS and RFS. A log-rank P<0.05 was considered to indicate a statistically significant difference.

To explore the role of the hub DEGs in EC invasion, an additional dataset, GSE21293 (17), was evaluated from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). This dataset contained gene expression data for invasive and non-invasive EC cells. Fragments per kilobase million values for ubiquitin conjugating enzyme E2 C (UBE2C), cyclin dependent kinase inhibitor 3 (CDKN3), CDC28 protein kinase regulatory subunit 2 (CKS2), kinesin family member 20A (KIF20A) and RAD51 associated protein 1 (RAD51AP1) genes were log₂ transformed and the resulting heatmap was constructed using HemI software (version 1.0) (18). Graph design and data anlysis were performed using GraphPad Prism 5.0 (GraphPad Software, Inc.). Student's t test was used to analyze difference between two groups. P<0.05 was considered to indicate a statistically significant difference.

To explore the expression levels of hub genes across different stages of EC, the Gene Expression Profiling Interactive Analysis tool (http://gepia.cancer-pku.cn/) was used to analyze the hub gene expression data from the TCGA EC dataset, which contained clinical staging information (19).

To analyze the differential expression of the prognostic hub genes in different histological subtypes of EC, the interactive website Ualcan (http://ualcan.path.uab.edu/analysis.html) was used to analyze the hub gene expression data from the TCGA EC dataset, which contained histological subtype information (20).

To identify DEGs between EC and normal tissues, the present study used two public gene expression databases. Initially, the TCGA database EC data (including 89 EAC samples and 11 normal esophageal samples) were used to identify DEGs. It was observed that samples from the same tissue type aggregated into a single group, and each tissue type exhibited one or more unique tissue-specific mRNA blocks. The top 50 genes were significantly upregulated in EC tissues compared with in normal tissues (Fig. 1A).

Subsequently, gene expression in three microarray datasets [GSE5364 (7), GSE20347 (8) and GSE23400 (9)] from the GEO database was explored. Among these, GSE5364 contained 16 EC and 13 normal samples, GSE20347 contained 17 EC and 17 normal samples, and GSE23400 contained 53 EC and 53 normal samples (Table I). Due to the small sample size, the criteria P<0.05 and |log FC| ≥1 were used to identify DEGs. A total of 1,576 DEGs were identified from the GSE5364 dataset, including 789 upregulated and 787 downregulated genes. In the GSE20347 dataset, 1,368 DEGs were identified; 632 genes were upregulated and 736 genes were downregulated. In the GSE23400 dataset, 522 DEGs were identified, including 255 upregulated genes and 267 downregulated genes. To avoid bias from individual studies, a Venn diagram was constructed to identify genes that were commonly upregulated and downregulated among the three datasets. The present study identified 426 DEGs that were significantly differentially expressed among all three datasets, of which 196 were significantly upregulated (Fig. 1B) and 230 were downregulated (Fig. 1C).

Based on the significance criteria of P<0.05 and |log FC| ≥2, a total of 2,000 DEGs were identified, including 1,156 upregulated genes and 844 downregulated genes in TCGA datasets (Fig. 1D-E). A comprehensive analysis of the data from TCGA and GEO was performed using Venn diagram analysis to identify the intersection of all DEGs. The present study identified 105 genes that were commonly significantly differentially expressed between EC and normal tissues in the two databases, of which 83 were significantly upregulated (Fig. 1D) and 22 were downregulated (Fig. 1E).

To investigate the effects of DEGs on the occurrence of EC at a functional level, the DAVID web tool was used to perform GO functional and KEGG pathway enrichment analysis of the DEGs. The enriched GO terms were divided into cellular components (CCs), biological processes (BPs) and molecular functions (MFs). The BPs associated with upregulated DEGs were ‘collagen catabolic process’, ‘cell division’, ‘extracellular matrix disassembly’, ‘mitotic nuclear division’ and ‘cell proliferation’. For CCs, the upregulated DEGs were mainly enriched in ‘midbody’, ‘spindle’, ‘proteinaceous extracellular matrix’, ‘spindle microtubule’ and ‘kinetochore’. MF analysis demonstrated that the upregulated DEGs were mainly enriched in ‘metalloendopeptidase activity’, ‘protein binding’, ‘ATP binding’, ‘serine-type endopeptidase activity’ and ‘protein kinase binding’. KEGG pathway analysis revealed that the upregulated DEGs were mainly enriched in ‘cell cycle’, ‘DNA replication’, ‘ECM-receptor interaction’, ‘oocyte meiosis’ and ‘progesterone-mediated oocyte maturation’ pathways (Fig. 2A).

The downregulated DEGs were mainly enriched in BPs, including ‘muscle contraction’, ‘complement activation, alternative pathway’ and ‘regulation of smooth muscle contraction’. The downregulated DEGs were significantly enriched in the CC terms ‘focal adhesion’, ‘extracellular exosome’, ‘myosin filament’ and ‘apical plasma membrane’. Enrichment analysis of MF terms revealed that downregulated DEGs were mainly enriched in ‘calmodulin binding’ and ‘structural constituent of muscle’. There was no significant KEGG pathway enrichment for the downregulated DEGs (Fig. 2B).

Cytoscape software was used to predict protein interactions among DEGs. The PPI network comprised 82 nodes and 982 edges (Fig. 3A). After analyzing the degree of connectivity in the PPI network, genes with a connectivity degree >40 were considered to be hub genes. There were 36 hub genes, and all of them were upregulated genes in EC (Table II).

Screening of the PPI network using MCODE identified two modules. In module 1, the MCODE score was 40.7, including 41 nodes and 814 edge connection lines (Fig. 3B). To investigate the effect of modules on the occurrence of EC at a more functional level, DEGs within the modules were classified into KEGG terms. The DEGs in module 1 were mainly enriched in KEGG pathways such as ‘cell cycle’, ‘DNA replication’, ‘oocyte meiosis’, ‘progesterone-mediated oocyte maturation’ and ‘HTLV–I infection’ (P<0.05; Fig. 3C). The MCODE score in module 2 was 8.4, including 16 nodes and 63 edge lines (Fig. 3D), and the KEGG pathway enrichment analysis of module 2 mostly identified terms such as ‘protein digestion and absorption’, ‘ECM-receptor interaction’, ‘focal adhesion’, ‘amoebiasis’, ‘PI3K-Akt signaling pathway’ and ‘platelet activation’ (P<0.05; Fig. 3E).

To determine whether these DEGs were associated with survival of patients with EC, the DEGs were submitted to the Kaplan-Meier plotter bioinformatics analysis platform. Kaplan-Meier plotter database data is mainly derived from GEO (Affymetrix microarray only), EGA and TCGA. All possible cutoff values between the lower and upper quartiles were computed, and the best performing threshold was used as a cutoff. The log-rank test was used to analyze OS (161 EC samples) and RFS (73 EC samples). The OS prognosis of patients with EC with high expression levels of UBE2C (Fig. 4A), CDKN3 (Fig. 4B), CKS2 (Fig. 4C), KIF20A (Fig. 4D) and RAD51AP1 (Fig. 4E) was worse than that of patients with low expression levels of these genes (P<0.05). There was no significant prognostic value identified for the other hub genes (all P>0.05). To further investigate the impact of these five hub genes on the prognosis of patients with EC, the association between the expression of these genes and RFS was analyzed. The results demonstrated that UBE2C expression (Fig. 4F) was not associated with RFS, and higher expression levels of CDKN3 (Fig. 4G) and CKS2 (Fig. 4H) were significantly associated with shorter RFS in patients with EC. However, the expression levels of KIF20A (Fig. 4I) and RAD51AP1 (Fig. 4J) were not associated with RFS in patients with EC.

Poor prognosis of EC is closely associated with early invasion and metastasis (21). To determine if the five prognostic hub genes mediated EC migration and invasion, differences in expression levels of these genes in invading and non-invading EC cells from a GEO microarray dataset (GSE21293) were analyzed (17). Primary esophageal cells were established from ESCC surgical specimens (n=35). The present study revealed that UBE2C, CDKN3, CKS2, KIF20A and RAD51AP1 were significantly differentially expressed between invading and non-invading EC cells (Fig. 5A). UBE2C, CDKN3, CKS2, KIF20A and RAD51AP1 were upregulated in the invading cells (n=12) compared with the non-invading cells using Student's t test (Fig. 5B). These results indicated that these hub genes may mediate the migration and invasion of EC.

Pathological stage is a major indicator of the progression of EC (3). The present study ascertained the association between expression levels of the five prognostic hub genes and pathological staging of EC. There were no significant differences identified for expression of UBE2C (Fig. 6A), CDKN3 (Fig. 6B) or CKS2 (Fig. 6C) in different pathological stages of EC (P>0.05). However, KIF20A expression was associated with different pathological stages of EC, and upregulation of KIF20A was closely associated with advanced EC stage (P<0.05; Fig. 6D). RAD51AP1 expression exhibited a non-significant trend of difference in different pathological stages of EC (P>0.05; Fig. 6E). Therefore, among these five hub genes, KIF20A expression may be used to predict the pathological stage of EC.

To determine if the expression of the five prognostic hub genes differed in different histological subtypes of EC, the expression levels of these genes in ESCC and EAC were analyzed based on the TCGA database. The expression levels of UBE2C (Fig. 7A), CDKN3 (Fig. 7B), CKS2 (Fig. 7C), KIF20A (Fig. 7D) and RAD51AP1 (Fig. 7E) in EAC and ESCC were significantly higher than those in normal tissues. The expression levels of UBE2C, CDKN3 and CKS2 exhibited no significant difference between different histological EC subtypes. Higher expression levels of KIF20A and RAD51AP1 were identified in ESCC compared with EAC. KIF20A and RAD51AP1 may represent more specific molecular targets for ESCC than EAC.