# Identification of Molecular Biomarkers and Key Pathways for Esophageal Carcinoma (EsC): A Bioinformatics Approach

Authors: Md. Rakibul Islam; Mohammad Khursheed Alam; Bikash Kumar Paul; Deepika Koundal; Atef Zaguia; Kawsar Ahmed
Journal: BioMed Research International
 (2022)

---

## Abstract

Esophageal carcinoma (EsC) is a member of the cancer group that occurs in the esophagus; globally, it is known as one of the fatal malignancies. In this study, we used gene expression analysis to identify molecular biomarkers to propose therapeutic targets for the development of novel drugs. We consider EsC associated four different microarray datasets from the gene expression omnibus database. Statistical analysis is performed using R language and identified a total of 1083 differentially expressed genes (DEGs) in which 380 are overexpressed and 703 are underexpressed. The functional study is performed with the identified DEGs to screen significant Gene Ontology (GO) terms and associated pathways using the Database for Annotation, Visualization, and Integrated Discovery repository (DAVID). The analysis revealed that the overexpressed DEGs are principally connected with the protein export, axon guidance pathway, and the downexpressed DEGs are principally connected with the L13a-mediated translational silencing of ceruloplasmin expression, formation of a pool of free 40S subunits pathway. The STRING database used to collect protein-protein interaction (PPI) network information and visualize it with the Cytoscape software. We found 10 hub genes from the PPI network considering three methods in which the interleukin 6 (IL6) gene is the top in all methods. From the PPI, we found that identified clusters are associated with the complex I biogenesis, ubiquitination and proteasome degradation, signaling by interleukins, and Notch-HLH transcription pathway. The identified biomarkers and pathways may play an important role in the future for developing drugs for the EsC.


---

## Body

## 1. Introduction

Esophageal carcinoma (EsC) is a member of the cancer group that occurs in the esophagus; globally, it is known as one of the fatal malignancies. In the year of 2018, EsC ranked as the ninth most common type of cancer with 572,000 new cases (3.72% of all types of cancer cases) and the sixth most common form of cancer in mortality with 509,000 deaths [1]. EsC remains an endemic disease in several parts of the world especially in third world countries [2]. Though the incidence rates of EsC are unstable worldwide with the highest rates of incidence were found in Africa and eastern Asia [1]. Gender-wise studies claimed that around 70% of EsC patients are male [1]. Drinking alcohol and smoking are listed as risk factors for esophageal squamous cell carcinoma in the United States [3]. Gastroesophageal reflux disease (GERD) and Barrett’s esophagus are connected with an increased risk of the development of EsC [4, 5]. Obesity also accounts as a risk factor of esophagus-related adenocarcinoma [6]. EsC remains a global concern for its lower survival rate, 5-year survival rates until now stayed less than 20% [7]. Though a huge improvement had occurred in the medical field over the last few decades, the median survival rates of EsC have been slightly grown in the last few years [8]. Most of the EsC cases are diagnosed in its latter stages for the lack of early clinical symptoms. Some common symptoms are accounted such as sudden weight loss, breastbone burn feel, chest pain, and dysphagia. Microarray gene expression profile and gene chip analysis have been hugely applied in the medical field [9]. Gene expression analysis helps to decode differentially expressed genes and molecular biomarkers using several techniques that may have a potential influence on cancer development [10]. Molecular biomarkers acted a significant role with an early diagnostic and prognostic value in cancer treatment. A few studies have been produced to identify molecular biomarkers for EsC. In a study, Dong et al. showed that Methyltransferase Like 7B can take part in the early detection of esophageal adenocarcinoma [11]. Wang et al. claimed that the MAPK1 gene showed abnormal expression which may contribute to the development of EsC [12]. EsC is one of the cancers that take lots of attention from the researchers but still not much known about its mechanism and progression. The increasing study of EsC-associated molecular biomarkers may provide a foundation for unique approaches in preventing, diagnosing, and treating EsC. In this study, we have conducted a comprehensive microarray-based genome-wide analysis to identify molecular signatures using bioinformatics methods and tools. The current study is started by collecting 4 EsC-associated microarray datasets. We identified differentially expressed genes (DEGs) from datasets. DEGs are presented to complete functional study and protein-protein interaction analysis. Significant clusters are identified from protein interaction networks. We also identified hub genes using connectivity value, maximum neighborhood component (MNC), and bottleneck methods.


## 2. Methodology

### 2.1. Microarray Data Collection

Many studies have been conducted on esophageal cancer to explore genetic biomarkers [13–15]. But there are very few numbers of comprehensive analyses on EsC so that the exact genetic mechanisms are remained unknown till now. To explore genetic biomarkers, we applied a comprehensive analysis in our current study. We used four different microarray datasets to complete this study. GSE93756, GSE94012, GSE104958, and GSE143822 datasets are selected from National Center for Biotechnology Information’s (NCBI) Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) [16]. GSE93756 dataset has four samples based on platform GPL21282 Phalanx Human OneArray Ver. 7 Release 1. GSE94012 dataset has six samples based on platform GPL15207 [PrimeView] Affymetrix Human Gene Expression Array. GSE104958 dataset has a total of 46 samples, and the dataset is based on platform GPL21185 Agilent-072363 SurePrint G3 Human GE v3 8x60K Microarray 039494 (Probe Name Version) [17]. GSE143822 dataset has eight samples, and it is based on platform GPL20844 Agilent-072363 SurePrint G3 Human GE v3 8x60K Microarray 039494. Step by step process of this study is demonstrated in Figure 1.Figure 1
Flow diagram of this study. This diagram explains we start our first step of this study from the GEO database; from the database, we select 4 datasets for statistical analysis and identify DEGs maintaining our cut-off filtration. After that, we categorize the identified DEGs according to their expression (upregulated and downregulated). After categorization, we implement function analysis and protein-protein interaction analysis, which are the two most key analyses of this study.


### 2.2. Data Processing and DEG Identification

Limma stands for linear models for microarray data, and most of the functionality of limma has been developed for microarray data. Using limma for microarray data processing is simple, and its result is mostly accurate. We used the limma package of the R language to convert the raw files of our selected four datasets [18]. The datasets are converted into gene expression measures for further analysis. To identify statistical significance of genes log 2 FC foldchange>1.50 for overexpression, log2FC<−1.50 for downexpression, and standard adjusted P value < 0.05 are applied [19, 20].


### 2.3. GO and Pathway Enrichment Analysis of DEGs

Gene Ontology (GO) analysis provides wide biological exploration outcomes for a single gene or gene set. In recent years, GO analysis is a crucial part of system biology-related studies. In another corner, pathway enrichment analysis assists in explore mechanistically insight between gene sets produced from the wide genome-scale analysis [21]. In this study, we used the Gene Ontology database to explore DEGs associated GO terms [22], and pathway analysis is conducted using Kyoto Encyclopedia of Genes and Genomes (KEGG) [23], REACTOME [24], BIOCARTA [25], and Biological Biochemical Image Database (BBID) [26] databases. The Database for Annotation, Visualization, and Integrated Discovery (DAVID, http://david.abcc.ncifcrf.gov/) is fruitful to gather all outcomes [27]. Statistical significance P value < 0.05 is maintained for identifying the final outcomes.


### 2.4. PPI Construction and Clustering Analysis

The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING,https://string-db.org/) repository is used to explore internal interactions between DEGs [28]. A high combine score>0.70 is used to validate the interactions. Open-source software Cytoscape [29] is used to generate the protein-protein interaction (PPI) networks. CytoHubba plugin is applied to get topological parameter value [30]. To identify clusters from PPI networks, we used the Molecular Complex Detection (MCODE) algorithm [31]. The MCODE plugin built-in parameter is used for the analysis degree cutoff=2, node score cutoff=0.2, k−core=2, and maximumdepth=100 is counted as a minimum criterion. The functional pathway analysis in the cluster is performed by using the REACTOME database.


## 2.1. Microarray Data Collection

Many studies have been conducted on esophageal cancer to explore genetic biomarkers [13–15]. But there are very few numbers of comprehensive analyses on EsC so that the exact genetic mechanisms are remained unknown till now. To explore genetic biomarkers, we applied a comprehensive analysis in our current study. We used four different microarray datasets to complete this study. GSE93756, GSE94012, GSE104958, and GSE143822 datasets are selected from National Center for Biotechnology Information’s (NCBI) Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) [16]. GSE93756 dataset has four samples based on platform GPL21282 Phalanx Human OneArray Ver. 7 Release 1. GSE94012 dataset has six samples based on platform GPL15207 [PrimeView] Affymetrix Human Gene Expression Array. GSE104958 dataset has a total of 46 samples, and the dataset is based on platform GPL21185 Agilent-072363 SurePrint G3 Human GE v3 8x60K Microarray 039494 (Probe Name Version) [17]. GSE143822 dataset has eight samples, and it is based on platform GPL20844 Agilent-072363 SurePrint G3 Human GE v3 8x60K Microarray 039494. Step by step process of this study is demonstrated in Figure 1.Figure 1
Flow diagram of this study. This diagram explains we start our first step of this study from the GEO database; from the database, we select 4 datasets for statistical analysis and identify DEGs maintaining our cut-off filtration. After that, we categorize the identified DEGs according to their expression (upregulated and downregulated). After categorization, we implement function analysis and protein-protein interaction analysis, which are the two most key analyses of this study.


## 2.2. Data Processing and DEG Identification

Limma stands for linear models for microarray data, and most of the functionality of limma has been developed for microarray data. Using limma for microarray data processing is simple, and its result is mostly accurate. We used the limma package of the R language to convert the raw files of our selected four datasets [18]. The datasets are converted into gene expression measures for further analysis. To identify statistical significance of genes log 2 FC foldchange>1.50 for overexpression, log2FC<−1.50 for downexpression, and standard adjusted P value < 0.05 are applied [19, 20].


## 2.3. GO and Pathway Enrichment Analysis of DEGs

Gene Ontology (GO) analysis provides wide biological exploration outcomes for a single gene or gene set. In recent years, GO analysis is a crucial part of system biology-related studies. In another corner, pathway enrichment analysis assists in explore mechanistically insight between gene sets produced from the wide genome-scale analysis [21]. In this study, we used the Gene Ontology database to explore DEGs associated GO terms [22], and pathway analysis is conducted using Kyoto Encyclopedia of Genes and Genomes (KEGG) [23], REACTOME [24], BIOCARTA [25], and Biological Biochemical Image Database (BBID) [26] databases. The Database for Annotation, Visualization, and Integrated Discovery (DAVID, http://david.abcc.ncifcrf.gov/) is fruitful to gather all outcomes [27]. Statistical significance P value < 0.05 is maintained for identifying the final outcomes.


## 2.4. PPI Construction and Clustering Analysis

The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING,https://string-db.org/) repository is used to explore internal interactions between DEGs [28]. A high combine score>0.70 is used to validate the interactions. Open-source software Cytoscape [29] is used to generate the protein-protein interaction (PPI) networks. CytoHubba plugin is applied to get topological parameter value [30]. To identify clusters from PPI networks, we used the Molecular Complex Detection (MCODE) algorithm [31]. The MCODE plugin built-in parameter is used for the analysis degree cutoff=2, node score cutoff=0.2, k−core=2, and maximumdepth=100 is counted as a minimum criterion. The functional pathway analysis in the cluster is performed by using the REACTOME database.


## 3. Result and Demonstration

### 3.1. DEG Screening

Initially, a total of 20102, 20085, 33762, and 32212 DEGs are identified from GSE93756, GSE94012, GSE104958, and GSE143822 datasets. After applying the minimum log (FC) andP value criterion, 5802, 5393, 5945, and 7024 DEGs are identified correspondingly. 380 upregulated and 703 downregulated DEGs are screened out in selected four datasets that are used for further analysis (Table 1). The top 10 upregulated and downregulated DEGs are shown in Table 2.Table 1
Dataset analysis details (a) before filtration and (b) after filtration.
Accession numberAmount of sampleUpregulated DEGsDownregulated DEGsTotal DEGs(A) Before logFC filtrationGSE937564 samples61971390520102GSE940126 samples44641562120085GSE10495846 samples222381152433762GSE1438228 samples145521766032212Overlapped100338184821(B) After logFC filtrationGSE937564 samples109447085802GSE940126 samples152038735393GSE10495846 samples112846175945GSE1438228 samples84161837024Overlapped3807031083Table 2
Top 10 (a) upregulated and (b) downregulated DEG name and LogFC value.
DEG symbolLogFC(A) Upregulated DEGsAFG3L28.456189CAMKK28.268455EIF4H6.039133SLC6A195.983653OR2L35.71944FUNDC2P25.401373KRT6B5.195899WRAP535.106136OR56A35.054809LINC014655.031665(B) Downregulated DEGsCOPS5-11.2045C3orf59-9.20752NOX6-8.2488RAB3B-7.23927LINC01279-6.98733NEDD5-6.46107USP26-6.41274TTLL9-6.22158NKD1-6.09627FCAR-6.01488


### 3.2. GO and Pathway Enrichment Analysis of DEGs

We applied functional analysis using the DAVID database to achieve further knowledge into the function of identified DEGs. The functional analysis reveals significant enriched GO terms and pathways of identified DEGs. The GO analysis explores that the overexpressed DEGs are mainly associated with protein ubiquitination, and regulation of cell cycle for biological process (BP); endoplasmic reticulum membrane and nucleoplasm for cellular component (CC); and protein binding, DNA binding for molecular function (MF) (Table3, Figure 2). On another chapter of GO analysis explores the downexpressed DEGs associated with the translational initiation and SRP-dependent cotranslational protein targeting to membrane for BP; extracellular matrix, and ribosome for CC; structural constituent of ribosome, and NADH dehydrogenase (ubiquinone) activity for MF (Table 4, Figure 3).Table 3
Gene Ontology analysis of upregulated DEGs using DAVID functional tools.
CategoryGO IDGO termCount%P valueBPGO:0045892Negative regulation of transcription, DNA-templated220.0359890.001478BPGO:0046513Ceramide biosynthetic process50.0081790.001541BPGO:0016567Protein ubiquitination170.027810.00305BPGO:0051726Regulation of cell cycle90.0147230.003954BPGO:0048013Ephrin receptor signaling pathway70.0114510.008317BPGO:0016477Cell migration100.0163590.008953BPGO:0000045Autophagosome assembly50.0081790.009555BPGO:0045893Positive regulation of transcription, DNA-templated200.0327170.009702BPGO:0055007Cardiac muscle cell differentiation40.0065430.017191BPGO:0051865Protein autoubiquitination50.0081790.018845CCGO:0005789Endoplasmic reticulum membrane370.0605271.84E-05CCGO:0005654Nucleoplasm830.1357768.75E-05CCGO:0005634Nucleus1360.2224777.26E-04CCGO:0005802Trans-Golgi network100.0163590.001509CCGO:0031965Nuclear membrane130.0212660.002026CCGO:0005622Intracellular390.0637980.013664CCGO:0005737Cytoplasm1230.2012110.015049CCGO:0005829Cytosol800.1308690.036397CCGO:0005671Ada2/Gcn5/Ada3 transcription activator complex30.0049080.038709CCGO:0000139Golgi membrane190.0310810.046953MFGO:0005515Protein binding2370.3876981.25E-08MFGO:0003677DNA binding560.0916085.22E-04MFGO:0003684Damaged DNA binding70.0114510.002026MFGO:0004842Ubiquitin-protein transferase activity160.0261740.004151MFGO:0061630Ubiquitin protein ligase activity110.0179940.00625MFGO:0005070SH3/SH2 adaptor activity60.0098150.006276MFGO:0008565Protein transporter activity60.0098150.017551MFGO:0017137Rab GTPase binding80.0130870.022939MFGO:0003676Nucleic acid binding310.0507120.02606MFGO:0032794GTPase activating protein binding30.0049080.03379∗GO: Gene Ontology; ∗BP: biological process; ∗CC: cellular component; ∗MF: molecular function.Figure 2
Gene Ontology analysis of upregulated DEGs using DAVID functional tools. Different colors of dots mean different categories of GO terms. The green-colored dot indicates biological process, the blue-colored dot indicates cellular component, and the red-colored dot defines molecular functions. Thex-axis indicates the |Log (P value)| of associated GO terms. y-axis indicates the GO term name. The size of a dot represents gene count.Table 4
Gene Ontology analysis of downregulated DEGs using DAVID functional tools.
CategoryGO IDTermCount%P valueBPGO:0006413Translational initiation370.0143474.26E-09BPGO:0006614SRP-dependent cotranslational protein targeting to membrane270.0104691.96E-07BPGO:0006412Translation500.0193883.60E-07BPGO:0019083Viral transcription290.0112456.66E-07BPGO:0000184Nuclear-transcribed mRNA catabolic process, nonsense-mediated decay300.0116337.51E-07BPGO:0048245Eosinophil chemotaxis80.0031021.64E-06BPGO:0007155Cell adhesion690.0267554.88E-05BPGO:0006364rRNA processing380.0147351.20E-04BPGO:0002548Monocyte chemotaxis130.0050412.75E-04BPGO:0007156Homophilic cell adhesion via plasma membrane adhesion molecules290.0112455.11E-04CategoryTermCount%P valueCCGO:0031012Extracellular matrix550.0213263.51E-07CCGO:0005840Ribosome370.0143474.78E-07CCGO:0022625Cytosolic large ribosomal subunit200.0077555.24E-06CCGO:0030424Axon400.015513.48E-05CCGO:0005578Proteinaceous extracellular matrix450.0174496.17E-05CCGO:0005788Endoplasmic reticulum lumen350.0135719.14E-05CCGO:0005747Mitochondrial respiratory chain complex I140.0054292.80E-04CCGO:0022627Cytosolic small ribosomal subunit130.0050418.53E-04CCGO:0098793Presynapse150.0058160.001358CCGO:0015935Small ribosomal subunit90.003490.00193CategoryTermCount%P valueMFGO:0003735Structural constituent of ribosome490.0191.85E-08MFGO:0008137NADH dehydrogenase (ubiquinone) activity130.0050410.001124MFGO:0044822Poly(A) RNA binding1340.0519590.001979MFGO:0008237Metallopeptidase activity170.0065920.002751MFGO:0005201Extracellular matrix structural constituent150.0058160.002881MFGO:0003723RNA binding700.0271430.004922MFGO:00475553′,5′-cyclic-GMP phosphodiesterase activity60.0023270.006619MFGO:0042056Chemoattractant activity80.0031020.009717MFGO:0001077Transcriptional activator activity, RNA polymerase II core promoter proximal region sequence-specific binding340.0131840.01087MFGO:0008009Chemokine activity110.0042650.012935Figure 3
Gene Ontology analysis of downregulated DEGs using DAVID functional tools. Different colors of dots mean different categories of GO terms. The green-colored dot indicates biological process, the blue-colored dot indicates cellular component, and the red-colored dot defines molecular functions. Thex-axis indicates the |Log (P value)| of associated GO terms. y-axis indicates the GO term name. The size of a dot represents gene count.We used four different databases to achieve the associated pathways more clearly. The pathway analysis revealed that the overexpressed DEGs are principally connected with the protein export, axon guidance, and RHO GTPases Activate Formins pathway (Table5(a), Figure 4); the downexpressed DEGs are principally connected with the L13a-mediated translational silencing of ceruloplasmin expression, formation of a pool of free 40S subunits, and GTP hydrolysis and joining of the 60S ribosomal subunit pathways (Table 5(b), Figure 5).Table 5
Pathway enrichment analysis of (a) upregulated and (b) downregulated DEGs using DAVID functional tools.
Pathway termBenjaminiP valueSource(A) UpregulatedProtein export0.1699239.40E-04KEGGAxon guidance0.2848190.00338KEGGRHO GTPases Activate Formins0.9592430.007813REACTOMEHATs acetylate histones0.8826820.010449REACTOMESphingolipid metabolism0.5834730.013182KEGGPathogenic Escherichia coli infection0.5804960.017396KEGGGolgi associated vesicle biogenesis0.9758180.026998REACTOMEErbB signaling pathway0.6715620.027725KEGGSphingolipid signaling pathway0.6370040.030241KEGGXBP1(S) activates chaperone genes0.9569160.030359REACTOMELysosome0.5935130.031324KEGGThe information-processing pathway at the IFN-beta enhancer0.9783740.034876BIOCARTAActivation of RAC10.9611510.039023REACTOMESignaling of hepatocyte growth factor receptor0.8909620.040208BIOCARTAEpithelial cell signaling in helicobacter pylori infection0.6548120.042066KEGG(B) DownregulatedL13a-mediated translational silencing of ceruloplasmin expression3.71E-074.17E-10REACTOMEFormation of a pool of free 40S subunits4.51E-075.06E-10REACTOMEGTP hydrolysis and joining of the 60S ribosomal subunit4.84E-075.43E-10REACTOMERibosome3.48E-071.26E-09KEGGPeptide chain elongation4.55E-065.11E-09REACTOMESelenocysteine synthesis1.04E-051.17E-08REACTOMEEukaryotic translation termination1.37E-051.53E-08REACTOMEViral mRNA translation1.99E-052.24E-08REACTOMENonsense mediated decay (NMD) independent of the exon junction complex (EJC)2.30E-052.58E-08REACTOMESRP-dependent cotranslational protein targeting to membrane3.08E-043.45E-07REACTOMENonsense mediated decay (NMD) enhanced by the exon junction complex (EJC)3.78E-044.24E-07REACTOMEFormation of the ternary complex, and subsequently, the 43S complex0.0082659.31E-06REACTOMERibosomal scanning and start codon recognition0.0124291.40E-05REACTOMETranslation initiation complex formation0.0124291.40E-05REACTOMEChemokine0.0018863.85E-05BBIDFigure 4
Bar plot diagram to demonstrate pathway analysis outcomes of upregulated DEGs. Different color of bars indicates the different database name. Thex-axis indicates the value of |log10 (P value)|, and y-axis indicates the pathway term name.Figure 5
Bar plot diagram to demonstrate pathway analysis outcomes of downregulated DEGs. Different color of bars indicates the different database name. Thex-axis indicates the value of |log10 (P value)|, and y-axis indicates the pathway term name.


### 3.3. PPI Construction and Hub Gene Identifications

Using the STRING database, we generated the PPI network and visualized with Cytoscape software. Constructed PPI network has 646 nodes and 2055 connections, including 172 upregulated DEGs and 474 downregulated DEGs (Figure6). Using CytoHubba plugin, we identified the top 10 hub genes from the PPI network including IL6, CDH1, NOTCH1, ATP5C1, BPTF, MRPS11, MRPS15, MRPL1, NDUFB7, and NDUFS5. CytoHubba plugin has 11 different methods to identify significant genes from the PPI network; in this study, we consider three methods including connectivity value (degree), maximum neighborhood component (MNC), and bottleneck to identify hub genes. In the PPI network, the IL6 gene has the highest number of degree value 68, MNC value 60, and bottleneck value 151 (Figure 7). The top 10 hub gene name and their rank based on three methods are screened in Table 6.Figure 6
PPI network using identified DEGs. Nodes represent DEGs, and edge represents the connection between DEGs. The network has 646 nodes and 2055 connections. Green nodes represent upregulated DEGs, and red nodes represent downregulated DEGs. Eclipse-shaped nodes indicate the hub genes of the network. Hub genes are explored using 3 combined methods.Figure 7
Bar plot diagram to represent the values of degree, MNC, and bottleneck for specific hub genes. The red bar indicates degree value, the blue bar indicates MNC value, and the black bar indicates bottleneck value. Thex-axis represents the gene name, and the y-axis represents numerical values of the corresponding method.Table 6
Rank of 10 hub genes based on degree, MNC, and bottleneck methods.
Gene nameRank degreeRank MNCRank bottleneckIL6111CDH1232NOTCH1323ATP5C1445BPTF5104MRPS11656MRPS15768MRPL1879NDUFB7987NDUFS510910


### 3.4. Clustering Analysis

Cluster analysis is conducted using the MCODE method. In this analysis, 11 clusters are identified where the number of nodes is greater than 5. We identified four significant clusters from the constructed PPI network. The most significant cluster is enriched with MCODE score 17.5 and node density 33; 2nd significant cluster has MCODE score 12 and node density 12; 3rd significant cluster has MCODE score 9.238 and node density 22; the 4th significant cluster has MCODE score 5 and node density 9. Pathway enrichment analysis explored that clusters are significantly enriched with the complex I biogenesis, mitochondrial translation termination, ubiquitination and proteasome degradation, signaling by interleukins, and Notch-HLH transcription pathway (Table7). Cluster outcomes with their associated pathways are shown in Figure 8.Table 7
Associated pathways of significant 4 clusters.
Pathway termsCount%P value(A) Cluster 1Complex I biogenesis1545.454556.48E-24Mitochondrial translation termination1545.454556.35E-21Mitochondrial translation initiation1545.454556.35E-21Mitochondrial translation elongation1545.454556.35E-21Respiratory electron transport1339.393943.91E-17(B) Cluster 2Ubiquitination and proteasome degradation120.5931785.78E-17(C) Cluster 3Interferon alpha/beta signaling90.31591.32E-13Signaling by interleukins30.10530.003542ISG15 antiviral mechanism30.10530.008028(D) Cluster 4B-WICH complex positively regulates rRNA expression30.1948059.95E-05Notch-HLH transcription pathway20.129870.002863Figure 8
Top 4 clusters and their associated pathways for (a) cluster 1, (b) cluster 2, (c) cluster 3, and (d) cluster 4. Hexagonal-shaped nodes present pathway name, and eclipse-shaped presents the gene name.
(a)(b)(c)(d)


## 3.1. DEG Screening

Initially, a total of 20102, 20085, 33762, and 32212 DEGs are identified from GSE93756, GSE94012, GSE104958, and GSE143822 datasets. After applying the minimum log (FC) andP value criterion, 5802, 5393, 5945, and 7024 DEGs are identified correspondingly. 380 upregulated and 703 downregulated DEGs are screened out in selected four datasets that are used for further analysis (Table 1). The top 10 upregulated and downregulated DEGs are shown in Table 2.Table 1
Dataset analysis details (a) before filtration and (b) after filtration.
Accession numberAmount of sampleUpregulated DEGsDownregulated DEGsTotal DEGs(A) Before logFC filtrationGSE937564 samples61971390520102GSE940126 samples44641562120085GSE10495846 samples222381152433762GSE1438228 samples145521766032212Overlapped100338184821(B) After logFC filtrationGSE937564 samples109447085802GSE940126 samples152038735393GSE10495846 samples112846175945GSE1438228 samples84161837024Overlapped3807031083Table 2
Top 10 (a) upregulated and (b) downregulated DEG name and LogFC value.
DEG symbolLogFC(A) Upregulated DEGsAFG3L28.456189CAMKK28.268455EIF4H6.039133SLC6A195.983653OR2L35.71944FUNDC2P25.401373KRT6B5.195899WRAP535.106136OR56A35.054809LINC014655.031665(B) Downregulated DEGsCOPS5-11.2045C3orf59-9.20752NOX6-8.2488RAB3B-7.23927LINC01279-6.98733NEDD5-6.46107USP26-6.41274TTLL9-6.22158NKD1-6.09627FCAR-6.01488


## 3.2. GO and Pathway Enrichment Analysis of DEGs

We applied functional analysis using the DAVID database to achieve further knowledge into the function of identified DEGs. The functional analysis reveals significant enriched GO terms and pathways of identified DEGs. The GO analysis explores that the overexpressed DEGs are mainly associated with protein ubiquitination, and regulation of cell cycle for biological process (BP); endoplasmic reticulum membrane and nucleoplasm for cellular component (CC); and protein binding, DNA binding for molecular function (MF) (Table3, Figure 2). On another chapter of GO analysis explores the downexpressed DEGs associated with the translational initiation and SRP-dependent cotranslational protein targeting to membrane for BP; extracellular matrix, and ribosome for CC; structural constituent of ribosome, and NADH dehydrogenase (ubiquinone) activity for MF (Table 4, Figure 3).Table 3
Gene Ontology analysis of upregulated DEGs using DAVID functional tools.
CategoryGO IDGO termCount%P valueBPGO:0045892Negative regulation of transcription, DNA-templated220.0359890.001478BPGO:0046513Ceramide biosynthetic process50.0081790.001541BPGO:0016567Protein ubiquitination170.027810.00305BPGO:0051726Regulation of cell cycle90.0147230.003954BPGO:0048013Ephrin receptor signaling pathway70.0114510.008317BPGO:0016477Cell migration100.0163590.008953BPGO:0000045Autophagosome assembly50.0081790.009555BPGO:0045893Positive regulation of transcription, DNA-templated200.0327170.009702BPGO:0055007Cardiac muscle cell differentiation40.0065430.017191BPGO:0051865Protein autoubiquitination50.0081790.018845CCGO:0005789Endoplasmic reticulum membrane370.0605271.84E-05CCGO:0005654Nucleoplasm830.1357768.75E-05CCGO:0005634Nucleus1360.2224777.26E-04CCGO:0005802Trans-Golgi network100.0163590.001509CCGO:0031965Nuclear membrane130.0212660.002026CCGO:0005622Intracellular390.0637980.013664CCGO:0005737Cytoplasm1230.2012110.015049CCGO:0005829Cytosol800.1308690.036397CCGO:0005671Ada2/Gcn5/Ada3 transcription activator complex30.0049080.038709CCGO:0000139Golgi membrane190.0310810.046953MFGO:0005515Protein binding2370.3876981.25E-08MFGO:0003677DNA binding560.0916085.22E-04MFGO:0003684Damaged DNA binding70.0114510.002026MFGO:0004842Ubiquitin-protein transferase activity160.0261740.004151MFGO:0061630Ubiquitin protein ligase activity110.0179940.00625MFGO:0005070SH3/SH2 adaptor activity60.0098150.006276MFGO:0008565Protein transporter activity60.0098150.017551MFGO:0017137Rab GTPase binding80.0130870.022939MFGO:0003676Nucleic acid binding310.0507120.02606MFGO:0032794GTPase activating protein binding30.0049080.03379∗GO: Gene Ontology; ∗BP: biological process; ∗CC: cellular component; ∗MF: molecular function.Figure 2
Gene Ontology analysis of upregulated DEGs using DAVID functional tools. Different colors of dots mean different categories of GO terms. The green-colored dot indicates biological process, the blue-colored dot indicates cellular component, and the red-colored dot defines molecular functions. Thex-axis indicates the |Log (P value)| of associated GO terms. y-axis indicates the GO term name. The size of a dot represents gene count.Table 4
Gene Ontology analysis of downregulated DEGs using DAVID functional tools.
CategoryGO IDTermCount%P valueBPGO:0006413Translational initiation370.0143474.26E-09BPGO:0006614SRP-dependent cotranslational protein targeting to membrane270.0104691.96E-07BPGO:0006412Translation500.0193883.60E-07BPGO:0019083Viral transcription290.0112456.66E-07BPGO:0000184Nuclear-transcribed mRNA catabolic process, nonsense-mediated decay300.0116337.51E-07BPGO:0048245Eosinophil chemotaxis80.0031021.64E-06BPGO:0007155Cell adhesion690.0267554.88E-05BPGO:0006364rRNA processing380.0147351.20E-04BPGO:0002548Monocyte chemotaxis130.0050412.75E-04BPGO:0007156Homophilic cell adhesion via plasma membrane adhesion molecules290.0112455.11E-04CategoryTermCount%P valueCCGO:0031012Extracellular matrix550.0213263.51E-07CCGO:0005840Ribosome370.0143474.78E-07CCGO:0022625Cytosolic large ribosomal subunit200.0077555.24E-06CCGO:0030424Axon400.015513.48E-05CCGO:0005578Proteinaceous extracellular matrix450.0174496.17E-05CCGO:0005788Endoplasmic reticulum lumen350.0135719.14E-05CCGO:0005747Mitochondrial respiratory chain complex I140.0054292.80E-04CCGO:0022627Cytosolic small ribosomal subunit130.0050418.53E-04CCGO:0098793Presynapse150.0058160.001358CCGO:0015935Small ribosomal subunit90.003490.00193CategoryTermCount%P valueMFGO:0003735Structural constituent of ribosome490.0191.85E-08MFGO:0008137NADH dehydrogenase (ubiquinone) activity130.0050410.001124MFGO:0044822Poly(A) RNA binding1340.0519590.001979MFGO:0008237Metallopeptidase activity170.0065920.002751MFGO:0005201Extracellular matrix structural constituent150.0058160.002881MFGO:0003723RNA binding700.0271430.004922MFGO:00475553′,5′-cyclic-GMP phosphodiesterase activity60.0023270.006619MFGO:0042056Chemoattractant activity80.0031020.009717MFGO:0001077Transcriptional activator activity, RNA polymerase II core promoter proximal region sequence-specific binding340.0131840.01087MFGO:0008009Chemokine activity110.0042650.012935Figure 3
Gene Ontology analysis of downregulated DEGs using DAVID functional tools. Different colors of dots mean different categories of GO terms. The green-colored dot indicates biological process, the blue-colored dot indicates cellular component, and the red-colored dot defines molecular functions. Thex-axis indicates the |Log (P value)| of associated GO terms. y-axis indicates the GO term name. The size of a dot represents gene count.We used four different databases to achieve the associated pathways more clearly. The pathway analysis revealed that the overexpressed DEGs are principally connected with the protein export, axon guidance, and RHO GTPases Activate Formins pathway (Table5(a), Figure 4); the downexpressed DEGs are principally connected with the L13a-mediated translational silencing of ceruloplasmin expression, formation of a pool of free 40S subunits, and GTP hydrolysis and joining of the 60S ribosomal subunit pathways (Table 5(b), Figure 5).Table 5
Pathway enrichment analysis of (a) upregulated and (b) downregulated DEGs using DAVID functional tools.
Pathway termBenjaminiP valueSource(A) UpregulatedProtein export0.1699239.40E-04KEGGAxon guidance0.2848190.00338KEGGRHO GTPases Activate Formins0.9592430.007813REACTOMEHATs acetylate histones0.8826820.010449REACTOMESphingolipid metabolism0.5834730.013182KEGGPathogenic Escherichia coli infection0.5804960.017396KEGGGolgi associated vesicle biogenesis0.9758180.026998REACTOMEErbB signaling pathway0.6715620.027725KEGGSphingolipid signaling pathway0.6370040.030241KEGGXBP1(S) activates chaperone genes0.9569160.030359REACTOMELysosome0.5935130.031324KEGGThe information-processing pathway at the IFN-beta enhancer0.9783740.034876BIOCARTAActivation of RAC10.9611510.039023REACTOMESignaling of hepatocyte growth factor receptor0.8909620.040208BIOCARTAEpithelial cell signaling in helicobacter pylori infection0.6548120.042066KEGG(B) DownregulatedL13a-mediated translational silencing of ceruloplasmin expression3.71E-074.17E-10REACTOMEFormation of a pool of free 40S subunits4.51E-075.06E-10REACTOMEGTP hydrolysis and joining of the 60S ribosomal subunit4.84E-075.43E-10REACTOMERibosome3.48E-071.26E-09KEGGPeptide chain elongation4.55E-065.11E-09REACTOMESelenocysteine synthesis1.04E-051.17E-08REACTOMEEukaryotic translation termination1.37E-051.53E-08REACTOMEViral mRNA translation1.99E-052.24E-08REACTOMENonsense mediated decay (NMD) independent of the exon junction complex (EJC)2.30E-052.58E-08REACTOMESRP-dependent cotranslational protein targeting to membrane3.08E-043.45E-07REACTOMENonsense mediated decay (NMD) enhanced by the exon junction complex (EJC)3.78E-044.24E-07REACTOMEFormation of the ternary complex, and subsequently, the 43S complex0.0082659.31E-06REACTOMERibosomal scanning and start codon recognition0.0124291.40E-05REACTOMETranslation initiation complex formation0.0124291.40E-05REACTOMEChemokine0.0018863.85E-05BBIDFigure 4
Bar plot diagram to demonstrate pathway analysis outcomes of upregulated DEGs. Different color of bars indicates the different database name. Thex-axis indicates the value of |log10 (P value)|, and y-axis indicates the pathway term name.Figure 5
Bar plot diagram to demonstrate pathway analysis outcomes of downregulated DEGs. Different color of bars indicates the different database name. Thex-axis indicates the value of |log10 (P value)|, and y-axis indicates the pathway term name.


## 3.3. PPI Construction and Hub Gene Identifications

Using the STRING database, we generated the PPI network and visualized with Cytoscape software. Constructed PPI network has 646 nodes and 2055 connections, including 172 upregulated DEGs and 474 downregulated DEGs (Figure6). Using CytoHubba plugin, we identified the top 10 hub genes from the PPI network including IL6, CDH1, NOTCH1, ATP5C1, BPTF, MRPS11, MRPS15, MRPL1, NDUFB7, and NDUFS5. CytoHubba plugin has 11 different methods to identify significant genes from the PPI network; in this study, we consider three methods including connectivity value (degree), maximum neighborhood component (MNC), and bottleneck to identify hub genes. In the PPI network, the IL6 gene has the highest number of degree value 68, MNC value 60, and bottleneck value 151 (Figure 7). The top 10 hub gene name and their rank based on three methods are screened in Table 6.Figure 6
PPI network using identified DEGs. Nodes represent DEGs, and edge represents the connection between DEGs. The network has 646 nodes and 2055 connections. Green nodes represent upregulated DEGs, and red nodes represent downregulated DEGs. Eclipse-shaped nodes indicate the hub genes of the network. Hub genes are explored using 3 combined methods.Figure 7
Bar plot diagram to represent the values of degree, MNC, and bottleneck for specific hub genes. The red bar indicates degree value, the blue bar indicates MNC value, and the black bar indicates bottleneck value. Thex-axis represents the gene name, and the y-axis represents numerical values of the corresponding method.Table 6
Rank of 10 hub genes based on degree, MNC, and bottleneck methods.
Gene nameRank degreeRank MNCRank bottleneckIL6111CDH1232NOTCH1323ATP5C1445BPTF5104MRPS11656MRPS15768MRPL1879NDUFB7987NDUFS510910


## 3.4. Clustering Analysis

Cluster analysis is conducted using the MCODE method. In this analysis, 11 clusters are identified where the number of nodes is greater than 5. We identified four significant clusters from the constructed PPI network. The most significant cluster is enriched with MCODE score 17.5 and node density 33; 2nd significant cluster has MCODE score 12 and node density 12; 3rd significant cluster has MCODE score 9.238 and node density 22; the 4th significant cluster has MCODE score 5 and node density 9. Pathway enrichment analysis explored that clusters are significantly enriched with the complex I biogenesis, mitochondrial translation termination, ubiquitination and proteasome degradation, signaling by interleukins, and Notch-HLH transcription pathway (Table7). Cluster outcomes with their associated pathways are shown in Figure 8.Table 7
Associated pathways of significant 4 clusters.
Pathway termsCount%P value(A) Cluster 1Complex I biogenesis1545.454556.48E-24Mitochondrial translation termination1545.454556.35E-21Mitochondrial translation initiation1545.454556.35E-21Mitochondrial translation elongation1545.454556.35E-21Respiratory electron transport1339.393943.91E-17(B) Cluster 2Ubiquitination and proteasome degradation120.5931785.78E-17(C) Cluster 3Interferon alpha/beta signaling90.31591.32E-13Signaling by interleukins30.10530.003542ISG15 antiviral mechanism30.10530.008028(D) Cluster 4B-WICH complex positively regulates rRNA expression30.1948059.95E-05Notch-HLH transcription pathway20.129870.002863Figure 8
Top 4 clusters and their associated pathways for (a) cluster 1, (b) cluster 2, (c) cluster 3, and (d) cluster 4. Hexagonal-shaped nodes present pathway name, and eclipse-shaped presents the gene name.
(a)(b)(c)(d)


## 4. Discussion

Globally EsC is considered one of the most deadly diseases for its fast development and base presage. Around 80% of EsC cases are recorded from less developed regions in the world [2]. In 2012 in China, EsC had listed the fifth common diagnosed cancer type and the fourth eminent cause of mortality [32]. It is urgent to understand the clinical epidemiology of EsC to develop medical treatment. In this study, we developed a microarray gene profile analysis to identify molecular signatures. EsC-associated four different datasets GSE93756, GSE94012, GSE104958, and GSE143822 are selected, and these datasets are analyzed with the limma package of R language. 380 upregulated and 703 downregulated DEGs are matched in all datasets following every criterion. These DEGs are applied to draw significant GO terms using the DAVID database. GO analysis shows that the upregulated DEGs are associated with protein ubiquitination, regulation of cell cycle, endoplasmic reticulum membrane, nucleoplasm, and protein binding. The downregulated DEGs are associated with translational initiation, SRP-dependent cotranslational protein targeting to membrane, extracellular matrix, ribosome, and structural constituent of ribosome. Cell cycle abnormalities had been indicated as a key factor of esophagus tumorigenesis [33, 34]. In 2017, Otto et al. claimed that the cell cycle protein may play a promising role in cancer therapy [35].In this study, PPI network is constructed by using identified DEGs. From the PPI network, we found 10 hub genes (IL6, CDH1, NOTCH1, ATP5C1, BPTF, MRPS11, MRPS15, MRPL1, NDUFB7, and NDUFS5) using three combined methods. Interleukin 6 (IL6) gene is a member of the Interleukin family, and it takes part in cell growth operation. IL6 can act as both a proinflammatory cytokine and an anti-inflammatory myokine, and it is associated with many types of cancer development [36]. A study showed that breast cancer cells produced IL6 as a core compound [37]. IL6 also listed as a therapeutic biomarker in renal cell carcinoma [38]. IL6 shows poor prognosis values in lung cancer patients [39]. IL6-associated signaling pathways also take part in cancer progression. Based on the above discussion, we can say that IL6 may play a significant role in EsC progression. Cadherin 1 (CDH1) gene is connected with protein-coding. CDH1 is associated with the cell proliferation pathway, which plays an important preface in cancer development [40]. Mutations of CDH1 protein marked as an increased risk factor for hereditary diffuse gastric cancer (HDRC) [41, 42].HDRC affected women to embrace a high risk of having breast cancer [43]. HDRC patients increased high risk of developing stomach cancer which is associated with the esophagus organ. Several characteristics indicate that CDH1 may take part in the development of EsC. NOTCH1 is known for encoding the NOTCH family of proteins. NOTCH1 plays a role in cell growth and proliferation, differentiation, and apoptosis. NOTCH1 is engaged in many types of cancer, including triple-negative breast cancer, leukemia, brain tumors, and many others. It influences apoptosis, proliferation, immune response, and the population of cancer stem cells [44]. Regarding the above discussion, we can assume NOTCH1 may impact EsC development. The Bromodomain PHD Finger Transcription Factor (BPTF) gene was found overexpressed and showed poor prognosis value in the tissue of lung adenocarcinoma [45]. A study from 2015 proposed BPTF as a novel target for anticancer therapy [46].In the PPI analysis section, we applied the MCODE method to identify clusters. Significant four clusters are identified, and pathway analysis is performed. Pathway analysis showed that the clusters are principally enriched with complex I biogenesis, mitochondrial translation termination, mitochondrial translation initiation, and interferon-alpha/beta signaling pathway. Mitochondrial biogenesis develops breast cancer tumors in the epithelial cell lines [47].The authors believe the outcomes of this study will make an impact on the biomarker identification of EsC. But more studies are required to prove the statement. Lack of tools and established laboratory, we could not verify our outcomes which is the limitation of this study. For future goals, we will use the outputs to explore microRNA biomarkers for EsC, which will give us deeper knowledge regarding EsC development.



---

*Source: 5908402-2022-01-12.xml*