Bioinformatics Analysis of Key Genes and circRNA-miRNA-mRNA Regulatory Network in Gastric Cancer

Gastric cancer (GC) is one of the most common malignancies in the world, with morbidity and mortality ranking second among all cancers. Accumulating evidences indicate that circular RNAs (circRNAs) are closely correlated with tumorigenesis. However, the mechanisms of circRNAs still remain unclear. This study is aimed at determining hub genes and circRNAs and analyzing their potential biological functions in GC. Expression profiles of mRNAs and circRNAs were downloaded from the Gene Expression Omnibus (GEO) data sets of GC and paracancer tissues. Differentially expressed genes (DEGs) and differentially expressed circRNAs (DE-circRNAs) were identified. The target miRNAs of DE-circRNAs and the bidirectional interaction between target miRNAs and DEGs were predicted. Functional analysis was performed, and the protein-protein interaction (PPI) network and the circRNA-miRNA-mRNA network were established. A total of 456 DEGs and 2 DE-circRNAs were identified with 3 mRNA expression profiles and 2 circRNA expression profiles. GO analysis indicated that DEGs were mainly enriched in extracellular matrix and cell adhesion, and KEGG confirmed that DEGs were mainly associated with focal adhesion, the PI3K-Akt signaling pathway, extracellular matrix- (ECM)- receptor interaction, and gastric acid secretion. 15 hub DEGs (BGN, COL1A1, COL1A2, FBN1, FN1, SPARC, SPP1, TIMP1, UBE2C, CCNB1, CD44, CXCL8, COL3A1, COL5A2, and THBS1) were identified from the PPI network. Furthermore, the survival analysis indicate that GC patients with a high expression of the following 9 hub DEGs, namely, BGN, COL1A1, COL1A2, FBN1, FN1, SPARC, SPP1, TIMP1, and UBE2C, had significantly worse overall survival. The circRNA-miRNA-mRNA network was constructed based on 1 circRNA, 15 miRNAs, and 45 DEGs. In addition, the 45 DEGs included 5 hub DEGs. These results suggested that hub DEGs and circRNAs could be implicated in the pathogenesis and development of GC. Our findings provide novel evidence on the circRNA-miRNA-mRNA network and lay the foundation for future research of circRNAs in GC.


Introduction
Gastric cancer (GC) is one of the most common malignancies worldwide, with the morbidity and mortality of GC ranking second among all cancers, especially in East Asia, Eastern Europe, and South America [1]. Although much progress has been made in the diagnosis and treatment of GC, the 5-year overall survival rate remains poor (approximately 20-25%) [2]. Therefore, there is an urgent need to study the mechanism underlying the occurrence and development of GC in order to achieve early diagnosis, effective treatment, and good prognosis for GC.
Bioinformatics analysis, including the use of microarray expression data sets, protein/gene-protein/gene interaction networks, and gene annotation, can be used to study cancer progression and identify potential therapeutic targets for development [3]. In addition, bioinformatics analysis methods can overcome the inconsistent results in the literature as a result of different sample sizes or microarray platforms in individual studies [4]. A large number of studies have used bioinformatics analysis to predict biomarkers for cancer treatments [5][6][7][8].
Circular RNA (circRNA) is a special type of endogenous noncoding RNA formed by reverse splicing of exon or intron cyclization. circRNAs exhibit stability, conservation, abundance, and expression-specific expression [9]. Recent studies have found that circRNAs are associated with many types of cancer, and several circRNAs have been identified as novel cancer biomarkers [10]. circRNAs are members of the competing endogenous RNA (ceRNA) family and act as regulators of miRNA activity; they can competitively inhibit the binding ability of miRNAs and their mRNA targets [11]. Although there are an increasing number of studies on cir-cRNAs, the biological functions of most circRNAs remain unclear. To this end, we constructed a circRNA-miRNA-mRNA network as a means to further evaluate the roles of dysregulated circRNAs and mRNAs in GC.
In this study, we identified 456 differentially expressed genes (DEGs) between 119 human GC tissues and paracancer tissues by analyzing three sets of mRNA expression profiles in a public Gene Expression Omnibus (GEO) data set. Next, functional enrichment analysis and pathway enrichment analysis were performed to explore the roles of these DEGs. Then, a protein-protein interaction (PPI) network was constructed and 15 hub DEGs and three important modules in the network were identified. To assess the prognostic value of these hub DEGs, we also performed the Kaplan-Meier analysis. In addition, two differentially expressed cir-cRNAs (DE-circRNAs) between GC tissues and paracancer tissues were identified by analyzing the two circRNA expression profiles. Finally, a circRNA-miRNA-mRNA network was successfully constructed.
2.3. Gene Ontology Analysis and KEGG Analysis. Gene ontology (GO) analysis is used to provide gene annotation terms [14], while the KEGG database is used for pathway enrichment analysis [15]. clusterProfiler V3.14.0 is an ontologybased R software package for statistical analysis and visualization of functional clusters of genomes or gene clusters [16]. In this study, the GO analysis and pathway enrichment analysis of mRNA were performed using the clusterProfiler package. In addition, a p value < 0.05 and a q value < 0.05 were chosen as the cutoff criteria for significant pathway terms.

PPI Network and circRNA-miRNA-mRNA Network
Construction. The STRING database is a search tool for the retrieval of interacting genes or proteins (https://string-db .org), which can then be used to establish a PPI network [17]. 456 DEGs were imported into the STRING database; an interaction score > 0:4 [18] was used as the extraction cutoff standard for the PPI pair. Then, Cytoscape_3.7.1 (https:// cytoscape.org) was used to visualize the PPI network [19]. The cytoHubba plug-in can be used to screen hub DEGs with the node degree. The MCODE Plug-in was used to filter important modules in the PPI network with a degree cutoff ≥ 2, node score cutoff = 0:2, K − core ≥ 2, and max:depth = 100 as the cutoff criteria [18]. In addition, the obtained circRNA-miRNA pairs and miRNA-mRNA pairs were combined to construct a circRNA-miRNA-mRNA network. The network was also visualized with Cytoscape.
2.5. Verification of Gene Expression. UALCAN (http://ualcan .path.uab.edu/) is a website for online analysis and mining of TCGA databases; it includes 34 normal stomach tissues and 415 GC tissues [20]. In this study, this database was used to verify the expression of hub DEGs in the PPI network. For this, p < 0:05 was considered statistically significant.
2.6. Survival Analysis. The Kaplan-Meier plotter (https:// kmplot.com/analysis/) is an online database that evaluates the prognostic value of biomarkers in breast cancer, ovarian cancer, lung cancer, GC, and other cancers [21]. In this study, the Kaplan-Meier plotter database was used to assess the overall survival of hub DEGs among 1065 samples of GC. The hazard ratio (HR) and corresponding 95% confidence intervals were calculated. p < 0:05 was considered statistically significant.

GO and KEGG Pathway Enrichment Analyses of DEGs.
The functions of the DEGs included three GO categories: molecular function (MF), biological process (BP), and cellular component (CC). The top 10 enriched GO terms showed that the DEGs were involved in an extracellular matrix, extracellular structure organization, collagen-containing extracellular matrix, extracellular matrix organization, extracellular matrix structural constituent, cell-substrate adhesion, skeletal system development, endoplasmic reticulum lumen, apical part of the cell, and regulation of body fluid levels ( Figure 2(a)). The KEGG analysis confirmed that the DEGs were mainly associated with protein digestion and absorption, focal adhesion, the PI3K-Akt signaling pathway, extracellular matrix-(ECM-) receptor interaction, gastric acid secretion, and so on (Figure 2(b)).

PPI Network and Modular
Analysis. 456 DEGs were imported into the STRING database to explore the interrelationships between the various genes. 456 DEGs were used to establish the PPI network using the Cytoscape software. The PPI network consisted of 369 nodes and 1570 edges (Figure 3(a)). The cytoHubba plug-in in Cytoscape was used to screen the top 15 genes with node degree indicating the hub DEGs from the PPI network, including BGN, COL1A1, COL1A2, FBN1, FN1, SPARC, SPP1, TIMP1, UBE2C, CCNB1, CD44, CXCL8, COL3A1, COL5A2, and THBS1 ( Figure 3(b) and Table 1). All were upregulated DEGs in the microarray expression profiles, and the mRNA expression levels of these genes were also significantly upregulated in GC tissues compared with normal stomach samples in     3.6. Construction of the circRNA-miRNA-mRNA Network. To better understand the role of DE-circRNAs in GC, a circRNA-miRNA-mRNA network was established. Using the Circular RNA Interactome online tool, there were 108 target miRNAs for hsa_circ_0001013 and 11 target miRNAs for hsa_circ_0021087. We selected the bidirectional search function in the mirDIP database and then entered 167 upregulated DEGs and 108 miRNAs that were targeted by upregulated hsa_circ_0001013. Ten classic databases were selected to predict the bidirectional relationship between miRNAs and upregulated DEGs. The same method was used to predict the targeting relationship between 289 downregulated DEGs and 11 target miRNAs of downregulated hsa_circ_ 0021087. The results showed that 15 miRNAs targeted 45 upregulated DEGs and also interacted with hsa_circ_ 0001013 (Supplementary Table 3

GO and KEGG Pathway Enrichment Analyses in
circRNA-miRNA-mRNA Network. To date, the functions of hsa_circ_0001013 have not been reported. Therefore, a functional forecast of hsa_circ_0001013 was performed according to the mRNA annotations in the hsa_circ_0001013-related ceRNA network. GO analysis indicated that the function prediction of the circRNA-miRNA-mRNA network was mainly related to ECM and adhesion (Figure 7(a)).    BioMed Research International analysis showed that the network was mainly enriched by the following terms: the PI3K-Akt signaling pathway, ECMreceptor interaction, focal adhesion, and human papillomavirus infection (Figure 7(b)).

Discussion
The tumorigenesis and development of GC is a relatively complex process because of the involvement of aberrations in gene expression regulatory networks. Traditionally, gene expression regulation analyses mainly focused on protein-coding genes (mRNAs), until the discovery of diseaserelated noncoding RNAs, including circRNAs. Still, little is known about the regulatory mechanisms of circRNAs in GC. To address this gap in understanding, this study used microarray data and performed bioinformatics analysis to identify disease-associated hub DEGs and the circRNA-miRNA-mRNA network. This study obtained 456 DEGs based on three microarray profiles, including 167 upregulated genes and 289 downregulated genes. GO analysis showed that the 456 DEGs were mainly related to ECM, extracellular tissue, and collagen-    9 BioMed Research International containing ECM. Similarly, the results of KEGG analysis showed that the upregulated DEGs were not only enriched in protein digestion and absorption and the PI3K-Akt signaling pathway, but were also involved in focal adhesion and ECM-receptor interaction. The most significant pathway for downregulated DEGs was the metabolism of xenobiotics by cytochrome P450 (CYP3A5, CYP2C9, CYP2S1, ADH1A, ADH1C, ADH7, AKR1C1, AKR7A3, SULT2A1, UGT2B15, ALDH3A1, and GSTA1). Cytochrome P450 (CYP) is present in extrahepatic tissues and plays a key role in target tissue metabolic activation of xenobiotic compounds, the detoxification of toxic compounds following CYP-catalyzed biotransformation, and activation of inert xenobiotics, including drugs, to become toxicants [24]. A recent study found that low expression of ADH1A, ADH1C, and ADH7 was significantly associated with the increased risk of mortality for GC patients receiving 5-fluorouracil-(5-FU-) based adjuvant chemotherapy [25]. However, how the deficient CYP enzymes contribute to oncogenicity in the development of GC is not yet clear. It is worth noting that both GO analysis and KEGG analysis revealed relationships with ECM. ECM is a key regulator of cell behavior and phenotype; it activates or suppresses distinct sets of intracellular signaling pathways to modulate downstream cellular decisions, such as focal adhesions and collagen remodeling [26]. In many different types of cancer, ECM is highly dysregulated and the loss of tissue ECM homeostasis and integrity is seen as a hallmark of cancer and the typical transitional events in progression and metastasis [27]. Our data showed that GC progression was powerfully influenced by ECM and ECMrelated events. Thus, remodeling of ECM in GC patients should be considered a potential therapeutic target.
To explore the molecular mechanism of GC, a GCrelated PPI network was constructed and 15 hub DEGs were identified. These hub DEGs were all overexpressed in GC tissues. The Kaplan-Meier plotter was used to evaluate the effects of the 15 hub genes on the survival of GC patients. The results showed that overexpression of BGN, COL1A1, COL1A2, FBN1, FN1, SPARC, SPP1, TIMP1, and UBE2C led to poor prognosis in GC. As the main components of ECM, SPARC (acid-secreted and cysteine-rich secreted protein), BGN (biglycan), FBN1 (fibrillin 1), SPP1 (secreted phosphoprotein 1), fibronectin (FN1), and collagen play important roles in cell proliferation, differentiation, migration, and metastasis in GC [28]. Typically, two COL1A1 chains pair with one COL1A2 chain to form a triple helix of type I collagen [29]. Studies have shown that COL1A1 and COL1A2 are more highly expressed in GC tissues than in normal tissues, and both are related to the invasion and proliferation of GC cells [30,31]. Type I collagen is required for collagen in bones, but it is also involved in ECM synthesis and in the promotion of changes to cell shape. The SPARC protein product has been associated with a variety of cancers, including GC [32]; it may affect GC metastasis by regulating the tumor microenvironment [33]. FBN1 is an extracellular matrix glycoprotein that can be used as a structural   component of calcium-bound microfibers. The expression of FBN1 is upregulated in GC tissues and cells; the silence of FBN1 can inhibit the proliferative, migratory, and invasive abilities of GC cells [34]. As a high molecular weight glycoprotein, recent evidence has shown that FN1 is associated with a variety of cancers, and it is also involved in the invasion and migration of GC [35]. BGN (biglycan) is a member of the leucine-rich small proteoglycan family; it can enhance the invasion, migration, and formation of endometrial cancer cells [36], and it can promote the invasion of GC cells by activating the FAK signaling pathway [37]. The protein encoded by SPP1, also known as osteopontin, facilitates the combination of mineralized bone matrix and osteoclasts and is overexpressed in a variety of cancers, including breast cancer, lung cancer, pancreatic cancer, and GC [38]. Overexpression of SPP1 leads to poor GC prognosis, and positive osteopontin has an important role in the prediction of blood-borne metastasis [39]. UBE2C encodes a member of the E2 ubiquitin-binding enzyme family, which is involved in the ubiquitin-proteasome system and plays important roles in cell mitotic withdrawal and cell cycle progression. Knocking down UBE2C can inhibit the occurrence and development of GC through the Wnt/β-catenin and PI3K/Akt signaling pathways [40]. UBE2C is also related to poor GC prognosis [41]. TIMP1, a member of the TIMP gene family, is a natural inhibitor of matrix metalloproteinase. TIMP1 plays important roles in cell proliferation, tumorigenesis, angiogenesis, and antiapoptosis [42]. It has been reported that TIMP1 can be overexpressed in human GC by relying on the NF-κB pathway and has the ability to regulate the proliferation of GC cells [43].
The current results revealed that overexpression of CCNB1, CD44, and CXCL8 was associated with significantly improved prognosis of GC patients. Studies have confirmed that CD44 is a nonkinase transmembrane glycoprotein that can promote tumor cell proliferation and invasion [44]. CD44 is also the most specific biomarker for the detection and isolation of oncogenic and chemoresistive cancer stem cells in noncardiac GC [45]. CXCL8 is an important member of the CXC chemokine family and plays an important role in the proliferation, migration, and activation of inflammatory systems in tumor cells [46]. CCNB1 (encoding cyclin B1) plays a vital role in cell mitosis. Overexpression of CCNB1 contributes to the proliferation of colorectal cancer cells [47], promotes the proliferation, migration, and invasion of bladder cancer [48], and also leads to poor prognosis of hepatocellular carcinoma [49]. CCNB1 can promote the proliferation and metastasis of GC by participating in the

13
BioMed Research International composition of the heterogeneous ribonucleoprotein-CCNB1/CENPF axis [50]. The overexpression of CCNB1, CD44, and CXCL8 in this study significantly improved the quality of life of GC patients. Therefore, the biological significance of these three genes in GC needs to be determined in further research. Taken together, these findings indicate that these nine hub DEGs are closely related to the prognosis of GC patients. As such, they may be potential prognostic indicators in GC.
Recent studies have found that circRNAs are associated with many types of cancer, and circRNAs have been proposed as novel cancer biomarkers [10]. The biological function of circRNAs is a "miRNA sponge"; circRNAs can competitively inhibit the binding ability of miRNAs and their mRNA targets [11]. There are an increasing number of studies on the role of circRNAs in GC, but the biological functions of most circRNAs remain unclear. In this study, we identified one upregulated DE-circRNA (hsa_circ_0001013) and constructed a circRNA-miRNA-mRNA network to estimate the function of hsa_circ_0001013 in GC. The results showed that hsa_circ_0001013 might play pivotal regulating roles in the gene expression of ECM and adhesion, and it also appears to be involved in the PI3K-Akt signaling pathway and ECM-receptor interaction in GC. Although there is currently no relevant research on the function of hsa_circ_ 0001013, the network revealed some important information. For instance, it was found that miR-182-5p, as a suppressed miRNA, improves migration and invasion of GC [51]. In our data, hsa_circ_0001013 was found to regulate the expression of FBN1, FN1, THBS1, and so on, through competing miRNA response elements (MREs) of miRNA-182-5p, which results in adhesion and metastasis in GC. In addition to miR-182-5p, studies have also shown that miR-758-3p and miR-145-5p are also downregulated in GC tissues and play important roles in promoting the migration and invasion of GC cells [52,53]. In our circRNA-miRNA-mRNA network, the hsa_circ_0001013/miR-145-5p/FN1 axis and the hsa_circ_ 0001013/miR-758-3p/FBN1 axis may also affect the development of GC. These data indicate that circRNA plays a key role in GC. It is worth noting that miR-1197, miR-323-3p, miR-507, and miR-330-5p have not been studied in GC, but studies have shown that they play key roles in other tumors. For example, downregulated miR-1197 can inhibit the progression of human non-small-cell lung cancer by upregulating HOXC11 [54]. Further, miR-323-3p can inhibit the apoptosis of human lung cancer cells and inhibit the invasion and metastasis of pancreatic ductal adenocarcinoma cells [55,56]. Moreover, miR-507 participates in the hsa_ circ_0005394/miR-507/E2F3 axis and affects the process of hepatocellular carcinoma [57]. miR-330-5p can affect the expression of ELK1 and thus affect the proliferation, migration, and invasion of colon cancer cells [58]. It can also affect the expression of KLK4 to change the development process of ovarian cancer [59]. The relationships in the circRNA-miRNA-mRNA network require further exploration. Nonetheless, this pioneering network might provide a novel understanding of GC. Further research on the associated functions of hsa_circ_0001013 is being carried out in our laboratory.

Conclusions
In summary, through bioinformatics analysis, we identified 15 hub DEGs, of which, nine hub DEGs were associated with poor prognosis of GC patients, including BGN, COL1A1, COL1A2, FBN1, FN1, SPARC, SPP1, TIMP1, and UBE2C. In addition, two DE-circRNAs were identified. In order to identify the regulatory mechanism of circRNAs in GC, we constructed a related circRNA-miRNA-mRNA regulatory network of hsa_circ_0001013. Although these results need to be further verified, these hub DEGs and previously unreported hsa_circ_0001013 may play crucial roles in GC and may provide new ideas for the diagnosis, prognosis, and therapeutic targeting of GC.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.