Bioinformatics Characterization of Candidate Genes Associated with Gene Network and miRNA Regulation in Esophageal Squamous Cell Carcinoma Patients

Featured Application

Esophageal squamous cell carcinoma (ESCC) is the most common and frequent epithelial malignancy of the esophageal diagnosed among esophageal carcinoma (ESCA) worldwide. Despite enormous studies, the effect of gene overexpression and the overall survival rate associated with a genes mutation is still poorly understood. A biomarker to predict survival is a critical need in ESCC treatment. We hypothesized that a top ten hub gene (HGsT10) could be used as a predictive biomarker for the treatment response in metastatic ESCC. We found that HGsT10 expressed highly at the tumor state, reflecting the tumor development and poor survival rate in ESCC.

Abstract

The present study aimed to identify potential therapeutic targets for esophageal squamous cell carcinoma (ESCC). The gene expression profile GSE161533 contained 84 samples, in that 28 tumor tissues and 28 normal tissues encoded as ESCC patients were retrieved from the Gene Expression Omnibus database. The obtained data were validated and screened for differentially expressed genes (DEGs) between normal and tumor tissues with the GEO2R tool. Next, the protein–protein network (PPI) was constructed using the (STRING 2.0) and reconstructed with Cytoscape 3.8.2, and the top ten hub genes (HGsT10) were predicted using the Maximal Clique Centrality (MCC) algorithm of the CytoHubba plugin. The identified hub genes were mapped in GSE161533, and their expression was determined and compared with The Cancer Genome Atlas (TCGA.) ESCC patient’s samples. The overall survival rate for HGsT10 wild and mutated types was analyzed with the Gene Expression Profiling Interactive Analysis2 (GEPIA2) server and UCSC Xena database. The functional and pathway enrichment analysis was performed using the WebGestalt database with the reference gene from lumina human ref 8.v3.0 version. The promoter methylation for the HGsT10 was identified using the UALCAN server. Additionally, the miRNA-HGsT10 regulatory network was constructed to identify the top ten hub miRNAs (miRT10). Finally, we identified the top ten novel driving genes from the DEGs of GSE161533 ESCC patient’s sample using a multi-omics approach. It may provide new insights into the diagnosis and treatment for the ESCC affected patients early in the future.

1. Introduction

Esophageal carcinoma (ESCA) is the eighth most frequent cancer globally in terms of incidence and fatality. ESCA pathogenesis is complicated by several genes and proteins with aberrant expression and activity, although the precise cause is unknown [1,2,3]. The etiology of this condition remains unclear, and there is no definitive serum marker [1,4]. It has been recognized that ESCA has various histological subtypes, the most common of which are esophageal adenocarcinoma (AC), esophageal squamous cell carcinoma (ESCC), and small cell carcinoma (SCC) [5]. The squamous cell carcinoma is the most frequent epithelial malignancy of the esophageal [6,7]. Despite continual advancements in the use of surgeries, radiation, chemotherapy, targeted therapy, immunotherapy, and anti-cancer drugs in the holistic treatment of cancer, the 5-year mortality rate is less than 10% globally [8,9,10,11]. The molecular mechanism underlying the incidence and progression of ESCC has yet to be fully understood [12]. Therefore, it is urgent to explore novel biomarkers for diagnosis, prognosis, and development of new anti-cancer drugs for ESCC susceptible individuals. Although, deep processing of publicly accessible genomic data to find new and unique gene predictive signals can help patients with prognostic classification and tailored treatment. Presently, one of the primary goals of ESCC research is to find essential genes or particular biomarkers. It influences the onset and progression of ESCC to identify genetic and epigenetic risk factors and better understand their possible underlying molecular mechanism [12,13,14].

Epigenetics is the study to describe any heritable change in gene expression that is not associated with a difference in the chromosomal DNA sequence [15,16,17]. The most prevalent epigenetic modification is DNA methylation, linked to gene expression regulation, fetal development, specific gene inactivation, genomic integrity, genomic imprinting, and cell proliferation [18,19,20]. DNA methylation of genes at the promoter region is linked to the silence of essential genes, including oncogenes and tumor suppressors, and is considered a cornerstone of many malignancies [21,22,23]. Even though several researchers have demonstrated that specific genes with aberrant DNA hyper or hypomethylation are involved in various stages of ESCC development, determining the complete profile and mechanisms of the interaction network remains problematic [24,25,26]. In this pertain, we have predicted methylation at the promoter region for the top ten hub genes (HGsT10).

The ability to accurately identify and characterize the metabolic pathway from the identified hub genes that regulate the biological processes related to ESCC is essential. It may be provided by implementing an extensive assessment of numerous datasets. This perspective identifies the biological process, molecular function, and metabolic pathway of hub genes/oncogenes using the WebGestalt toolkit to understand the underlying interactive mechanisms extensively. Likewise, Liu et al. predicted the prognostic marker and hub genes by analyzing the hepatocellular and gynecological carcinoma TCGA samples. The study’s outcome stated that it might be helpful in personalized therapy to decide medical practice [27,28]. Feng et al. also described that identifying signature genes in ovarian cancer patients’ samples is vital as potential therapeutic targets [29]. Wang et al. determined that the vital role of hub genes in the progression of head and neck cancer might be used to invent the novel diagnostic method and treatment [30]. However, no extensive investigation has been conducted to determine the mutational essentiality and its impact on the overall survival rate of ESCC patients [31,32]. We integrated genome-wide DNA methylation and gene expression data to find DNA methylation indicators for early ESCC diagnosis.

The present study found the top ten hub genes and their essentiality in tumor growth by implementing publicly available a gene expression omnibus (GEO) dataset, GSE161533, which contains DEGs across 28 tumors and 28 normal tissues of the esophagus from ESCC patients. At the first top 250 DEG’s, they were predicted from the GEO dataset GSE161533 based on the p-value GEO2R analyzer. We constructed a protein–protein interaction (PPI) network for the top 250 DEGs using STRING 2.0 with hierarchical clustering and functional enrichment analyses of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. We identified HGsT10 with the help of CytoHubba. Likewise, we have analyzed the finding essentiality in expression at various stages, mutation and its effect on overall survival rate and promoter methylation with its histology.

Furthermore, the miRNAs that targeted the HGsT10 were also predicted to play central roles in cancer development via the regulation of its target genes. On the other hand, mutation-based studies often focus on a single gene or cancer type and do not contain extensive assessments of possible gene interactions that might predict survival. The overall survival analysis for HGsT10 and its variant was performed using the Gene Expression Profiling Interactive Analysis 2 (GEPIA2) and University of California, Santa Cruz (UCSC) Xena databases. Our study combined the multi-omics analysis platform, which included gene expression profiling microarrays, comparative analysis with various stages, functional enrichment, promoter methylation, wild and mutated type overall survival analysis, miRNA regulation. We suggested that our data provided detailed biological information about differentially expressed novel hub genes. These findings may be used to improve our understanding of ESCC to derive novel diagnostic methods and anti-cancer drugs against ESCC.

2. Materials and Methods

2.1. Affymetrix Microarray Data

The gene expression profile dataset no. GSE161533 deposited in the Gene Expression Omnibus (GEO) database (https://0-www-ncbi-nlm-nih-gov.brum.beds.ac.uk/geo/) (accessed on 17 December 2021) [33] by Qiu et al., based on the GPL570 platform (Affymetrix Gene Chip Human Genome U133 plus 2.0 Array), was subjected to bioinformatics analysis in the present study [34]. The dataset contained a total of 84 samples, including 28 esophageal squamous cell carcinoma tissues and 28 normal paired esophageal tissues, and 28 para-tumor esophageal tissues obtained from ESCC patients.

2.2. DEG’s Processing Using GEO2R

The GEO2R back-end transforms and analyses GEO data using well-known Bioconductor [35] R tools and then delivers the findings as a table of genes sorted by relevance displayed with GEO Profile visualizations. DEGs between ESCC tissues and normal esophageal tissues were identified via GEO2R online tools with |logFC| > 2 and adjusted p-value < 0.05 with default false discovery rate for multi testing. The top 250 genes were extracted from the whole dataset for further analysis. The DEGs of genes were designated as downregulated ones with log FC < 0, whereas upregulated ones were designated with log FC > 0.

2.3. PPI Network Construction to Identify Hub Genes

The gene network for the anticipated top 250 DEGs from the previous analysis was constructed using the STRING 2.0 database with a clustering coefficient of 0.395 and a PPI enrichment p-value of <1.0 × 10⁻¹⁶ [36]. The network was enriched functionally with Kyoto Encyclopedia of Genes and Genomes (KEGG), Orthology [37], and Gene Ontology [38]. The functionally enriched network was obtained and reconstructed with the CytoScape 3.8.2 with a clustering coefficient of 0.098. The hub genes were identified from the reconstructed network using the CytoHubba module for the highly interacted top 10 genes using the Maximal Clique Centrality (MCC) algorithm with the clustering coefficient of 0.898 [39,40].

2.4. Analysis and Validation of HGsT10 Expression and Its Various ESCC Stages

The predicted HGs10 expressions were analyzed in the ESCC patient’s sample data GSE161533 and then compared with TCGA (https://portal.gdc.cancer.gov) (accessed on 16 September 2021) samples of ESCC patients. The TCGA. samples were accessed via GEPIA2 (http://gepia.cancer-pku.cn/) (accessed on 16 September 2021) tool with the cutoff of |Log2FC| = 1.0 and p-value = 0.01 matched with TCGA normal samples and GTEx data [41]. The TCGA sample size for the normal is 286, and the ESCC tumor is 182. The HGsT10 expression was analyzed in the ESCC major stages may be implied to predict different stages of ESCC patients.

2.5. The Overall Survivability and Promoter Methylation

The overall survival (OS) rate of HGsT10 was analyzed. The TCGA patient samples were differentiated into two groups (high vs. low expression) and assessed by a Kaplan–Meier survival plot, with a log-rank p-value of 0.5 generated from GEPIA2. The survival rate was measured by,

$$\widehat{S} \left(t\right) ={{\displaystyle \prod }}_{i:ti \le t}\left(1-\frac{di}{ni}\right)$$

where, $\widehat{S} \left(t\right)$ = survival function; t_i = a time when at least one event happened; d_i = the number of deaths; and n_i = the individuals survived up to time [41]. The promoter methylation under the tumor histological condition was analyzed for the HGsT10 genes using the UALCAN server (http://ualcan.path.uab.edu) (Accessed: 10 December 2021). UALCAN is an online resource for exploring cancer transcriptome data [42]. TCGA data may investigate gene expression, promoter methylation, correlation, and prognosis across cancers.

2.6. Comparison of HGsT10 Gene Expression and Somatic Mutation Identification

From the University of California Santa Cruz (UCSC) Xena database (https://xenabrowser.net/datapages/) (Accessed: 10 December 2021) [43], the HTSeq counts of RNA-seq data for ESCC including 251 tumor samples were obtained (Fragments Per Kilobase of transcript per Million mapped reads upper quartile) with its log-transformed value as follow,

$$log2\left(fpkm-uq+1\right)$$

and analyzed. The overall survival data of the mutated HGsT10 gene for the ESCA patients were obtained. The Kaplan–Meier analysis with log-rank test was implemented to compare overall survival across different samples of ESCA.

2.7. HGsT10 Interaction with Other Genes and miRNA

The miRNAs are majorly involved in the pathophysiology of all types of human tumors. It was expressed between tumors and normal tissues with various differential rates in multiple cancers such as lung, glioblast, hepatocellular carcinomas, pancreatic, pituitary adenomas, cervical, brain, etc. However, the interaction between miRNAs and oncogenes remains unknown in ESCC. Here, the study constructed the miRNA- HGsT10 network to know the miRNA regulation in ESCC. The study used the miRNet 2.0 databases [44] to obtain the raw data and then reconstructed using CytoScape 3.8.2 with the help of the CytoHubba module ranked by the clustering coefficient algorithm and identified the top 10 hub miRNA (HmiRT10) for the predicted HGsT10 [38,39]. The network parameter was analyzed with the Network Analyzer plugin [45].

2.8. Functional and Pathway Enrichment Analysis for the Predicted HGsT10 Genes

DEGs from the two datasets were intersected to identify appropriate differential genes, which were then used to investigate the proper functioning of particular genes. Gene Ontology (GO) and KEGG enrichment analysis were used to explore molecular functions, biological pathways, and cellular components using the WebGestalt online tool (http://www.webgestalt.org) (Accessed: 17 December 2021) [46,47]. It uses the hypergeometric test to evaluate the significance of enrichment. The basic parameters were set to include: Homo sapiens as the organism of interest, over-representation analysis (ORA) as the method of claim, with the reference gene from lumina human ref 8.v3.and gene ID type as the gene symbol. The reference genes mapped with 24,770 Entrez gene IDs and 15,372 IDs were annotated to the selected functional categories of the predicted HGsT10. The significance of the enrichment ratio was calculated using the formula,

$$P={{\displaystyle \sum }}_{i=k}^n\frac{\left(\begin{array}{cc}m& -J\\ n& -i\end{array}\right)\left(\begin{array}{c}j\\ i\end{array}\right)}{\left(\begin{array}{c}m\\ n\end{array}\right)}$$

where, P = significance of enrichment; m = genes in the reference gene set (B); n—genes in the interesting gene set (A); Significant differences were set as FDR < 0.05.

3. Results

3.1. The GEO Dataset Validation

GSE161533 dataset was downloaded from the Gene Expression Omnibus (GEO) database and validated with GEO2R before starting the analysis. A volcano plot displays statistical significance (−log10 p-value) versus magnitude of change (log2 fold change), and the adjusted p-value set as <0.05 is useful for visualizing differentially expressed genes in normal and tumor samples (Figure 1A).

The sample groups are clean and separated significantly, which implements the Uniform Manifold Approximation and Projection method (UMAP) (Figure 1B). The downloaded dataset was analyzed to estimate the normalization generated using an R boxplot with its median-cantered values shown in Figure 1C. It describes that the selected data are normalized and cross-comparable. The results achieved with GEO2R in this study suggest that the selected samples can be used for further analysis.

3.2. Identification of DEGs and HGsT10 and PPI Network Construction for ESCC Tumor and Normal Tissues

Our present study has 28 ESCC tumor tissues and 28 normal esophageal tissues. Via GEO2R online tools, we extracted 250 DEGs from GSE161533. Results showed that a total of 113 DEGs were downregulated (logFC < 0), and 117 DEGs were upregulated genes (logFC > 0) in both tumor and normal tissues noted in Tables S1–S3 (Supplementary File S1). The extracted genes were constructed as a PPI network using the STRING 2.0 database, designed to predict gene–gene interactions. The STRING database was used to build a network model, and the genes with a score of ≥0.4 were chosen. The resulted network consisted of 195 nodes and 292 edges and an average local clustering coefficient of 0.395. The constructed network was reconstructed with the CytoScape 3.8.2. and HGsT10 were identified using the CytoHubba module (Figure 2A). It used the Maximal Clique Centrality (MCC) algorithm, the most effective method of finding hub nodes. The top ten genes ranked with MCC scores were deliberated as hub genes in this annotation. The hub genes selected from the PPI network using the MCC algorithm of the CytoHubba plugin are shown in Figure 2B.

According to the MCC sores, the top ten highest-scoring genes, namely matrix metallopeptidase 9 (MMP9), matrix metallopeptidase 3 (MMP3), matrix metallopeptidase 1 (MMP1), C-X-C motif chemokine ligand 8 (CXCL8), secreted phosphoprotein 1 (SPP1), urokinase plasminogen activator (PLAU), matrix metallopeptidase 13 (MMP13), collagen type I alpha 1 chain (COL1A1), collagen type I alpha two chains (COL1A2), and plasminogen activator, urokinase receptor (PLAUR), were identified as the hub genes.

3.3. Verification of the HGsT10 Expression Patterns

After the ten hub genes (MMP9, MMP3, MMP1, CXCL8, SPP1, PLAU, MMP13, COL1A1, COL1A2, and PLAUR) were screened out by the CytoHubba plugin, we verified the expression levels of the HGsT10 among the ESCC patient’s sample with the GEO dataset GSE161533. The HGsT10 expression was demonstrated among the ESCC patient’s sample of the TCGA database (Figure 3).

As shown in Figure 3A the MMP9 expression was noted maximum (198.822 TPM) in normal tissues (sample id GSM4909563). The samples GSM4909612, and GSM4909628 showed the expression of 2204.43, and 2231.53 TPM, respectively. Likewise, MMP3 (normal tissue sample max (NTSM), GSM4909557 = 9.8356; tumor tissues sample max (TTSM), GSM4909617 = 987.431), MMP1 (NTSM, GSM4909563 = 128.084; TTSM, GSM4909619 = 3143.52), CXCL8 (NTSM, GSM4909578 = 153.367; TTSM, GSM4909619 = 3034.75), SPP1 (NTSM, GSM4909566 = 106.629; TTSM, GSM4909612 = 2348.91), PLAU (NTSM, GSM4909569 = 133.366; TTSM, GSM4909632 = 1186.85), MMP13 (NTSM, GSM4909568 = 9.7126; TTSM, GSM4909629 = 313.765), COL1A1 (NTSM, GSM4909553 = 38.4753; TTSM, GSM4909622 = 1844.38), COL1A2 (NTSM, GSM4909566 = 2420.06; TTSM, GSM4909622 = 8155.97), and PLAUR (NTSM, GSM4909566 = 110.395; TTSM, GSM4909636 = 488.016 TPM) was upregulated in tumor tissues samples aberrantly, which might be a major cause in the ESCC tumor development. Furthermore, Figure 3B, describes that all of the HGsT10 was significantly upregulated in ESCC patient’s samples compared with normal subjects.

3.4. The Comparison of Expression of HGsT10 in the ESCC Patients with Other Types of EC

The identified gene expression data were represented as a heatmap, which is the common method of visualizing gene expression data (Figure 4A). The heatmap combined with clustering methods includes the identified HGsT10 genes and 160 samples together in the gene expression pattern. The expression pattern is displayed in a grid, and each row represents a gene, and each column represents a sample. The color red represents upregulated genes, blue represents downregulated genes, and white represents an unchanged expression.

Only two patients were affected by cystic, mucinous, and serous neoplasm (CMSN); hence the study excluded the expression comparison of HGsT10 in CMSN. The gene MMP9 was more upregulated in squamous cell neoplasms (SCN) than adenomas and adenocarcinomas (AAC) with the lower quartile (Q1) = −0.444, upper quartile (Q3) = 0.871 and median = 0.219; MMP3 expressed highly in AAC with Q1 = −0.571, Q3 = 0.752 and median = 0.168; MMP1 expressed majorly in SCN with Q1 = −0.705, Q3 = 2.25, and median = 0.279; CXCL8 expression slightly increased in AAC with Q1 = −0.551, Q3 = 0.702, and median 0.137; SPP1 was highly expressed in SCN with Q1 = −0.279, Q3 = 1.10, and median = 0.470; PLAU expressed highly in SCN with Q1 = −1.44, Q3 = 2.36, and median = 0.487; MMP13 expression is highly observed in SCN with Q1 = 0.0822, Q3 = 1.28, and median 0.473; COL1A1 highly expresses in SCN with Q1 = −0.272, Q3 = 1.13, and median = 0.314; Likewise, COL1A2 was expressed highly in SCN with Q1 = −0.561, Q3 = 1.13, and median 0.356; PLAUR expression was increased in AAC with Q1 = −0.477, Q3 = 1.08, and median = 0.532. It describes that genes MMP9, MMP1, SPP1, PLAU, MMP13, COL1A1, and COL1A2 expression might be increased in the advanced stage of ESCC tumor patients.

3.5. Correlation between HGST10 Genes Expression and Tumor Stage in Patients with ESCC Carcinoma

We have analyzed the tumor stage in ESCC patients and represented this in Figure 5. Using TCGA data, the identified HGsT10 gene expression was analyzed by main pathological stages I, II, III, and IV, to reveal the importance of genes for the concerns in the diagnosis and treatment of ESCC patients. The analyses demonstrated that the expression of MMP9, SPP1, PLAU, MMP13, COL1A1, and COL1A2 with the Pr-value 0.01, 0.000951, 0.0491, 0.0237, 0.0217, and 0.00258 (Pr-value < 0.05), respectively. The major stage plot annotation revealed their high expression levels significantly associated with the advanced ESCC stage. The remaining genes such as MMP3, MMP1, CXCL8, and PLAUR with the Pr-value of 0.421, 0.981, 0.369, and 0.596 are not significantly associated with the advanced stage of ESCC, respectively. In addition, OS analyses of the predicted HGsT10 genes were performed by the GEPIA2 database (Figure S1) to investigate the hub genes’ predictive values in the ESCC patients. The relationship of HGsT10 overexpression with prognosis in all 182 ESCC patients was also analyzed. ESCC patients were divided into low (91) and high (91) patients and identified the overall survival rate in months based on the gene expression median value. The results showed that patients with a high expression level of HGsT10 genes have a lower overall survival rate. It suggests that the inhibition of the expression level of the anticipated HGsT10 might reduce the tumor growth of ESCC patients.

3.6. The Validation of HGsT10 Expression by Promoter Methylation in ESCC Patients

The predicted HGsT10 promoter methylation under tumor histology was analyzed in ESCC patients using UALCAN server (Figure 6).

The results showed that the methylation was reduced in ESCC more than the normal (N) and adenocarcinoma (AC) for the following genes, MMP3 (N = 0.84; AC = 0.73; Squamous cell carcinoma (SCC) = 0.69) (Figure 6B), PLAU (N = 0.18; AC = 0.17; SCC = 0.17) (Figure 6F) and COL1A2 (N = 0.20; AC = 0.39; SCC = 0.17) (Figure 6I). The resulted median value of the methylated genes indicates higher expression in ESCC patients’ samples. Hence, the hypomethylated genes such as MMP3, PLAU, and COL1A2 could lead to cell proliferation, dysregulation of metabolic pathways, uncontrollable activation of tumor-relevant signaling networks, and cell cycle dysregulation. A convincing argument for a tumor-driving role of DNA methylation can be made when the same genes are also frequently mutated in cancer. Identification of somatic mutation and mRNA gene expression is crucial, which typically affects the survival outcomes of ESCA patients. We have also ensured the predicted survival outcomes are not generated randomly with null models while maintaining the correlations among variables. The predicted KM plot for the somatic mutation and its effect on the overall survival rate of the ESCA patients are shown in Figure S2 (Supplementary File S2). The mutation that occurs in the HGsT10 drastically decreases the patient’s survival rate. However, the mutation occurs only in smaller samples and reveals targetable combinations of mRNA expression events. Notably, COL1A2 and MMP9 mutated in six (1 Splice, 3 Missense/Inframe, 1 Silent, and 1 Intron/R.N.A.) and five (2 Missense/Inframe, 2Silent, and 1Intron/R.N.A.) genes samples resulted in the decreased survival rate with the p-value of 0.9460 and 0.06011. To our knowledge, this is the first systematic analysis of the association between genes wild-type and mutated type effects in the overall survivability of the ESCA patients.

3.7. The Predicted HmiRT10 from the Constructed miRNA-HGsT10 Regulatory Network

The miRNA-HGsT10 genes network was constructed to identify the greatest number of possible miRNAs involved in the pathogenesis of ESCC (Figure 7A).

The constructed network targeted 367 miRNA nodes and 597 edges. The top ten hub miRNAs (HmiRT10) were predicted with the constructed miRNA-HGsT10 regulatory network. The identified HmiRT10 are hsa-let-7b-5p, hsa-let-7f-5p, hsa-let-7c-5p, hsa-let-7d-5p, hsa-mir-19a-3p, hsa-mir-21-5p, hsa-let-7e-5p, hsa-let-7a-5p, hsa-mir-16-5p, and hsa-mir-19b-3p (Figure 7B). The parameter for the constructed HmiRT10 was tabulated (

Table 1. The network parameter for the constructed PPI network for the top 250 genes, predicted HGsT10, miRNA-HGsT10, predicted HmiRT10 were analyzed using the Network Analyzer plugin and tabulated.

S.No.	Name of the Parameter	Top 250 Genes (GSE161533 Samples)	Predicted HGs T10	miRNA Regulatory Network	Predicted HmiR T10
1.	Number of Nodes	129	10	377	10
2.	Number of edges	262	38	597	30
3.	Avg.number of neighbors	4.339	7.600	3.167	3.333
4.	Network diameter	13	2	4	4
5.	Network radius	7	1	3	3
6.	Characteristic path length	4.196	1.156	3.170	2.261
7.	Clustering coefficient	0.215	0.898	0.000	0.000
8.	Network density	0.037	0.844	0.008	0.196
9.	Network heterogeneity	0.180	0.169	3.880	0.640

).

3.8. Functional and Pathway Enrichment Analysis

The illumina_humanht_12_v3 was used as the reference to map 24,770 Entrez gene IDs, and 15,372 IDs are annotated to the selected functional categories for HGsT10 genes to enrich the pathway. The biological function was predicted with the KEGG ontology to know the functional integrity of the HGsT10 (Figure 8). Notably, most of the genes are involved in cell communication, metabolic process, stimulus-response, and biological regulation. Likewise, most of the genes located in the extracellular matrix is a key determinant of tumor development and prognosis and are involved in the protein-binding activity.

The enriched pathway significantly resulted in group one HGsT10 (CXCL8, MMP1 and, MMP9) and group two HGsT10 (CXCL8, MMP1, MMP13, MMP3 and, MMP9) are involved in the bladder cancer and IL-17 signaling pathway response with the enrichment ratio of 52.972 and 37.973. The predicted biological process (Figure 8A), cellular component (Figure 8B), and molecular function (Figure 8C). The enriched and annotated pathway was depicted in Figure 8D as a bar chart.

The resulting enriched pathway for the mapped HGsT10 genes is noted in

Table 2. The pathway enrichment analysis resulted in the following parameter for the predicted HGsT10 genes.

Group	Pathway Enriched	Gene Symbol	FDR	p-Value	Expected Value	Enrichment Ratio
1.	Bladder cancer	CXCL8, MMP1 and MMP9	0.0021221	0.000019649	0.056633	52.972
2.	IL-17 signaling pathway	CXCL8, MMP1, MMP13, MMP3, MMP9	0.000027535	8.4985 × 10−8	0.13167	37.973
3.	Relaxin signaling pathway	COL1A1, COL1A2, MMP1, MMP13, MMP9	0.000074171	4.5785 × 10−7	0.18406	27.165
4.	ECM-receptor interaction	COL1A1, COL1A2, SPP1	0.0098986	0.00015276	0.11185	26.822
5.	Rheumatoid arthritis	CXCL8, MMP1, MMP3	0.010116	0.00020353	0.12318	24.355
6,	Amoebiasis	COL1A1, COL1A2, CXCL8	0.010116	0.00024023	0.13026	23.032
7.	Prostate cancer	MMP3, MMP9, PLAU	0.010116	0.00026419	0.13450	22.304
8.	AGE-RAGE signaling pathway in diabetic complications	COL1A1, COL1A2, CXCL8	0.010116	0.00028100	0.13734	21.844
9.	Transcriptional misregulation in cancer	CXCL8, MMP3, MMP9, PLAU	0.0057635	0.000071154	0.25060	15.962
10.	Focal adhesion	COL1A1, COL1A2, SPP1	0.061616	0.0020919	0.27325	10.979

and Table S5 (Supplementary File S1).

4. Discussion

ESCA is the most common cause of cancer death and starts in the esophageal mucosa’s epithelial tissues. ESCA is the world’s eighth most prevalent cancer and the sixth most lethal worldwide. The most frequent histologic subtype of ESCA in ESCC. Most of the patients affected by EC are discovered at an advanced stage and are typically incurable. Despite extensive research into the mechanisms causing ESCC, reliable biomarkers for ESCC diagnosis, prognosis, and clinical targeting are still lacking. The mechanisms regulating ESCC have still not been completely understood [48,49]. In the present study, GEO dataset GSE161533 with high quality was selected to identify the hub genes associated with ESCC and their associated biological pathways by integrated multi-omics analysis. We chose the top 250 DEGs consisting of 136 upregulated genes, and 113 downregulated genes were identified, and they were used to construct the PPI network. Finally, we identified the HGsT10 significantly on the MCC scores in the PPI network.

The expression of HGsT10 genes (namely, MMP9, MMP3, MMP1, CXCL8, SPP1, PLAU, MMP13, COL1A1, COL1A2, and PLAUR) was identified in the GSE161533 sample and validated with the TGCA dataset. The expression rate of HGsT10 was increased for all the predicted hub genes. Notably, MMP1, SPP1, and MMP13 log-transformed TPM values are highly increased compared to other HGsT10 genes. The overexpression of MMP1 promoted tumor growth and metastasis by activating the PI3K/AKT pathway [50]. Furthermore, it makes a significant contribution to the extracellular matrix remodeling that occurs during tumor progression and metastasis [51]. Likewise, SPP1 plays a crucial role in the tumor growth of ESCC and can be used as a molecular biomarker in the early detection and diagnosis [52]. Similarly, the upregulation of MMP13 by activating the ERK1/2 pathway promoted tumor progression at metastasis and transformed squamous epithelial cells across various tissues [53,54,55].

The expression of HGsT10 genes based on the type of diseases was predicted with the GDC TCGA ESCA patients, including 160 patient samples. The expression level of MMP9, MMP1, SPP1, PLAU, MMP13, COL1A1 and COL1A2 genes was more highly increased in the squamous cell carcinoma than other ESCA types and worsened the overall survival rate of the patients. Most of the predicted HGsT10 genes are overexpressed early, such as MMP9, MMP1, CXCL8, SPP1, PLAU, MMP13, COL1A1, and PLAUR may be used as significant biomarkers to diagnose and treat the ESCC patients at an early stage. The overall survival rate of ESCA patients was drastically decreased with the increased expression of the predicted HGsT10 genes, especially CXCL8, SPP1, PLAU, MMP13, and PLAUR. Jiao et al. and Sedighi et al. stated that the MMP13 concentration in serum was used to diagnose and treat ESCC patients early [56,57]. Likewise, MMP-9, MMP1, and MMP13 serum concentrations were linked to clinical stages, and tumor size may decrease the overall survival rate of the ESCC patients [58,59,60]. The upregulation of CXCL8 in the microenvironment may lead to ESCC tumor development. Its regulation is triggered by the activation of the Akt and Erk1/2 pathways. It was accelerating the cancer cell movement and penetration. It also affects the overall survival rate of the patients drastically [61,62,63,64]. PLAU is abundantly elevated in ESCC, which triggers the formation of fibroblasts to inflammatory cancer-associated fibroblasts by activating the uPAR/Akt/NF-κB/IL8 pathway, which leads to a poor prognosis of ESCC patients [65].

The overexpression of MMP1 promoted tumor growth and metastasis by activating the PI3K/AKT pathway [50]. Furthermore, COL1A1 showed a poor survival rate associated with PI3K/Akt/mTOR, p53, apoptosis, and cell cycle pathway [66]. These findings imply that further research into the involvement of collagen family members in ESCC is needed in the future to advance ESCC identification and therapy [67]. Furthermore, the somatic mutations were identified and warranted meaningful association with decreased overall survival rate compared with wild-type. The mutated genes such as MMP9, PLAU, MMP13, and COL1A1 all drastically reduced the overall survival of ESCA patients. Booka et al. also identified that the somatic mutations in NFE2L2, CDKN2A and, KMT2D in ESCC significantly diminish overall survivability compared with its wild-type [68]. Our results suggest, the mutations, on the other hand, are found in only a few samples. Even with its overexpression rather than the mutation, the overall survival rate is reduced significantly.

The promoter methylation identified for the predicted HGsT10 genes based on the disease type by utilizing the TCGA ESCA samples consisted of 201 samples. The hypo-methylation was observed in MMP3, MMP1, CXCL8, SPP1, PLAU, MMP13, COL1A1, COL1A2, and PLAUR in ESCC patient’s samples when compared with normal samples. Our promoter methylation results suggested that hypo-methylated nine genes among HGsT10 are overexpressed in ESCC patient samples. Our findings suggest that the hypo-methylation at the promoter region directly correlated with overexpression, which leads to tumor proliferation and poorer prognosis of ESCC patients, as noted in the overall survival analysis of HGsT10. Li et al. and Cui et al. described that aberrant expression of overexpressed in the lung cancer 1 (OLC1) gene and the phospholipase C epsilon 1 (PLCE1) and nuclear factor-kappaB (NF-kappaB)-related proteins highly promote the ESCC tumorigenesis [69,70]. Likewise, Iwaya et al. described that the downregulation of the ST6GALNAC1 gene might be a possible causative gene for ESCC tumor proliferation due to its hypermethylation [71].

The MicroRNAs (miRNAs) regulate gene expressions in tumor and tumor-adjacent tissue [72]. Lorusso et al. explained that miRNAs could facilitate communication between tumor cells and their surroundings and control tumor immunity by directly regulating genes associated with immune stimulation or inactivation [73]. Here, we analyzed prognostic miRNAs of crucial histopathological subtypes ESCC on the predicted HgsT10 genes from the GSE161533 dataset. The constructed miRNA-HGsT10 network consisted of 367 miRNAs. The top ten hub miRNAs (HmiRT10) were identified, which regulates HGsT10 (namely hsa-let-7b-5p, hsa-let-7f-5p, hsa-let-7c-5p, hsa-let-7d-5p, hsa-mir-19a-3p, hsa-mir-21-5p, hsa-let-7e-5p, hsa-let-7a-5p, hsa-mir-16-5p and hsa-mir-19b-3p). The miRNA hsa-let-7b-5p, hsa-let-7f-5p, hsa-let-7c-5p, hsa-let-7d-5p and hsa-let-7a-5p regulates the target genes COL1A1, COL1A2 and CXCL8. The miRNA hsa-let-7e-5p regulates CXCL8, COL1A2, and MMP9. Let-7 family miRNAs are known to serve as tumor suppressor genes, and decreased expression of these miRNAs is linked to a poor prognosis in cancer patients [74,75,76,77].

The miRNA hsa-mir-21-5p regulates CXCL8 and MMP9 with the decreased overall survival rate of the ESCA patients in the TCGA dataset. Furthermore, the resulting mutations for CXCL8 and MMP9 clearly stated that the ESCC tumor progression was directly correlated with overexpression. Zheng et al., also stated that hsa-miR-21-5p upregulates the MMP9 gene in cardioembolic stroke patients [78]. Likewise, the miRNA hsa-mir-16-5p upregulated the HGsT10 such as CXCL8, COL1A1, MMP3, PLAUR, SPP1, and PLAU, which might affect the overall survival rate of the ESCC patients. The miRNA hsa-mir-19b-3p regulates PLAU and COL1A2, and hsa-mir-19a-3p regulates COL1A2 and PLAU, respectively. Jiang et al. identified that the miR-19b-3p significantly upregulated and led to colon cancer progression in 211 colon samples [79]. Similarly, overexpression of miR-19a was linked to a poor prognosis in various human cancers, including ESCC [80]. Our results revealed that the predicted HmiRT10 could be used as independent predictive biomarkers for the appropriate identification and treatment of ESCC.

GO analysis of the HGsT10 demonstrated that they were significantly enriched in the extracellular matrix (GO:0031012), biological processes of cell communication (GO:0007154), and protein binding (GO:0005515). The extracellular matrix has been linked to various elements of ESCA advancement, with this compartment being linked to several characteristics of EC genesis and progression [81]. The resulting pathway enrichment analysis evidenced that the HGsT10 genes are majorly involved in bladder cancer, IL-17 signaling pathway, and ECM-receptor interaction in rheumatoid arthritis (RA). Ahn et al. analyzed the association between breast cancer and RA using mendelian randomization (MR) analysis and reported that the underlying mechanism between these two diseases is still unclear [82].

5. Conclusions

Our study established a unique top ten signature genes and predicted the overall survival of ESCC, which can help in clinical decision-making for individual treatment by quantifying the expression of HGsT10 using quantitative real-time PCR (qPCR) technique. Accordingly, medications have to provide for the ESCC patients to save their lives. Additionally, independent miRNA datasets acquired from miRNet, the identified let-7 family and other miRNAs, mir-21-5p, miR-21-5p, mir-16-5p, mir-19b-3p, and mir-19a-3p were identified as independent prognostic miRNAs for ESCC overexpression of the predicted miRNAs, which might increase expression of HGsT10 within the GSE161533 and GDC TCGA samples. Our findings could also be clinically valuable in predicting the prognosis of ESCC in the advanced tumor stage. However, further experiments are warranted to confirm the miRNA utility in the expression of HGsT10, which is involved in the ESCC progression and development. Furthermore, further validation, such as small molecule discovery and transdifferentiation, are needed to develop appropriate medications to treat ESCC cancer patients even at the advanced stage.