Esophageal squamous cell carcinoma: Integrated bioinformatics analysis for differential gene expression with identification of hub genes and lncRNA.

Background: Esophageal squamous cell carcinoma (ESCC) is a typical Gastro-Intestinal (GI) tract neoplasm. This study was conducted to know the Differential Expressed Genes (DEGs) profile of ESCC along with hub gene screening, lncRNA identification, and drug-genes interactions.
Methods: GSE161533, GSE20347, GSE45670 microarray datasets were retrieved from the NCBI Gene Expression Omnibus (GEO) database. GEO2R was used for the DEGs identification, whereas GO (Gene Ontology) and KEGG enrichment analysis were performed in DAVID. PPI network constructed using STRING and visualized with Cytoscape app with the help of MCODE. The top ten connectivity genes were selected as hub genes-further survival analysis was performed in the Kaplan-Meier plotter. Moreover, Boxplot, pathological stage plots were constructed using GEPIA (Gene Expression Profiling Interactive Analysis). The methylation heatmap assembled in the DiseaseMeth version 2.0. lncRNA (Long non-coding RNA) was identified comparing the list of genes in HUGO, and Gene-drug interactions were accumulated from the DgiDB platform.
Results: This experiment showed 16 upregulated, and 59 downregulated DEGs shared among the three datasets. Biological process analysis showed significant terms such as extracellular matrix disassembly and collagen catabolism. The extracellular region was detected as the most crucial cellular compartment. Notably, metalloen dopeptidease and serine-type endopeptidase activity showed significant molecular functions term. In contrast, transcriptional misregulation was a highly substantial KEGG pathway. Kaplan-Meier plotter showed higher expression of CXCL8, SPP1, MMP13, CXCL1, and TOP2A have a significant impact on the overall survival of the patients. Nine out of ten hub genes have significantly different expression levels than normal and cancer tissues. HYMAI was the only lncRNA commonly expressed upregulated among the three datasets. Drug-gene interaction showed multiple genes have no drug options exist till now.

Introduction

Cancer is the cause of a significant number of mortality worldwide. Cancer can be defined as uncontrolled cell growth, and almost every tissue can be affected by this disease [2]. Esophageal Squamous Cell Carcinoma (ESCC) is considered one of the most common gastrointestinal (GI) neoplasms worldwide [1]. In most cases of ESCC, symptoms did not show early, resulting in a higher death rate due to limited treatment regimens in the late stages [11]. Several risk factors were identified for ESCC, including alcohol consumption, tobacco products consumption, smoking, lower fiber intake, etc. Upper Body Mass Index (BMI) and micronutrient deficiency were also hazardous [4].

Due to the present day's advancement of genomic techniques such as microarray analysis and high throughput analysis, genes associated with ESCC are now a topic of interest to discover the specific genes with their expression correlation to the tumor [14]. Differential gene expression and their related activity are the fundamental way to understand the mechanism of disease advancement. Many genes and their co-expression were regularly identified for the progression of ESCC worldwide. Some of the genes showed significant expression results in corresponding survival analyses. CDK1 and TOP2A were analyzed as the critical genes for ESCC neoplasm by Yang and his group (W. [35]. Whereas CDCA5 was considered the crucial gene for the prognosis of ESCC by another author [33]. CFLAR, LAMA5, ITGA6, ITGB4, and SDC4 genes were also validated for ESCC progression (L. [36]. Whereas, there were also few Long Non-coding RNA (lncRNA) identified, which might impact ESCC pathogenesis, and the BANCR gene was identified by another author [26]. HCG22 was also detected as a lncRNA by a previous author (X. [15].

Many publicly available microarray and high throughput genomic data are available, but there is a lack of bioinformatics analysis and correlation of the disease occurrence. That analysis can quickly identify the potential genes associated with cancer or tumor. The following study used three microarray datasets to place the common differential gene expression with functional identification of Esophageal Squamous Cell Carcinoma (ESCC) hub genes.

Materials and Method

Data source: Three microarray datasets (GSE161533, GSE20347, and GSE45670) were collected from Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/). 28 normal tissue and 28 tumor tissue of esophageal squamous cell carcinoma data were from the GSE161533 dataset. Besides, the GSE20347 dataset contributed 17 samples for both the normal and tumor tissue. At the same time, the GSE45670 dataset has ten normal samples with 28 ESCC samples.

Differential gene expression: GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r/) was used for the differential gene expression analysis (J. [16]. All the statistical values were auto-selected by the GEO. For the selection of differentially expressed, the genes should qualify the cutoff criteria of Log Fold Change (|logFC|) ≥ 1.5 and the adj. p-value ≤ 0.05. All datasets were analyzed with identical selection criteria. HemI 2.0 - Heatmap Illustrator 2.0 software (http://hemi.biocuckoo.cn:81/) was used to compare the three datasets' expression values [7].

Annotation of DEGs: DAVID platform (https://david.ncifcrf.gov/) was used to annotate the 75 common, regulated DEGs [18]. Upregulated and downregulated genes GO analysis were taken into account, with the p-value less than 0.01 considered significant. Biological process (BP), Cellular Compartment (CC), and Molecular Function (MF) were analyzed. At the same time, KEGG pathway analysis was selected with similar p-values.

Protein-protein interaction (PPI) network of DEGs: STRING (https://string-db.org/) is the online bioinformatics tool to ascertain the hub gene and examine the interactions between the genes [37]. The interaction score and the maximum number of interactions were >0.4 and 10, respectively.

Selection of the hub genes: The highest connectivity of the correlated genes was calculated from the Cytoscape software (https://cytoscape.org/), an available bioinformatic analytical tool. MCODE app used for the selection of hub gene with following criteria- MCODE score >5, degree cut off = 2, node score cut off = 0.2, Max depth = 100, k-score = 2 [41]. The top ten connected genes were considered as hub genes for this study.

Survival analysis: Kaplan-Meier survival analysis was performed in the online analytical tool (https://kmplot.com/analysis/) [28]. Overall survival analysis was performed according to the collected gastric cancer database of the website.

Box plot and pathological stage plot expression comparison: The GEPIA (http://gepia.cancer-pku.cn/index.html) tool was used to construct and compare the boxplots. Overall survival method was taken into consideration along with median group cutoff (50%), Hazards Ratio (HR), and 95% confidence interval [17]. |LogFC| Cutoff value was 1.5. Hub gene pathological stage plot was also constructed to compare the stages of the diseases.

Disease methylation: Normal and diseased methylation data compared in the DiseaseMeth Version 2 (http://bio-bigdata.hrbmu.edu.cn/diseasemeth/index.html) [32]. All default criteria were selected for the analysis.

lncRNA identification: Differentially expressed lncRNA identified compared to the approved lncRNA list from the HUGO database (https://www.genenames.org/) [25]. |LogFC| > 0.5 and adj p-value less than 0.05 were considered as significant during the identification.

Drugs and genes interactions: The selected 75 genes were analyzed according to their interactions with the currently available and approved drugs. DGIdb (https://www.dgidb.org/) is an open, accessible public repository for the identification of drugs and genes interactions [9]. The exchanges were visualized through Cytoscape software.

Results

Identification of the DEGs: The number of upregulated and downregulated genes among the three datasets were shown in the Venn diagram (Fig. 1A and B). There were 506, 549, and 1768 significant DEGs for the GSE161533, GSE20347, and GSE45670 datasets with 16 upregulated and 59 downregulated genes in common. Fig. 2 depicts the comparative values of the three datasets. All the gene names are given in Table 1.

Fig. 1

DEGs among the three datasets (Left - Upregulated; Right - Downregulated).

Fig. 2

Gene expression values among the three datasets.

Table 1

Names of the common DEGs from the datasets.

Upregulated: HSPB8, TFAP2B, HLF, SORBS2, EMP1, ABCA8, EDN3, FCER1A, GALNT12, HPGD, GPD1L, ID4, LIMCH1, COL14A1, COBL, ADH1B

Downregulated: BID, PLAU, KRT17///JUP, MMP13, IGF2BP3, IGF2BP2, MMP3SERPINH1, SOX4, HOXA10-HOXA9///MIR196B///HOXA9, CDC25B, GINS1, IFI6, ATAD2, MAGEA6///MAGEA3, UBE2C, LAMP3, PTHLH, HOXC10, MLF1, TOP2A, MAGEA6, CXCL1, CDC6, MFAP2, MIR8071-2///MIR8071-1///IGHV4-31///IGHM///IGHG2///IGHG1, APOBEC3B, SNX10, CXCL8, KIF14, ODC1, CXCL10, LAMB3, DTL, ANO1, LOX, MEST, LAMC2, MAGEA11, HOXA10, MAGEA12, CCL20, SPP1, COL10A1, SHOX2, IGFBP3, TFRC, HOXB7, RBP1, MMP10, CDKN3, MMP1, APOC1, MMP12, HMGA2, ISG15, ECT2, RPL39L, MMP11

GO, and KEGG enrichment analysis: Common shared differentially expressed genes; DAVID analysis showed significant terms in Table 2.

Table 2

GO and KEGG analysis for the DEGs.

Category	Term	P-Value	Genes
GOTERM_BP_DIRECT	GO:0022617∼extracellular matrix disassembly	6.35E-10	MMP12, MMP11, MMP13, LAMB3, MMP1, MMP3, SPP1, LAMC2, MMP10
GOTERM_BP_DIRECT	GO:0030574∼collagen catabolic process	2.01E-07	MMP12, MMP11, MMP13, MMP1, MMP3, COL10A1, MMP10
GOTERM_BP_DIRECT	GO:0008284∼positive regulation of cell proliferation	1.48E-05	CXCL10, TFAP2B, EDN3, ODC1, ID4, KIF14, LAMC2, PTHLH, SOX4, HOXC10, CDC25B
GOTERM_BP_DIRECT	GO:0030198∼extracellular matrix organization	1.29E-04	LAMB3, LOX, COL14A1, MFAP2, SPP1, COL10A1, LAMC2
GOTERM_BP_DIRECT	GO:0001501∼skeletal system development	2.16E-04	HOXA10, SHOX2, COL10A1, PTHLH, SOX4, HOXC10
GOTERM_BP_DIRECT	GO:0032467∼positive regulation of cytokinesis	4.28E-04	KIF14, CDC6, ECT2, CDC25B
GOTERM_BP_DIRECT	GO:0030199∼collagen fibril organization	5.01E-04	MMP11, LOX, COL14A1, SERPINH1
GOTERM_BP_DIRECT	GO:0008285∼negative regulation of cell proliferation	9.90E-04	TFAP2B, CXCL8, IGFBP3, CXCL1, CDC6, PTHLH, SOX4, CDKN3
GOTERM_BP_DIRECT	GO:0032461∼positive regulation of protein oligomerization	0.001384	MMP1, MMP3, BID
GOTERM_BP_DIRECT	GO:0006935∼chemotaxis	0.001428	CXCL10, CXCL8, PLAU, CCL20, CXCL1
GOTERM_BP_DIRECT	GO:0070098∼chemokine-mediated signaling pathway	0.002859	CXCL10, CXCL8, CCL20, CXCL1
GOTERM_BP_DIRECT	GO:0007275∼multicellular organism development	0.004632	HOXA10, MMP11, ANO1, HLF, EDN3, HMGA2, EMP1, HOXB7
GOTERM_BP_DIRECT	GO:0008544∼epidermis development	0.004749	LAMB3, EMP1, LAMC2, PTHLH
GOTERM_BP_DIRECT	GO:0043065∼positive regulation of apoptotic process	0.006892	TOP2A, IGFBP3, HMGA2, ECT2, BID, SOX4
GOTERM_BP_DIRECT	GO:0042769∼DNA damage response, detection of DNA damage	0.009053	HMGA2, DTL, SOX4
GOTERM_BP_DIRECT	GO:0006508∼proteolysis	0.014588	MMP12, MMP11, MMP13, PLAU, MMP1, MMP3, MMP10
GOTERM_BP_DIRECT	GO:0045944∼positive regulation of transcription from RNA polymerase II promoter	0.015501	TOP2A, HOXA10, CXCL10, TFAP2B, HLF, ATAD2, SHOX2, ID4, HMGA2, SOX4
GOTERM_BP_DIRECT	GO:0097070∼ductus arteriosus closure	0.019794	TFAP2B, HPGD
GOTERM_BP_DIRECT	GO:0042035∼regulation of cytokine biosynthetic process	0.027603	IGF2BP3, IGF2BP2
GOTERM_BP_DIRECT	GO:0060326∼cell chemotaxis	0.02776	CXCL10, CCL20, CXCL1
GOTERM_BP_DIRECT	GO:0030593∼neutrophil chemotaxis	0.028555	CXCL8, EDN3, CCL20
GOTERM_BP_DIRECT	GO:0043154∼negative regulation of cysteine-type endopeptidase activity involved in apoptotic process	0.030994	TFAP2B, LAMP3, IFI6
GOTERM_BP_DIRECT	GO:0045444∼fat cell differentiation	0.03437	TFAP2B, ID4, HMGA2
GOTERM_BP_DIRECT	GO:0032330∼regulation of chondrocyte differentiation	0.03535	SHOX2, PTHLH
GOTERM_BP_DIRECT	GO:0030071∼regulation of mitotic metaphase/anaphase transition	0.03535	UBE2C, CDC6
GOTERM_BP_DIRECT	GO:0042127∼regulation of cell proliferation	0.037708	CXCL10, TFRC, PLAU, BID
GOTERM_BP_DIRECT	GO:2000406∼positive regulation of T cell migration	0.039201	CXCL10, CCL20
GOTERM_BP_DIRECT	GO:0001558∼regulation of cell growth	0.040602	TFRC, IGFBP3, KIF14
GOTERM_BP_DIRECT	GO:0009952∼anterior/posterior pattern specification	0.040602	HOXA10, HOXB7, HOXC10
GOTERM_BP_DIRECT	GO:0043085∼positive regulation of catalytic activity	0.041525	IGFBP3, APOC1, CXCL1
GOTERM_BP_DIRECT	GO:0031581∼hemidesmosome assembly	0.046858	LAMB3, LAMC2
GOTERM_BP_DIRECT	GO:0021846∼cell proliferation in forebrain	0.046858	KIF14, HMGA2
GOTERM_CC_DIRECT	GO:0005576∼extracellular region	5.78E-09	CXCL8, LAMB3, EDN3, TFRC, COL14A1, MMP1, CCL20, IGFBP3, MMP3, ISG15, CXCL1, LAMC2, PTHLH, MMP10, MMP12, MMP11, CXCL10, MMP13, LOX, PLAU, MFAP2, APOC1, SPP1, COL10A1
GOTERM_CC_DIRECT	GO:0005578∼proteinaceous extracellular matrix	6.07E-06	MMP12, MMP11, MMP13, LOX, MMP1, COL14A1, MMP3, COL10A1, MMP10
GOTERM_CC_DIRECT	GO:0005615∼extracellular space	2.21E-05	CXCL8, EDN3, TFRC, COL14A1, CCL20, IGFBP3, MMP3, CXCL1, LAMC2, PTHLH, MMP10, CXCL10, MMP13, LOX, PLAU, SERPINH1, SPP1
GOTERM_CC_DIRECT	GO:0005581∼collagen trimer	2.22E-05	MMP13, LOX, MMP1, COL14A1, SERPINH1, COL10A1
GOTERM_CC_DIRECT	GO:0048471∼perinuclear region of cytoplasm	0.007649	TFRC, LAMP3, ODC1, SPP1, COBL, SORBS2, LAMC2, CDKN3
GOTERM_CC_DIRECT	GO:0031012∼extracellular matrix	0.023505	MMP11, MMP13, MMP1, COL14A1, MMP10
GOTERM_CC_DIRECT	GO:0005654∼nucleoplasm	0.042373	GINS1, TOP2A, HPGD, UBE2C, HSPB8, ATAD2, HMGA2, ISG15, CDC6, PTHLH, HOXC10, CDC25B, RBP1, MAGEA11, HOXB7, DTL, SOX4
GOTERM_MF_DIRECT	GO:0004222∼metalloendopeptidase activity	8.53E-05	MMP12, MMP11, MMP13, MMP1, MMP3, MMP10
GOTERM_MF_DIRECT	GO:0004252∼serine-type endopeptidase activity	5.17E-04	MMP12, MMP11, MMP13, PLAU, MMP1, MMP3, MMP10
GOTERM_MF_DIRECT	GO:0008009∼chemokine activity	9.66E-04	CXCL10, CXCL8, CCL20, CXCL1
GOTERM_MF_DIRECT	GO:0001077∼transcriptional activator activity, RNA polymerase II core promoter proximal region sequence-specific binding	0.014362	HOXA10, TFAP2B, HLF, HMGA2, SOX4
GOTERM_MF_DIRECT	GO:0004175∼endopeptidase activity	0.019447	MMP12, MMP1, MMP3
GOTERM_MF_DIRECT	GO:0045182∼translation regulator activity	0.023582	IGF2BP3, IGF2BP2
GOTERM_MF_DIRECT	GO:0005518∼collagen binding	0.023691	MMP13, COL14A1, SERPINH1
GOTERM_MF_DIRECT	GO:0048027∼mRNA 5'-UTR binding	0.035167	IGF2BP3, IGF2BP2
KEGG_PATHWAY	hsa05202: Transcriptional misregulation in cancer	5.57E-06	HOXA10, CXCL8, HPGD, PLAU, IGFBP3, MMP3, HMGA2, MLF1
KEGG_PATHWAY	hsa05323: Rheumatoid arthritis	0.0064	CXCL8, MMP1, CCL20, MMP3
KEGG_PATHWAY	hsa04668: TNF signaling pathway	0.01095	CXCL10, CCL20, MMP3, CXCL1
KEGG_PATHWAY	hsa04622: RIG-I-like receptor signaling pathway	0.036968	CXCL10, CXCL8, ISG15
KEGG_PATHWAY	hsa04062: Chemokine signaling pathway	0.046265	CXCL10, CXCL8, CCL20, CXCL1

PPI network analysis and hub gene screening: 91 nodes and 120 edges were analyzed during the investigation. Genes having more than or equal ten connectivity were considered as hub genes. The selected gene list is given in Table 3. Fig. 3A, Fig. 3BA-3E showed the entire gene network and subsequent clusters.

Table 3

Hub genes list with corresponding connectivity.

Gene Name	Connectivity
CXCL8	16
MMP3	14
MMP1	13
SPP1	13
MMP13	12
UBE2C	11
PLAU	10
CXCL1	10
KIF14	10
TOP2A	10

Fig. 3A

Overall PPI network of the common upregulated and downregulated genes.

Fig. 3B

3D: Four cluster of genes from the overall networks after analysis (3B- Left; 3C- Middle left; 3D- Middle right; 3D- Right).

Survival analysis of the hub genes: Kaplan-Meier survival plot analysis showed the hub genes significantly affect the Overall Survival (OS) of the patient affected with gastric cancer (https://kmplot.com/analysis/). CXCL8, SPP1, MMP13, CXCL1, TOP2A were genes that significantly impacted the survivability of the patients (Fig. 4).

Fig. 4

Overall survival analysis of the hub genes.

Expression plots of the hub genes: Box plots showed all the genes had significant expression levels compared to the standard and cancer tissues, except UBE2C (Fig. 5).

Fig. 5

Boxplots of the hub genes (* mark showed significant difference).

Whereas pathological stage plot analysis of the hub genes showed different expression levels according to their stages (Fig. 6).

Fig. 6

Pathological stage plots of the hub genes.

DNA methylation: The resulting heatmap contains methylation data of 8 transcripts from 324 samples of 450k. In the heatmap, rows represent transcripts, and columns represent samples (green color represents standard profiles, black represents disease profiles) (Fig. 7).

Fig. 7

DNA methylation profile of the hub genes.

lncRNA identification: Three datasets shared only one lncRNA during the screening procedure. Fig. 8 illustrates the lncRNA numbers among the three studies datasets as a Venn diagram. All the names are given in Table 4.

Fig. 8

lncRNA among the three sets of GEO data.

Table 4

lncRNA names of the datasets.

Gene Sets	Common lncRNA	Names
GSE161533GSE20347GSE45670	1	HYMAI
GSE161533GSE20347	1	SNHG17
GSE161533GSE45670	40	TMPO-AS1, PGM5-AS1, BBOX1-AS1, PSMD6-AS2, KCNMB2-AS1, ZFAS1, HCG11, LINC00702, LCAL1, LINC00491, NR2F2-AS1, SLC8A1-AS1, CARMN, ADAMTS9-AS1, LINC01315, MIR100HG, MIAT, PWAR6, LINC00467, PCBP1-AS1, ZNF667-AS1, POU6F2-AS2, CADM3-AS1, LINC00240, FENDRR, TP73-AS1, WDFY3-AS2, LINC01082, LINC00472, MAGI2-AS3, ZNF790-AS1, HLA-F-AS1, ADAMTS9-AS2, ATP2A1-AS1, LINC00938, LINC01140, PDZRN3-AS1, KLF3-AS1, HAGLROS, LINC00844

Drug-genes interactions: Interaction showed that not every DEG has specific drugs to choose from in the current therapeutic practice. According to the findings of this study, very few genes have targeted medication (Fig. 9).

Fig. 9

Interaction between drugs and genes (Oval shapes represent drugs, Hexagonal shapes represent genes).

Discussion

GO enrichment analysis of this study found several terms associated with this neoplasm. This study detects extracellular matrix disassembly as the most significant biological process related to ESCC approved by the previously published article [24]. Collagen catabolic process, positive regulation of cell proliferation, extracellular matrix organizations were also detected as important biological process terms reported by previous authors [10,19,20,27].

The most important cellular component identified was the extracellular region aligned with the findings of the previous studies [29,31]. According to published articles, some other terms on the same functional dysregulation were proteinaceous extracellular matrix, extracellular matrix, collagen trimer, etc. that significantly impacted the ESCC oncogenesis [22,40]. Apart from the metallo endopeptidase activity, several other molecular functions, including chemokine activity, transcriptional activator activity, translation regulatory activity, were included in this study. Previous studies reported similar KEGG pathways linked to the ESCC occurrence [3,13,30].

This study identified five pathways for ESCC pathogenesis. Transcriptional misregulation in cancer was the most significant pathway, a similar route acknowledged by previous authors [6,12,23,39].

PPI networks showed four gene clusters in this study. CXCL8 was the highest corresponding gene which might have a crucial effect on ESCC metastasis. Previous authors detected this gene as the primary cause for ESCC occurrence [5,42]. MMP3, MMP1, and MMP13 gene expression related to cancer occurrence is explained by previous studies [8,38]. TOP2A and UBE2C were also detected by the published articles for cancer progression [21,34].

lncRNA identified in this study was not identified by the previous studies for ESCC occurrence. However, no previous reports about the HYMAI involvements, but all the three datasets shared this lncRNA. This might have an impact on the ESCC progression in a novel way. Several lncRNA identified were shared between 2 datasets, which might have a more significant effect on ESCC.

Drug interactions indicated that all the significant genes targeted therapy were unavailable until now. There needed much more research to precisely target the genes that cause the ESCC.

Ethical approval

This study has no invasive methods involvement regarding ethical approval requirements.

Declaration of competing interest

There is no conflict of interest about the article.